TITLE: Converting stem locations from lat-long to XY coordinates
DATE: 2024-06-12
AUTHOR: John L. Godlee
====================================================================


In 2018 I wrote some code to convert stem locations within square
plots from lat-long coordinates into XY grid coordinates. Recently
I was asked to do the opposite, to convert XY grid coordinates into
lat-long coordinates, so a colleague could combine multiple
adjacent plots into one larger plot with a shared coordinate
system. I've pasted the code for this procedure below:

   # Convert XY grid coordinates to Lat-long coordinates
   # John L. Godlee ([email protected])
   # Last updated: 2024-06-12

   # Packages
   library(dplyr)
   library(sf)
   library(NISTunits)
   library(geosphere)
   library(ggplot2)

   # Define functions

   #' Get valid UTM zone from latitude and longitude in WGS84
decimal degrees
   #'
   #' @param x vector of longitude coordinates in decimal degrees
   #' @param y vector of latitude coordinate in decimal degrees
   #'
   #' @return Vector of UTM zones for each latitude-longitude pair
   #'
   #' @export
   #'
   latLong2UTM <- function(x, y) {
     unlist(lapply(1:length(x), function(z) {
       paste((floor((as.numeric(x[z]) + 180) / 6) %% 60) + 1,
         ifelse(as.numeric(y[z]) < 0, "S", "N"),
         sep = "")
     }))
   }

   #' Generate a valid UTM WGS84 proj4string given a UTM zone
   #'
   #' @param x character vector defining UTM zones
   #'
   #' @return UTM proj4string character vector
   #'
   #' @export
   #'
   UTMProj4 <- function(x){
     unlist(lapply(1:length(x), function(y) {
       paste0(
         "+proj=utm +zone=",
         gsub("[A-z]", "", as.character(x[y])),
         ifelse(gsub("[0-9]", "", as.character(x[y])) == "S", "
+south", ""),
         " +ellps=WGS84")
     }))
   }

   #' Perform a rotation on an sf object
   #'
   #' @param a angle to rotate, in radians
   #'
   #' @keywords internal
   #' @noRd
   #'
   rot <- function(a) {
     matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
   }

   # Define UTM zone for test polygon location
   utm_crs <- "+proj=utm +zone=33 +south +ellps=WGS84"

   # Define origin corner (SW) of test polygon
   orig <- c(479000, 8319000)

   # Create dataframe of polygon corners (100x100 m)
   poly_df <- data.frame(
     longitude = c(orig[1], orig[1], orig[1] + 100, orig[1] + 100,
orig[1]),
     latitude = c(orig[2], orig[2] + 100, orig[2] + 100, orig[2],
orig[2]))

   # Convert dataframe of corners to a polygon
   poly_ex <- st_convex_hull(summarise(st_as_sf(poly_df,
         coords = c("longitude", "latitude"), crs = utm_crs)))

   # Get the centroid of the polygon and rotate a bit
   poly_ex_cent <- st_centroid(st_geometry(poly_ex))

   rot_ex <- sample(seq(0, 360, 1), 1)

   poly_ex_rot <- st_sf(geometry = (st_geometry(poly_ex) -
poly_ex_cent) *
     rot(NISTdegTOradian(rot_ex)) * 1 + poly_ex_cent)

   st_crs(poly_ex_rot) <- utm_crs

   # Transform rotated polygon to UTM
   poly_ex_wgs <- st_transform(poly_ex_rot, 4326)
   poly_ex_wgs$plot_id <- "A"

   # Sample points within the rotated polygon
   xy_range <- seq(0, 100, 0.5)
   s_ex <- data.frame(
     plot_id = "A",
     x_grid = sample(xy_range, 50),
     y_grid = sample(xy_range, 50))

   # Rename polygon and stem data to test code below.
   poly_data <- poly_ex_wgs
   x <- s_ex

   # Add ID column
   x$id <- seq_len(nrow(x))

   # Get UTM zone of polygon centroid
   poly_cent <- st_coordinates(st_centroid(poly_data))
   utm_string <- UTMProj4(latLong2UTM(poly_cent[1], poly_cent[2]))

   # Convert polygon to UTM
   poly_utm <- st_transform(poly_data, utm_string)

   # Convert UTM polygon to points
   points_utm <- st_cast(poly_utm, "POINT", warn = FALSE)[1:4,]

   # Extract coordinates as dataframe
   coords_utm <- as.data.frame(st_coordinates(points_utm)[1:4,1:2])

   # Define points to match corners to
   match_point <- st_sfc(st_point(
       x = c(mean(coords_utm$X) - 1000, mean(coords_utm$Y) -
1000)))
   other_point <- st_sfc(st_point(
       x = c(mean(coords_utm$X) + 1000, mean(coords_utm$Y) -
1000)))

   # Set CRS
   st_crs(other_point) <- utm_string
   st_crs(match_point) <- utm_string

   # Get sw and se corner
   sw_corner <- points_utm[st_nearest_feature(match_point,
points_utm),]
   se_corner <- points_utm[st_nearest_feature(other_point,
points_utm),]

   # Convert to WGS for geosphere compatibility
   sw_wgs <- st_coordinates(st_transform(sw_corner, 4326))
   se_wgs <- st_coordinates(st_transform(se_corner, 4326))

   # Find rotation along SW,SE edge
   xy_bearing <- bearing(sw_wgs, se_wgs)

   # Get centroid of polygon in UTM
   cent <- suppressWarnings(st_centroid(st_geometry(poly_utm)))

   # Get location of sw corner
   sw_cent <- st_geometry(sw_corner)

   # Convert angle to radians for rotation
   angle <- NISTdegTOradian(90 - xy_bearing)

   # Rotate polygon so SW-SE is flat
   poly_rot <- (st_geometry(poly_utm) - cent) * rot(angle) * 1 +
cent
   st_crs(poly_rot) <- st_crs(poly_utm)

   # Add SW UTM X and Y to XY grid coordinates
   points_rot <- st_cast(poly_rot, "POINT", warn = FALSE)[1:4]
   sw_corner_rot <- points_rot[st_nearest_feature(match_point,
points_rot)]
   se_corner_rot <- points_rot[st_nearest_feature(other_point,
points_rot)]
   x$x_grid_utm <- x$x_grid + st_coordinates(sw_corner_rot)[1]
   x$y_grid_utm <- x$y_grid + st_coordinates(sw_corner_rot)[2]

   x_sf <- st_as_sf(
     x[,c("id", "x_grid_utm", "y_grid_utm")],
     coords = c("x_grid_utm", "y_grid_utm"),
     na.fail = FALSE, crs = st_crs(poly_rot))

   # Rotate coordinates back
   x_sf_rot <- (st_geometry(x_sf) - cent) * rot(-angle) * 1 + cent
   st_crs(x_sf_rot) <- st_crs(x_sf)

   # Transform stems to WGS84
   x_wgs <- st_transform(x_sf_rot, 4326)
   x_out <- as.data.frame(cbind(x$id, st_coordinates(x_wgs)))
   names(x_out) <- c("id", "longitude", "latitude")

   # Clean dataframe
   out <- x_out[order(x_out$id),c("longitude", "latitude")]
   out$id <- NULL

   # Create plot to illustrate code above
   ggplot() +
    geom_sf(data = poly_utm) +
    geom_sf(data = sw_corner, colour = "red") +
    geom_sf(data = se_corner, colour = "blue") +
    geom_sf(data = cent, colour = "cyan") +
    geom_sf(data = poly_rot, colour = "green", fill = NA) +
    geom_sf(data = x_sf) +
    geom_sf(data = x_sf_rot, colour = "pink")

   # This dataframe contains lat-long coordinates for each row of x
   out

 ![Demonstration of the code to converrt lat-long coordinates to
XY
coordinates.](https://johngodlee.xyz/img_full/xy_latlong/xy_illus.pn
g)