Search code examples
rrgdalr-sp

How to transform Longitude/Latitude degree to UTM? R. rgdal, sp


My data is in latitude and longitude projection in the form of degrees, minutes and seconds. This has been and stored in text form (character-class, Ignore the row.names).

> pasaporte
         Latitud        Longitud
4  13°50¨52" sur 73°45¨12" oeste
36    13°01¨ sur    75°05¨ oeste
46 13°09¨26" sur 74°13¨22" oeste

By means of the answer provided in the mailing list, I've transformed the data to decimal form...

> pasaporte
    Latitud Longitud
4  13.84778 73.75333
36 13.01667 75.08333
46 13.15722 74.22278

After that, I've transformed the data.frame to a SpatialPoints object.

xy <- data.frame(cbind(round(pasaporte[, 'Longitud'], 5), round(pasaporte[, 'Latitud'], 5))) #Rounded the decimals out of doubt of it interfering later
xy <- SpatialPoints(coords = xy,
                proj4string =  CRS('+proj=longlat +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs'))

Then proceeded to transform the projection via spTransform

xy_utm <- spTransform(xy, CRS('+proj=utm +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84:0,0,0'))

I've converted the first coordinate in the example, using this website (make shure to provide the negative since the coordinates refer to South and West), and provides the following coordinates: 633726.46 mE 8468757.51 mN Zone 18L, which is correct. In contrast the sp object after the transformation is:

SpatialPoints:
           X1       X2
[1,] 12051333 14804033
[2,] 12569894 14792671
[3,] 12301330 14684870
Coordinate Reference System (CRS) arguments: +proj=utm +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0 

The spTransform manual specifically mentions that the metadata should be correctly provided. My data uses the WGS84 ellipsoid and lies in the 18L zone. To be more precise the epsg:32718. Where am I'm going wrong?


Solution

  • I would recommend to "project" (i.e. multiply by -1) the points in the -180.0000, -90.0000, 180.0000, 90.0000 bounds, as they are lying in the southern and western hemispheres.

    library(sp)
    
    pasaporte <- structure(list(Latitud = c(13.84778, 13.01667, 13.15722), 
                                Longitud = c(73.75333, 75.08333, 74.22278)),
                           .Names = c("Latitud", "Longitud"), class = "data.frame", 
                           row.names = c("4", "36", "46"))
    
    pasaporte <- pasaporte * -1
    
    pasaporte
    #      Latitud  Longitud
    # 4  -13.84778 -73.75333
    # 36 -13.01667 -75.08333
    # 46 -13.15722 -74.22278
    
    # Conversion into SpatialPoints
    coordinates(pasaporte) <- ~Longitud+Latitud
    # Setting default projection
    proj4string(pasaporte) <- CRS('+init=epsg:4326')
    
    pasaporte
    # SpatialPoints:
    #     Longitud   Latitud
    # 4  -73.75333 -13.84778
    # 36 -75.08333 -13.01667
    # 46 -74.22278 -13.15722
    # Coordinate Reference System (CRS) arguments: +init=epsg:4326
    # +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
    
    # Projection
    # Be sure to have rgdal first installed.
    xy_utm <- spTransform(pasaporte, CRS('+init=epsg:32718'))
    xy_utm
    # SpatialPoints:
    #    Longitud Latitud
    # 4  634726.5 8468758
    # 36 490964.2 8561019
    # 46 584231.8 8545348
    # Coordinate Reference System (CRS) arguments: +init=epsg:32718 +proj=utm
    # +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84
    # +towgs84=0,0,0