Search code examples
rcoordinatesgeospatialr-sf

Converting UTM to decimal degrees with st_transform (a southern hemisphere error)


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?


Solution

  • 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