Converting latitude and longitude points to UTM

Question

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

Answer 1

In your question you are not clear whether you already read in your data set into a data.frame or matrix. So I assume in the following you have your data set in a text file:

# read in data
dataset = read.table("dataset.txt", header=T)

# ... or use example data
dataset = read.table(text="ID X Y
1 118 10
2 119 50
3 100 12
4 101 12", header=T, sep=" ")

# create a matrix from columns X & Y and use project as in the question
project(as.matrix(dataset[,c("X","Y")]), "+proj=utm +zone=51 ellps=WGS84")
#             [,1]    [,2]
# [1,]   -48636.65 1109577
# [2,]   213372.05 5546301
# ...

Update:

The comments suggest that the problem comes from applying project() to data.frame. project() does not work on data.frames since it checks for is.numeric(). Therefore, you need to convert data to a matrix as in my example above. If you want to stick to your code that uses cbind() you have to do the following:

 X <- dd[,"X"]
 Y <- dd[,"Y"]
 xy <- cbind(X,Y) 

The difference between dd["X"] and dd[,"X"] is that latter will not return a data.frame and as a consequence cbind() will yield a matrix and not a data.frame.

Answer 2

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 

Answer 3

Regarding the example, the default UTM Zone for the given coordinates is 50. It is not advised to project coordinates into far away zones. You can check the conversion using the NCAT tool from NOAA.

The code below uses the sf package for doing the conversion.

library(sf)
library(tidyverse)

# Coordinate examples with expected UTM values
coord_sample <- data.frame(
  "Northing" = c(1105578.589, 5540547.370),
  "Easting" = c(609600.773, 643329.124),
  "Latitude" = c(10, 50),
  "Longitude" = c(118, 119))

#' Find UTM EPSG code from Latitude and Longitude coordinates (EPSG 4326 WGS84)
#' (vectorised)
#' Source: https://geocompr.robinlovelace.net/reproj-geo-data.html
#' Source: https://gis.stackexchange.com/questions/13291/computing-utm-zone-from-lat-long-point
LatLonToUTMEPSGCode <- function(lat, lon) {

  zone_number <- (floor((lon + 180) / 6) %% 60) + 1

  # Special zones for Norway
  cond_32 <- lat >= 56.0 & lat < 64.0 & lon >= 3.0 & lon < 12.0
  zone_number[cond_32] <- 32

  # Special zones for Svalbard
  cond_lat <- lat >= 72.0 & lat < 84.0

  cond_31 <- cond_lat & lon >= 0.0 & lon <  9.0
  zone_number[cond_31] <- 31

  cond_33 <- cond_lat & lon >= 9.0 & lon < 21.0
  zone_number[cond_33] <- 33

  cond_35 <- cond_lat & lon >= 21.0 & lon < 33.0
  zone_number[cond_35] <- 35

  cond_37 <- cond_lat & lon >= 33.0 & lon < 42.0
  zone_number[cond_37] <- 37

  # EPSG code
  utm <- zone_number
  utm[lat > 0] <- utm[lat > 0] + 32600
  utm[lat <= 0] <- utm[lat <= 0] + 32700

  return(utm)
}

sf_sample <- sf::st_as_sf(coord_sample, coords = c("Longitude", "Latitude"),
                          crs = 4326)

sf_sample %>%
  do(cbind(., st_coordinates(.))) %>%
  mutate(EPSG = LatLonToUTMEPSGCode(lat = Y, lon = X)) %>%
  group_by(EPSG) %>%
  do(cbind(as.data.frame(.) %>% select(Northing, Easting),
           st_coordinates(st_transform(., crs = head(.$EPSG, 1))))) %>%
  ungroup()

You can check the conversion by comparing with the expected values:

# A tibble: 2 x 5
   EPSG Northing Easting      X       Y
  <dbl>    <dbl>   <dbl>  <dbl>   <dbl>
1 32650  1105579  609601 609601 1105579
2 32650  5540547  643329 643329 5540547

Answer 4

Based on the code above, I also added a version of zone and hemisphere detection (that solves conversion problems, as described in some comments) + shorthand for CRS string and conversion back to WSG86:

library(dplyr)
library(sp)
library(rgdal)
library(tibble)

find_UTM_zone <- function(longitude, latitude) {

 # Special zones for Svalbard and Norway
    if (latitude >= 72.0 && latitude < 84.0 ) 
        if (longitude >= 0.0  && longitude <  9.0) 
              return(31);
    if (longitude >= 9.0  && longitude < 21.0)
          return(33)
    if (longitude >= 21.0 && longitude < 33.0)
          return(35)
    if (longitude >= 33.0 && longitude < 42.0) 
          return(37)

    (floor((longitude + 180) / 6) %% 60) + 1
}


find_UTM_hemisphere <- function(latitude) {

    ifelse(latitude > 0, "north", "south")
}

# returns a DF containing the UTM values, the zone and the hemisphere
longlat_to_UTM <- function(long, lat, units = 'm') {

    df <- data.frame(
        id = seq_along(long), 
        x = long, 
        y = lat
    )
    sp::coordinates(df) <- c("x", "y")

    hemisphere <- find_UTM_hemisphere(lat)
    zone <- find_UTM_zone(long, lat)

    sp::proj4string(df) <- sp::CRS("+init=epsg:4326") 
    CRSstring <- paste0(
        "+proj=utm +zone=", zone,
        " +ellps=WGS84",
        " +", hemisphere,
        " +units=", units)
    if (dplyr::n_distinct(CRSstring) > 1L) 
        stop("multiple zone/hemisphere detected")

    res <- sp::spTransform(df, sp::CRS(CRSstring[1L])) %>%
        tibble::as_data_frame() %>%
        dplyr::mutate(
            zone = zone,
            hemisphere = hemisphere
        )

    res
}

UTM_to_longlat <- function(utm_df, zone, hemisphere) {

    CRSstring <- paste0("+proj=utm +zone=", zone, " +", hemisphere)
    utmcoor <- sp::SpatialPoints(utm_df, proj4string = sp::CRS(CRSstring))
    longlatcoor <- sp::spTransform(utmcoor, sp::CRS("+init=epsg:4326"))
    tibble::as_data_frame(longlatcoor)
}

Where CRS("+init=epsg:4326") is shorthand for CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0").

Finding UTM zone reference: http://www.igorexchange.com/node/927

Answer 5

Converting latitude and longitude points to UTM

library(sp)
library(rgdal)

#Function
LongLatToUTM<-function(x,y,zone){
 xy <- data.frame(ID = 1:length(x), X = x, Y = y)
 coordinates(xy) <- c("X", "Y")
 proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")  ## for example
 res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
 return(as.data.frame(res))
}

# Example
x<-c( -94.99729,-94.99726,-94.99457,-94.99458,-94.99729)
y<-c( 29.17112, 29.17107, 29.17273, 29.17278, 29.17112)
LongLatToUTM(x,y,15)