library(tidyverse)
library(sf)
x <- 500000
y <- 4649776
zone_num <- "33"
ellipsoid <- "WGS84"
res <- data.frame(x, y) %>%
sf::st_as_sf(coords = c("x", "y"),
crs = paste0("+proj=utm +zone=", zone_num)) %>%
sf::st_transform(crs = paste0("+proj=longlat +datum=", ellipsoid)) %>%
sf::st_coordinates() %>%
data.frame() %>%
dplyr::rename("lon" = "X", "lat" = "Y");res
# Test case 2: Southern hemisphere
x2 <- 500000
y2 <- (4649776 +10000000)-1 # +10000000)-1 tp adjust for southern hemisphere latitude
res <- data.frame(x2, y2) %>%
sf::st_as_sf(coords = c("x2", "y2"),
crs = paste0("+proj=utm +zone=", zone_num)) %>%
sf::st_transform(crs = paste0("+proj=longlat +datum=", ellipsoid)) %>%
sf::st_coordinates() %>%
data.frame() %>%
dplyr::rename("lon" = "X", "lat" = "Y");res
For the first the output is correct lon: 15 lat: 42
For the southern hemisphere it incorrect as it should be lon: 15 lat: -48.30
but the output is this: lon: -165 lat: 48.2686
I don't understand why this is happening, the longitude should be unchanged, and the latitude is nearly 0.04 degrees off. Any idea why sf is having trouble transforming?
A couple of issues here. The main one is that your proj4string
for the southern hemisphere should have the +south
parameter.
The second issue is that the Northing range in UTM is between 0 and around 10 million. You are adding 10 million to your original Northing (y
) value so you end up with a number outside the valid range.
If you fix these issues you will get the same lon
value as in the northern hemisphere and the appropriate lat
value:
x2 <- x
y2 <- (y + 5e6) - 1 # only adding 5m
data.frame(x2, y2) %>%
sf::st_as_sf(
coords = c("x2", "y2"),
crs = paste0("+proj=utm +zone=", zone_num, " +south") # Note " + south"
) %>%
sf::st_transform(crs = paste0("+proj=longlat +datum=", ellipsoid)) %>%
sf::st_coordinates() %>%
data.frame() %>%
dplyr::rename("lon" = "X", "lat" = "Y")
# lon lat
# 1 15 -3.168563