Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Converting latitude and longitude points to UTM

I found a fairly simple example of how to do this but I cant get it to work for me. I'm pretty new to R

library(rgdal) 
xy <- cbind(c(118, 119), c(10, 50)) 
project(xy, "+proj=utm +zone=51 ellps=WGS84") 
          [,1]    [,2] 
[1,] -48636.65 1109577 
[2,] 213372.05 5546301

But this is with example numbers. I have thousands of coordinates I have to transform and I cant figure out how to get them from my table to into this script

My data set has 3 columns, ID, X, and Y. How can I transform them using this equation? I've been stuck on this for weeks

like image 681
Colin Avatar asked Sep 05 '13 15:09

Colin


1 Answers

To ensure that appropriate projection metadata are at every step associated with the coordinates, I'd suggest converting the points to a SpatialPointsDataFrame object as soon as possible.

See ?"SpatialPointsDataFrame-class" for more on how to convert simple data.frames or matrices to SpatialPointsDataFrame objects.

library(sp)
library(rgdal)

xy <- data.frame(ID = 1:2, X = c(118, 119), Y = c(10, 50))
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  ## for example

res <- spTransform(xy, CRS("+proj=utm +zone=51 ellps=WGS84"))
res
#            coordinates ID
# 1 (-48636.65, 1109577)  1
# 2    (213372, 5546301)  2

## For a SpatialPoints object rather than a SpatialPointsDataFrame, just do: 
as(res, "SpatialPoints")
# SpatialPoints:
#              x       y
# [1,] -48636.65 1109577
# [2,] 213372.05 5546301
# Coordinate Reference System (CRS) arguments: +proj=utm +zone=51
# +ellps=WGS84 
like image 65
Josh O'Brien Avatar answered Oct 01 '22 17:10

Josh O'Brien