UTM Zone conversion (#139)

* add utm_convert function

* document

* fix column extraction

* Add ability to keep xy cols and use a single zone

* Add tests, fix single zone

* Incorporate WGS84 and NAD83, update tests and docs

* Use rbind instead of dplyr::bind_rows, update NEWS

* add script for getting UTM zone/epsg lookup table

* deal with other invalid zones

* fix tests

* fix typo

Co-authored-by: Sam Albers <sam.albers@gov.bc.ca>

* change args to easting and northing

* update pkgdown with utm_convert entry

---------

Co-authored-by: Sam Albers <sam.albers@gov.bc.ca>

NAMESPACE

@@ -53,6 +53,7 @@ export(show_cached_files)
 export(summarize_raster_list)
 export(transform_bc_albers)
 export(tsa)
+export(utm_convert)
 export(vrt_files)
 export(vrt_info)
 export(water_districts)

NEWS.md

@@ -1,5 +1,8 @@
 # bcmaps (development version)
 
+* Added function `utm_convert()` to convert tabular data with X and Y coordinates
+  in (possibly multiple) UTM zones to a single CRS.
+
 # bcmaps 2.1.0
 
 * Added function `cded_terra()`

R/utm-convert.R

@@ -0,0 +1,103 @@
+#' Convert a data.frame of UTM coordinates to an sf object with a single CRS
+#'
+#' This can operate on a data frame containing coordinates from multiple UTM zones
+#' with a column denoting the zone, or a single zone for the full dataset.
+#'
+#' It supports data collected in either the NAD83 or WGS84 ellipsoid in
+#' the Northern hemisphere
+#'
+#' @param x data.frame containing UTM coordinates, with a zone column
+#' @param easting the name of the 'easting' column
+#' @param northing the name of the 'northing' column
+#' @param zone the name of the 'zone' column, or a single value if
+#'    the data are all in one UTM zone
+#' @param crs target CRS. Default BC Albers (EPSG:3005)
+#' @param datum The datum of the source data. `"NAD83"` (Default) or `"WGS84"`
+#' @param xycols should the X and Y columns be appended to the output? `TRUE` or `FALSE`
+#'
+#' @return sf object in the chosen CRS
+#' @export
+#' @examples
+#' # Data with multiple zones, and a column denoting the zone
+#' df <- data.frame(
+#'   animalid = c("a", "b", "c"),
+#'   zone = c(10, 11, 11),
+#'   easting = c(500000, 800000, 700000),
+#'   northing = c(5000000, 3000000, 1000000)
+#' )
+#' utm_convert(df, easting = "easting", northing = "northing", zone = "zone")
+#'
+#' # Data all in one zone, specify a single zone:
+#' df <- data.frame(
+#'   animalid = c("a", "b"),
+#'   easting = c(500000, 800000),
+#'   northing = c(5000000, 3000000)
+#' )
+#' utm_convert(df, easting = "easting", northing = "northing", zone = 11)
+utm_convert <- function(x, easting, northing, zone, crs = "EPSG:3005",
+                        datum = c("NAD83", "WGS84"), xycols = TRUE) {
+
+  one_zone <- !zone %in% names(x) && is.numeric(format_zone(zone))
+
+  datum <- match.arg(datum)
+
+  if (!one_zone && !zone %in% names(x)) {
+    stop(zone, " is not a column in x", call. = FALSE)
+  }
+  if (!easting %in% names(x)) {
+    stop(easting, " is not a column in x", call. = FALSE)
+  }
+  if (!northing %in% names(x)) {
+    stop(northing, " is not a column in x", call. = FALSE)
+  }
+
+  if (one_zone) {
+    res <- convert_from_zone(x, zone, easting, northing, crs, datum, xycols)
+    return(cbind(res, x[, setdiff(names(x), names(res))]))
+  }
+
+  x_split <- split(x, x[zone])
+
+  x_split <- lapply(x_split, function(z) {
+    zone <- z[[zone]][1]
+    convert_from_zone(z, zone, easting, northing, crs, datum, xycols)
+  })
+
+  res <- do.call("rbind", x_split)
+  cbind(res, x[, setdiff(names(x), names(res))])
+}
+
+convert_from_zone <- function(x, zone, easting, northing, crs, datum, xycols) {
+  epsg <- lookup_epsg_code(zone, datum)
+  x <- sf::st_as_sf(x, coords = c(easting, northing), crs = epsg)
+  res <- sf::st_transform(x, crs = crs)
+  if (xycols) {
+    res <- cbind(res, sf::st_coordinates(res))
+  }
+  res
+}
+
+lookup_epsg_code <- function(zone, datum) {
+  zone <- format_zone(zone)
+  if (!zone %in% utm_zone_lookup[["zone_numeric"]]) {
+    stop("Invalid zone: ", zone, call. = FALSE)
+  }
+  utm_zone_lookup[["epsg_code"]][
+    utm_zone_lookup[["zone_numeric"]] == zone &
+      utm_zone_lookup[["datum"]] == datum
+    ]
+}
+
+format_zone <- function(x) {
+  if (is.numeric(x)) return(x)
+
+  if (grepl("[Ss]$", x)) {
+    stop("It looks like your UTM zones are in the Southern hemisphere", call. = FALSE)
+  }
+
+  ret <- suppressWarnings(as.numeric(gsub("[NnUu$]", "", x)))
+  if (is.na(ret)) {
+    stop("Invalid zone(s): ", x, call. = FALSE)
+  }
+  ret
+}

_pkgdown.yml

@@ -33,6 +33,7 @@ reference:
     - bc_area
     - bc_bbox
     - transform_bc_albers
+    - utm_convert
     - combine_nr_rd
     - bc_neighbours
     - bec_colours

data-raw/create_internal_data.R

@@ -13,11 +13,19 @@
 ## Code to prepare the layers_df internal data
 layers_df <- dplyr::arrange(readr::read_csv("data-raw/layers_df.csv"), layer_name)
 
+utm_zone_lookup <- readr::read_csv("data-raw/utm-zone-lookup.csv")
 
 source("data-raw/mapsheets_250K/process_mapsheet_250K.R")
 
 # layers_df <- layers_df[!is.na(layers_df$record),]
 # layers_df <- layers_df[!layers_df$local,]
 
-usethis::use_data(layers_df, mapsheets_250K_data, mapsheets_50K_data,
-                  overwrite = TRUE, internal = TRUE, compress = "xz")
+usethis::use_data(
+  layers_df,
+  mapsheets_250K_data,
+  mapsheets_50K_data,
+  utm_zone_lookup,
+  overwrite = TRUE,
+  internal = TRUE,
+  compress = "xz"
+)

data-raw/utm-zone-lookup.R

@@ -0,0 +1,20 @@
+# Don't need to run this again, but this is how we got the utm zone-lookup table
+
+library(rvest)
+library(janitor)
+library(dplyr)
+
+html <- read_html("https://docs.up42.com/data/reference/utm")
+
+nad83_lookup <- html_table(html)[[3]] |>
+  clean_names() |>
+  mutate(zone_numeric = as.numeric(gsub("N", "", utm_zone)),
+         datum = "NAD83")
+
+wgs84_lookup <- html_table(html)[[1]] |>
+  clean_names() |>
+  mutate(zone_numeric = as.numeric(gsub("N", "", utm_zone)),
+         datum = "WGS84")
+
+bind_rows(nad83_lookup, wgs84_lookup) |>
+  readr::write_csv("data-raw/zone-lookup.csv")

data-raw/utm-zone-lookup.csv

@@ -0,0 +1,84 @@
+utm_zone,coordinate_reference_system,epsg_code,zone_numeric,datum
+1N,NAD83 / UTM zone 1N,EPSG:26901,1,NAD83
+2N,NAD83 / UTM zone 2N,EPSG:26902,2,NAD83
+3N,NAD83 / UTM zone 3N,EPSG:26903,3,NAD83
+4N,NAD83 / UTM zone 4N,EPSG:26904,4,NAD83
+5N,NAD83 / UTM zone 5N,EPSG:26905,5,NAD83
+6N,NAD83 / UTM zone 6N,EPSG:26906,6,NAD83
+7N,NAD83 / UTM zone 7N,EPSG:26907,7,NAD83
+8N,NAD83 / UTM zone 8N,EPSG:26908,8,NAD83
+9N,NAD83 / UTM zone 9N,EPSG:26909,9,NAD83
+10N,NAD83 / UTM zone 10N,EPSG:26910,10,NAD83
+11N,NAD83 / UTM zone 11N,EPSG:26911,11,NAD83
+12N,NAD83 / UTM zone 12N,EPSG:26912,12,NAD83
+13N,NAD83 / UTM zone 13N,EPSG:26913,13,NAD83
+14N,NAD83 / UTM zone 14N,EPSG:26914,14,NAD83
+15N,NAD83 / UTM zone 15N,EPSG:26915,15,NAD83
+16N,NAD83 / UTM zone 16N,EPSG:26916,16,NAD83
+17N,NAD83 / UTM zone 17N,EPSG:26917,17,NAD83
+18N,NAD83 / UTM zone 18N,EPSG:26918,18,NAD83
+19N,NAD83 / UTM zone 19N,EPSG:26919,19,NAD83
+20N,NAD83 / UTM zone 20N,EPSG:26920,20,NAD83
+21N,NAD83 / UTM zone 21N,EPSG:26921,21,NAD83
+22N,NAD83 / UTM zone 22N,EPSG:26922,22,NAD83
+23N,NAD83 / UTM zone 23N,EPSG:26923,23,NAD83
+1N,WGS84 / UTM zone 1N,EPSG:32601,1,WGS84
+2N,WGS84 / UTM zone 2N,EPSG:32602,2,WGS84
+3N,WGS84 / UTM zone 3N,EPSG:32603,3,WGS84
+4N,WGS84 / UTM zone 4N,EPSG:32604,4,WGS84
+5N,WGS84 / UTM zone 5N,EPSG:32605,5,WGS84
+6N,WGS84 / UTM zone 6N,EPSG:32606,6,WGS84
+7N,WGS84 / UTM zone 7N,EPSG:32607,7,WGS84
+8N,WGS84 / UTM zone 8N,EPSG:32608,8,WGS84
+9N,WGS84 / UTM zone 9N,EPSG:32609,9,WGS84
+10N,WGS84 / UTM zone 10N,EPSG:32610,10,WGS84
+11N,WGS84 / UTM zone 11N,EPSG:32611,11,WGS84
+12N,WGS84 / UTM zone 12N,EPSG:32612,12,WGS84
+13N,WGS84 / UTM zone 13N,EPSG:32613,13,WGS84
+14N,WGS84 / UTM zone 14N,EPSG:32614,14,WGS84
+15N,WGS84 / UTM zone 15N,EPSG:32615,15,WGS84
+16N,WGS84 / UTM zone 16N,EPSG:32616,16,WGS84
+17N,WGS84 / UTM zone 17N,EPSG:32617,17,WGS84
+18N,WGS84 / UTM zone 18N,EPSG:32618,18,WGS84
+19N,WGS84 / UTM zone 19N,EPSG:32619,19,WGS84
+20N,WGS84 / UTM zone 20N,EPSG:32620,20,WGS84
+21N,WGS84 / UTM zone 21N,EPSG:32621,21,WGS84
+22N,WGS84 / UTM zone 22N,EPSG:32622,22,WGS84
+23N,WGS84 / UTM zone 23N,EPSG:32623,23,WGS84
+24N,WGS84 / UTM zone 24N,EPSG:32624,24,WGS84
+25N,WGS84 / UTM zone 25N,EPSG:32625,25,WGS84
+26N,WGS84 / UTM zone 26N,EPSG:32626,26,WGS84
+27N,WGS84 / UTM zone 27N,EPSG:32627,27,WGS84
+28N,WGS84 / UTM zone 28N,EPSG:32628,28,WGS84
+29N,WGS84 / UTM zone 29N,EPSG:32629,29,WGS84
+30N,WGS84 / UTM zone 30N,EPSG:32630,30,WGS84
+31N,WGS84 / UTM zone 31N,EPSG:32631,31,WGS84
+32N,WGS84 / UTM zone 32N,EPSG:32632,32,WGS84
+33N,WGS84 / UTM zone 33N,EPSG:32633,33,WGS84
+34N,WGS84 / UTM zone 34N,EPSG:32634,34,WGS84
+35N,WGS84 / UTM zone 35N,EPSG:32635,35,WGS84
+36N,WGS84 / UTM zone 36N,EPSG:32636,36,WGS84
+37N,WGS84 / UTM zone 37N,EPSG:32637,37,WGS84
+38N,WGS84 / UTM zone 38N,EPSG:32638,38,WGS84
+39N,WGS84 / UTM zone 39N,EPSG:32639,39,WGS84
+40N,WGS84 / UTM zone 40N,EPSG:32640,40,WGS84
+41N,WGS84 / UTM zone 41N,EPSG:32641,41,WGS84
+42N,WGS84 / UTM zone 42N,EPSG:32642,42,WGS84
+43N,WGS84 / UTM zone 43N,EPSG:32643,43,WGS84
+44N,WGS84 / UTM zone 44N,EPSG:32644,44,WGS84
+45N,WGS84 / UTM zone 45N,EPSG:32645,45,WGS84
+46N,WGS84 / UTM zone 46N,EPSG:32646,46,WGS84
+47N,WGS84 / UTM zone 47N,EPSG:32647,47,WGS84
+48N,WGS84 / UTM zone 48N,EPSG:32648,48,WGS84
+49N,WGS84 / UTM zone 49N,EPSG:32649,49,WGS84
+50N,WGS84 / UTM zone 50N,EPSG:32650,50,WGS84
+51N,WGS84 / UTM zone 51N,EPSG:32651,51,WGS84
+52N,WGS84 / UTM zone 52N,EPSG:32652,52,WGS84
+53N,WGS84 / UTM zone 53N,EPSG:32653,53,WGS84
+54N,WGS84 / UTM zone 54N,EPSG:32654,54,WGS84
+55N,WGS84 / UTM zone 55N,EPSG:32655,55,WGS84
+56N,WGS84 / UTM zone 56N,EPSG:32656,56,WGS84
+57N,WGS84 / UTM zone 57N,EPSG:32657,57,WGS84
+58N,WGS84 / UTM zone 58N,EPSG:32658,58,WGS84
+59N,WGS84 / UTM zone 59N,EPSG:32659,59,WGS84
+60N,WGS84 / UTM zone 60N,EPSG:32660,60,WGS84