Search code examples
rspatiallatitude-longitudeutm

Convert lat/long to UTM (zone 27)

I am struggling to convert my lat and long columns to UTM (zone 27). I found the below function in another post but I believe I am running into an issue due to not all my sites not having spatial data. When I run the below function it does not change my lat long. When I run specific lines I get an error on the first line of the function for the coordinates().

#data

brazil_legacy_data <- structure(list(site = c("1 water", "2 water", "3 water", "4 water", 
"5 water", "6 water", "7 water", "8 water", "9 water", "10 water", 
"11 water", "lago angelim water", "lago anta water", "lago minutal water", 
"lago sede water", "poça temporária water", "trilha water", 
"1 land", "2 land", "3 land", "4 land", "5 land", "6 land", "7 land", 
"8 land", "9 land", "10 land", "11 land", "lago angelim land", 
"lago anta land", "lago minutal land", "lago sede land", "poça temporária land", 
"trilha land"), lat = c(-23.340813, -23.342987, -23.343694, -23.340738, 
-23.332042, -23.324259, -23.327646, -23.313991, NA, -23.386967, 
NA, NA, NA, NA, NA, NA, NA, -23.340744, -23.343771, -23.343741, 
-23.341032, -23.331908, -23.324572, -23.326528, -23.31427, NA, 
-23.387007, NA, NA, NA, NA, NA, NA, NA), long = c(-45.149414, 
-45.144503, -45.129222, -45.117734, -45.099698, -45.089126, -45.133222, 
-45.120015, NA, -45.187918, NA, NA, NA, NA, NA, NA, NA, -45.150486, 
-45.144043, -45.128388, -45.117835, -45.097041, -45.094422, -45.132481, 
-45.124575, NA, -45.186547, NA, NA, NA, NA, NA, NA, NA), region = c("santa virginia", 
"santa virginia", "santa virginia", "santa virginia", "santa virginia", 
"santa virginia", "santa virginia", "santa virginia", "santa virginia", 
"santa virginia", "santa virginia", "santa virginia", "santa virginia", 
"santa virginia", "santa virginia", "santa virginia", "santa virginia", 
"santa virginia", "santa virginia", "santa virginia", "santa virginia", 
"santa virginia", "santa virginia", "santa virginia", "santa virginia", 
"santa virginia", "santa virginia", "santa virginia", "santa virginia", 
"santa virginia", "santa virginia", "santa virginia", "santa virginia", 
"santa virginia"), location = c("brazil", "brazil", "brazil", 
"brazil", "brazil", "brazil", "brazil", "brazil", "brazil", "brazil", 
"brazil", "brazil", "brazil", "brazil", "brazil", "brazil", "brazil", 
"brazil", "brazil", "brazil", "brazil", "brazil", "brazil", "brazil", 
"brazil", "brazil", "brazil", "brazil", "brazil", "brazil", "brazil", 
"brazil", "brazil", "brazil")), row.names = c(NA, -34L), class = c("tbl_df", 
"tbl", "data.frame"))


# convert to utm
long2utm <- function(long)  {
  (floor((long + 180)/6) %% 60) +1
}

LongLatToUTM <- function(brazil_legacy_sites){
    ## Args: df, data frame must have x and y columns. Should be from same UTM zone.
    ## Create a spatial dataframe
    coordinates(brazil_legacy_sites) <- ~lat+long
    proj4string(brazil_legacy_sites) <- CRS("+proj=longlat +datum=WGS84")  

    ## Get zones for all the points in the data frame. 
    ## Stop if more than one zone is present. 
    ## You can write your own code for handling cases where your 
    ## data comes from different UTM zones.

    zone <- long2UTM(brazil_legacy_sites$lat)
    if (length(unique(zone)) > 1) stop("values from different UTM zones")
    zone <- unique(zone)

    ## Change CRS of the spatial data frame and convert to data frame
    res <- spTransform(brazil_legacy_sites, CRS(paste0("+proj=utm +zone=", zone, "+datum=WGS84")))
    return(as.data.frame(res))
}
Solution

  • With these data you only need zone 23

    b <- as.data.frame(brazil_legacy_data)
long2utm(b$long) |> unique()
#[1] 23 NA

    So I will skip the loop

    library(terra)
v <- vect(b, c("long", "lat"), crs="+proj=longlat")
u <- terra::project(v, "+proj=utm +zone=23")
utm <- crds(u)

head(utm)
#            x        y
#[1,] 484726.5 -2581256
#[2,] 485228.8 -2581497
#[3,] 486790.9 -2581573
#[4,] 487964.9 -2581245
#[5,] 489808.0 -2580281
#[6,] 490888.2 -2579419

    All together

    i <- which(!is.na(b$long))
b[i, c("x", "y")] <- utm

head(b)
#     site       lat      long         region location        x        y
#1 1 water -23.34081 -45.14941 santa virginia   brazil 484726.5 -2581256
#2 2 water -23.34299 -45.14450 santa virginia   brazil 485228.8 -2581497
#3 3 water -23.34369 -45.12922 santa virginia   brazil 486790.9 -2581573
#4 4 water -23.34074 -45.11773 santa virginia   brazil 487964.9 -2581245
#5 5 water -23.33204 -45.09970 santa virginia   brazil 489808.0 -2580281
#6 6 water -23.32426 -45.08913 santa virginia   brazil 490888.2 -2579419

    This was not working for you because you had "rgdal" loaded after "terra" and then you get

    rgdal::project(v, "+proj=utm +zone=23")
#Error in rgdal::project(v, "+proj=utm +zone=23") : xy not numeric