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 (johngodlee@gmail.com)
# 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.
