How to convert projected coordinates back to lat and long
Asked Answered
C

2

5

I am working with lat long coordinates in R. First I need to transform my coordinates to zone 17 and then add some buffer to them. For that I used next code (the data used is in orig_coords ):

library(sp)
library(rgdal)
library(rgeos)

orig_coords <- data.frame(lat=-0.1848787,lon=-78.48179)

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 +south=T +zone=",zone," ellps=WGS84",sep='')))
  
  #return(as.data.frame(res))
  return(res)
  
}
#Res
res <- LongLatToUTM(orig_coords$lon,orig_coords$lat,17)
#Add buffer
pc2km <- gBuffer( res, width=2000, byid=TRUE )
#Bounding box
bb <- pc2km@bbox

With that I got the bounding box which contains the projected coordinates with the buffer added to them:

bb
        min       max
x  778303.1  782303.1
y 9977545.5 9981545.5

My question is how can I transform the min and max coordinates from bb back to latitude and longitude? Many thanks.

Chiliarch answered 16/6 at 22:35 Comment(2)
I don't have a version of R that I can easily install the required packages on to test this, but you could try projecting your buffer object back to WGS84 before creating your bb object. Something like pc2km_WGS84 <- spTransform(pc2km, CRS("+proj=longlat +datum=WGS84")) and then bb <- pc2km_WGS84@bbox. Not sure if that syntax is correct as I have not really used those packages much, but the workflow seems logical.Nostoc
It is likely worth contemplating updating your geo tools to sf, with gdal and geos installed at system level as both rgdal and rgeos have both been retired/archived.Markowitz
P
4

Here is the answer using sf:

orig_coords <- data.frame(lat=-0.1848787,lon=-78.48179)

library(sf)

# Create the sf object for the points
orig_sf <- st_sfc(st_point(c(orig_coords$lon, orig_coords$lat)), crs = "epsg:4326")
> orig_sf
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -78.48179 ymin: -0.1848787 xmax: -78.48179 ymax: -0.1848787
Geodetic CRS:  WGS 84
POINT (-78.48179 -0.1848787)

# Transform to 32617 based on "+proj=utm +south=T +zone=17 +ellps=WGS84" and https://epsg.io/32717
utm_sf <- st_transform(orig_sf, crs = "epsg:32717")
> utm_sf
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 780303.1 ymin: 9979545 xmax: 780303.1 ymax: 9979545
Projected CRS: WGS 84 / UTM zone 17S
POINT (780303.1 9979545)

# Create buffer
utm_buff_sf <- st_buffer(utm_sf, 2000)

# Print buffer coordinates. You can also see them in the object itself
st_bbox(utm_buff_sf)
     xmin      ymin      xmax      ymax 
 778303.1 9977545.5  782303.1 9981545.5 

# Transform back to lat & long
latlon_sf <- st_transform(utm_buff_sf, crs = "epsg:4326")
st_bbox(latlon_sf)
       xmin        ymin        xmax        ymax 
-78.4997462  -0.2029557 -78.4638340  -0.1668017 

And this is the answer with terra:

library(terra)

# Could also use vect(<sf object>), if the sf object already exists.
orig_te <- vect(orig_coords, crs = "epsg:4326")
> orig_te
 class       : SpatVector 
 geometry    : points 
 dimensions  : 1, 0  (geometries, attributes)
 extent      : -78.48179, -78.48179, -0.1848787, -0.1848787  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 

# "Transform" 
utm_te <- project(orig_te, y = "epsg:32717")

# Buffer
buff_te <- buffer(utm_te, 2000)

# Get bounding box
bbox_te <- ext(buff_te)
bbox_te
SpatExtent : 778303.140207565, 782303.140207565, 9977545.46221348, 9981545.46221348 (xmin, xmax, ymin, ymax)

# Transform back to lat & long
latlon_te <- project(buff_te, y = "epsg:4326")
ext(latlon_te)
SpatExtent : -78.4997462383223, -78.4638340108377, -0.202955733632546, -0.16680166596275 (xmin, xmax, ymin, ymax)

The coordinates of the projected buffer's bounding box using either sf::st_transform() or terra::project() match those of the coordinates supplied with your custom function.

Prank answered 10/8 at 15:12 Comment(1)
Amazing, let's wait a couple of days and I will accept this answer. Many thanks! The results are the same as the original code!Chiliarch
C
5

Your data

orig_coords <- data.frame(lat=-0.1848787,lon=-78.48179)

The solution with terra (same as John Polo's but with more compact notation)

library(terra)
zone <- 17
vect(orig_coords, crs="+proj=longlat") |>
    project(paste0("+proj=utm +south +zone=", zone)) |> 
    buffer(2000) |> project("+proj=longlat") |> 
    ext() |> as.vector()

#       xmin        xmax        ymin        ymax 
#-78.4997462 -78.4638340  -0.2029557  -0.1668017 

But I would suggest a more direct approach. You can skip the coordinate transformations and compute the buffer directly on the lon/lat coordinates:

vect(orig_coords, crs="+proj=longlat") |> 
    buffer(2000) |> ext() |> as.vector()

#       xmin        xmax        ymin        ymax 
#-78.4997564 -78.4638236  -0.2029661  -0.1667913 

The difference between these two methods should be small, as long as you use the correct UTM zone and the buffer is small. But the direct approach is more accurate as it does not not suffer from the distortion introduced by the projection; and it is simpler as you do not need to specify a UTM zone.

Czarevna answered 13/8 at 9:12 Comment(0)
P
4

Here is the answer using sf:

orig_coords <- data.frame(lat=-0.1848787,lon=-78.48179)

library(sf)

# Create the sf object for the points
orig_sf <- st_sfc(st_point(c(orig_coords$lon, orig_coords$lat)), crs = "epsg:4326")
> orig_sf
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -78.48179 ymin: -0.1848787 xmax: -78.48179 ymax: -0.1848787
Geodetic CRS:  WGS 84
POINT (-78.48179 -0.1848787)

# Transform to 32617 based on "+proj=utm +south=T +zone=17 +ellps=WGS84" and https://epsg.io/32717
utm_sf <- st_transform(orig_sf, crs = "epsg:32717")
> utm_sf
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 780303.1 ymin: 9979545 xmax: 780303.1 ymax: 9979545
Projected CRS: WGS 84 / UTM zone 17S
POINT (780303.1 9979545)

# Create buffer
utm_buff_sf <- st_buffer(utm_sf, 2000)

# Print buffer coordinates. You can also see them in the object itself
st_bbox(utm_buff_sf)
     xmin      ymin      xmax      ymax 
 778303.1 9977545.5  782303.1 9981545.5 

# Transform back to lat & long
latlon_sf <- st_transform(utm_buff_sf, crs = "epsg:4326")
st_bbox(latlon_sf)
       xmin        ymin        xmax        ymax 
-78.4997462  -0.2029557 -78.4638340  -0.1668017 

And this is the answer with terra:

library(terra)

# Could also use vect(<sf object>), if the sf object already exists.
orig_te <- vect(orig_coords, crs = "epsg:4326")
> orig_te
 class       : SpatVector 
 geometry    : points 
 dimensions  : 1, 0  (geometries, attributes)
 extent      : -78.48179, -78.48179, -0.1848787, -0.1848787  (xmin, xmax, ymin, ymax)
 coord. ref. : lon/lat WGS 84 (EPSG:4326) 

# "Transform" 
utm_te <- project(orig_te, y = "epsg:32717")

# Buffer
buff_te <- buffer(utm_te, 2000)

# Get bounding box
bbox_te <- ext(buff_te)
bbox_te
SpatExtent : 778303.140207565, 782303.140207565, 9977545.46221348, 9981545.46221348 (xmin, xmax, ymin, ymax)

# Transform back to lat & long
latlon_te <- project(buff_te, y = "epsg:4326")
ext(latlon_te)
SpatExtent : -78.4997462383223, -78.4638340108377, -0.202955733632546, -0.16680166596275 (xmin, xmax, ymin, ymax)

The coordinates of the projected buffer's bounding box using either sf::st_transform() or terra::project() match those of the coordinates supplied with your custom function.

Prank answered 10/8 at 15:12 Comment(1)
Amazing, let's wait a couple of days and I will accept this answer. Many thanks! The results are the same as the original code!Chiliarch

© 2022 - 2024 — McMap. All rights reserved.