[GIS] Converting decimal degrees (longitude & latitude) to UTM WGS84 in R? (Southern hemisphere)

coordinate systemcoordinatesprojrrgdal

I'm trying to convert some longitude and latitude coordinates to UTM WGS84 projection. I'm using rgdal, although I can accommodate to other package (in R). I read the data and converted it to a SpatialPoints object:

coordinates(res) <- ~Lat_dec + Lon_dec
proj4string(res) <- CRS('+proj=longlat +datum=WGS84')

my points:

res
SpatialPoints:
   Lon_dec   Lat_dec
35 -78.60972 -9.073333
36 -78.61806 -9.078611
37 -78.60889 -9.093611
38 -78.59167 -9.088889
39 -78.57917 -9.100833
40 -78.56833 -9.112778
41 -78.59278 -9.113333
42 -78.61139 -9.119167
43 -78.57778 -9.122778
44 -78.56889 -9.135833
45 -78.59278 -9.137222
46 -78.60222 -9.160278
47 -78.57472 -9.152778
48 -78.59583 -9.152778
Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

The following error is no longer produced since swapping the coordinates to negative values (south america)
Now, when I try to transform the points using spTransform(res, CRS('+proj=utm +zone=17 ellps=WGS84')) I get an error and a warning:

non finite transformation detected:
  Lat_dec   Lon_dec                     
  9.073333 78.609722       Inf       Inf 
Error in spTransform(res, CRS("+proj=utm +zone=17 ellps=WGS84")) : 
failure in points 1
 In addition: Warning message:
 In spTransform(res, CRS("+proj=utm +zone=17 ellps=WGS84")) :
   14 projected point(s) not finite

When I try to convert using project(coordinates(res), '+proj=utm +zone=17 ellps=WGS84') I get the projection:

   [,1] [,2]
 [1,]  Inf  Inf
 [2,]  Inf  Inf
 [3,]  Inf  Inf
 [4,]  Inf  Inf
 [5,]  Inf  Inf
 [6,]  Inf  Inf
 [7,]  Inf  Inf
 [8,]  Inf  Inf
 [9,]  Inf  Inf
[10,]  Inf  Inf
[11,]  Inf  Inf
[12,]  Inf  Inf
[13,]  Inf  Inf
[14,]  Inf  Inf
Warning message:
In project(coordinates(res), "+proj=utm +zone=17 ellps=WGS84") :
  14 projected point(s) not finite

Changed according to the suggestions below and now this is produced

> project(coordinates(res), '+proj=utm +zone=17 ellps=WGS84')
          [,1]     [,2]
 [1,] 762744.8 -1003824
 [2,] 761824.4 -1004402
 [3,] 762821.6 -1006069
 [4,] 764719.2 -1005559
 [5,] 766085.1 -1006890
 [6,] 767267.6 -1008219
 [7,] 764579.0 -1008263
 [8,] 762528.0 -1008895
 [9,] 766221.6 -1009319
[10,] 767189.3 -1010770
[11,] 764561.3 -1010906
[12,] 763505.8 -1013450
[13,] 766535.2 -1012641
[14,] 764213.9 -1012625

Which are undeniably wrong. The latitude should be positive and near 899####.

Best Answer

Add a +south parameter to the utm PROJ.4 string. That will automatically add a false northing value of 10000000 meters.

"+proj=utm +zone=17 +south=T ellps=WGS84"

Note: I can't find the "=T" in the PROJ.4 documentation on the "+south" parameter, but two people agree that it works when using R.

Related Solutions

[GIS] Slight translation with coordinate-transformation from Gauss-Kruger to WGS 84. Possible Reasons

Your Gauss-Krueger projection uses +datum=potsdam. Up to 2012, this was hard coded in proj4 to a very unprecise value using a 3-parameter-transformation.

You find more exact values for 7-parameter transformations in this topic:

http://forum.openstreetmap.org/viewtopic.php?id=12723

There is an even better ntv2-grid transformation available here (take the binary), that has to be in the same folder as your application and data, unlsss you specify full pathnames.

To compare the different possible values, I made this test batch file:

echo epsg31467-epsg4326 >out.txt
cs2cs +init=epsg:31467 +to +init=epsg:4326 31467.txt >>out.txt

echo proj-Definition epsg >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs  +to +init=epsg:4326 31467.txt >>out.txt

echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:4326 31467.txt >>out.txt

echo proj-definition nadgrid >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +x_0=3500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=./BETA2007.gsb +wktext +to +init=epsg:4326 31467.txt >>out.txt

echo epsg31467-epsg3785 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:3785 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:3785 31467.txt >>out.txt

echo proj-definition nadgrid >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +x_0=3500000 +y_0=0 +k=1.000000 +ellps=bessel +units=m +nadgrids=./BETA2007.gsb +wktext +to +init=epsg:3785 31467.txt >>out.txt

echo epsg31467-epsg3857 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:3857 31467.txt >>out.txt
echo proj-definition Qgis >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +init=epsg:3857 31467.txt >>out.txt

echo epsg31467-epsg900913 >>out.txt
cs2cs +init=epsg:31467 +to +init=epsg:900913 31467.txt >>out.txt
echo epsg31467-proj900913 >>out.txt
cs2cs +init=epsg:31467 +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +over +no_defs 31467.txt >>out.txt
echo proj31467-proj900913 >>out.txt
cs2cs +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel +towgs84=582,105,414,1.04,0.35,-3.08,8.3 +units=m +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +over +no_defs 31467.txt >>out.txt

with any sample Gauss-Krüger lon-lat coordinate pair in 31467.txt

[GIS] Problem points into a shape polygon

You have to transform your polygon (shapefile) Coordinate Reference System (CRS) with the function spTransform from the sp package. You can simply do this:

polygon.t <- spTransform(polygon, CRSobj = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

I also give you a complete reproducible example below:

Get example data:

# Load libraries
library('rgdal')
library('raster')

# Get SpatialPolygonsDataFrame object example
polygon <- getData('GADM', country = 'URY', level = 0)

# Plot data
plot(polygon)

# Make sample points 
sample.spdf <- spsample(polygon, n = 200, type = "random")
sample.spdf <- spsample(sample.spdf, n = 200, type = "random") # get points outside polygon
sample.spdf <- SpatialPointsDataFrame(sample.spdf, data = data.frame("ID" = 1:200, "V2" = sample(x = 1:1000, size = 200))) # add random data

Plot data:

# Plot points and polygon
plot(polygon)
plot(sample.spdf, add = TRUE, pch = 19, cex = 1)
box()

# Interactive map with Mapview
library(mapview)
mapview(polygon) + mapview(sample.spdf)

Spatial Query: points in polygon with different CRS:

# Apply transformation to sample.spdf
sample.spdf.t <- spTransform(sample.spdf, CRSobj = "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

# Ask layers CRS
projection(sample.spdf.t) # [1] "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
projection(polygon) # [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

pt.in.poly.sp <- sp::over(sample.spdf.t, polygon, fn = NULL) 
# Error: identicalCRS(x, y) is not TRUE
pt.in.poly.rgeos <- rgeos::gIntersection(sample.spdf.t, polygon) 
# spgeom1 and spgeom2 have different proj4 strings

# Plot points transformed + polygon
plot(polygon)
plot(sampleSP.t, add = TRUE, pch = 19, cex = 0.5)

Spatial Query: points in polygon with same CRS:

# Ask layers CRS
projection(sample.spdf) # [1] [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
projection(polygon) # [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

pt.in.poly.sp <- sp::over(sample.spdf, polygon, fn = NULL)
pt.in.poly.rgeos <- rgeos::gIntersection(sample.spdf, polygon)

# Plot points in polygon
plot(polygon)
plot(sample.spdf, pch = 19, cex = 1, add = TRUE)
plot(pt.in.poly.rgeos, pch = 19, cex = 0.5, col = 'green', add = TRUE)
box()

# Interactive map with Mapview
mapview(polygon) + mapview(sample.spdf) + mapview(pt.in.poly.rgeos)