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

They

It's going to be hard to help you if you don't give us some of the numbers that don't work (plus possibly some description of the format in which they're stored). Also, do you know that all of your lat-long coordinates fall with a single UTM zone? – Hephzibah

Well its not so much that the numbers don't work because Im not putting in numbers. I need the script to read it from a table of thousands of coordinates. I just dont know how to do that. dd <- read.csv(file.choose(), header = T) X <- dd["X"] Y <- dd["Y"] xy <- cbind(c(X), c(Y)) project(xy, "+proj=utm +zone=51 ellps=WGS84") Not even sure if that equation makes sense but it doesnt run. "xy not numeric" Also I dont know if all my coordinates fall within a single UTM zone. I hate sounding like an idiot but this is new to me but I dont know what that means – They

What kind of table? (Is it stored in a text file? A csv? A relational database? An Excel file? etc. etc. etc.) It sounds like you're really at this point asking how to read data into R. In which case, try using Google or searching in the SO search bar to the above right and then ask again if you can't figure it out. – Hephzibah

I no how to read in data. Its a csv file. The example above is example numbers. My question is after i read in X (longitude) and Y (latitude) from my table. How does that data fit into the example script – They

Answer 1

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 

Candelaria

What is the variable "SP" in the spTransform() call? If I change that to "xy", I don't get the same output as shown below that call. – Moulmein

@Moulmein -- Oops, my bad. That SP was left over from my initial answer, since edited, which had involved a SpatialPoints object rather than a SpatialPointsDataFrame. For your future reference, you can view earlier versions of SO answers by clicking on the "edited ... ago" link directly below an edited answer. Thanks for catching that. – Hephzibah

Answer 2

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)

Cutcliffe

Answer 3

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.

Doubtful

Thank you. This does seem to have produced something which is a relief! Is the output in meters? Most of my coordinates now look like -9596387 -5670257, -9597204 -5666367, etc. Does that seem right? Some people have mentioned about UTM zones. From google i'd suggest appears to be spread across 39 l, 39 K and 38 k if that makes sense. Does that have to be taken into account in the script. As i said I just copied an example i found online. Im sorry if i sound like an idiot, im just trying to get a grasp on this. Thanks again – They

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

Stanislaus

Answer 5

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

Egon