An R function to convert an sf/sfc object to an appropriate utm or laea projection.

to-generic-projected.R

#' Project an sf/sfc object to a generic projected coordinate system
#' @param x an sf or sfc object
#' @param proj a character vector. The projection to use. One of "laea", "aeqd",
#' "utm", "pconic", or "eqdc".
#' @param ellps a character vector. The ellipsoid to use. Select from
#' `sf_proj_info(type = "ellps")`.
#' @param opts a character vector. Additional proj options to pass to the
#' proj string. see details for more information.
#' @return an sf or sfc object
#' @details For further info about the available "generic" projects see:
#' for utm: \url{https://proj.org/en/9.4/operations/projections/utm.html}
#' for laea: \url{https://proj.org/en/9.4/operations/projections/laea.html}
#' for aeqd: \url{https://proj.org/en/9.4/operations/projections/aeqd.html}
#' for pconic: \url{https://proj.org/en/9.4/operations/projections/pconic.html}
#' for eqdc: \url{https://proj.org/en/9.4/operations/projections/eqdc.html}
to_generic_projected <- function(
    x,
    proj = c("laea", "aeqd", "utm", "pconic", "eqdc"),
    ellps = "WGS84",
    no_defs = TRUE,
    opts = "",
    return_as = c("sf", "proj4", "wkt")) {
  # arg assertions
  if (!rlang::is_true(rlang::inherits_any(x, c("sf", "sfc")))) {
    rlang::abort(
      "`x` must be either an sf object or an sfc object",
      class = "to_generic_projected_incorrect_input"
    )
  }

  proj <- rlang::arg_match(proj)
  ellps <- rlang::arg_match(ellps, sf::sf_proj_info(type = "ellps")$name)

  if (!rlang::is_logical(no_defs)) {
    rlang::abort("`no_defs` must be a logical value",
      class = "to_generic_projected_incorrect_input"
    )
  }

  if (!rlang::is_character(opts) && !nchar(opts)) {
    rlang::abort("`opts` must be a character vector",
      class = "to_generic_projected_incorrect_input"
    )
  }

  return_as <- rlang::arg_match(return_as)

  # get centroid in latlong
  cent_coor <- sf::sf_project(
    sf::st_crs(x), "EPSG:4326",
    sf::st_coordinates(sf::st_centroid(sf::st_geometry(x)))
  )

  # configure proj args
  n_or_s <- ifelse(cent_coor[2] == 0, "",
    ifelse(cent_coor[2] > 0, "+north", "+south")
  )

  no_defs <- ifelse(no_defs, "+no_defs", "")

  if (proj %in% c("pconic", "eqdc")) {
    bounds <- sf::st_bbox(sf::st_transform(x, 4326))
    lat_1 <- bounds$ymax
    lat_2 <- bounds$ymin
  }

  # construct proj4 string
  prj <- trimws(switch(proj,
    laea = glue::glue(
      "+proj=laea +lon_0={cent_coor[1]} +lat_0={cent_coor[2]}",
      "+ellps={ellps} {no_defs}",
      opts,
      .sep = " "
    ),
    utm = glue::glue(
      "+proj=utm +zone={round((180 + cent_coor[1]) / 6)} {n_or_s}",
      "+ellps={ellps} {no_defs}",
      opts,
      .sep = " "
    ),
    aeqd = glue::glue(
      "+proj=aeqd +lon_0={cent_coor[1]} +lat_0={cent_coor[2]}",
      "+ellps={ellps} {no_defs}",
      opts,
      .sep = " "
    ),
    pconic = glue::glue(
      "+proj=pconic +lon_0={cent_coor[1]} +lat_0={cent_coor[2]}",
      "+lat_1={lat_1} +lat_2={lat_2}",
      "+ellps={ellps} {no_defs}",
      opts,
      .sep = " "
    ),
    eqdc = glue::glue(
      "+proj=eqdc +lon_0={cent_coor[1]}",
      "+lat_1={lat_1} +lat_2={lat_2}",
      "+ellps={ellps} {no_defs}",
      opts,
      .sep = " "
    )
  ))

  switch(return_as,
    sf = sf::st_transform(x, prj),
    proj4 = prj,
    wkt = sf::st_crs(prj)$wkt
  )
}

yyy-geo-bounds.R

#' @title get administritive outlines for a country
#' @description using the geoBoundaires API,
#' download the geojson outline of a country
#' @param country character vector of country names
#' @param admin_level character vector of admin levels to download
#' @return sf object of the outlines
#' @details check out the documentation for the geoboundaries API at:
#' geoBoundaries.org
#'
geo_bounds <- function(country, admin_level = c("ADM0", "ADM1", "ADM2")) {
  country <- countrycode::countrycode(country,
    origin = "country.name",
    destination = "iso3c"
  )
  url <- paste(
    "https://www.geoboundaries.org/api/current/gbOpen/",
    country, admin_level[1],
    sep = "/"
  )
  get <- httr::GET(url)
  cont <- httr::content(get, as = "parsed")
  area <- sf::read_sf(cont$gjDownloadURL)
  return(area)
}

zzz-main.R

source("to-generic-projected.R")
source("yyy-geo-bounds.R")


uk <- geo_bounds("United Kingdom")
x <- sf::st_geometry(uk)

plot_compare <- function(x, wgs_datum = FALSE) {
  plot(x,
    axes = TRUE, graticule = TRUE, main = "latlong",
    col_graticule = "#76c092"
  )
  p <- lapply(c("laea", "utm", "aeqd", "pconic", "eqdc"), \(proj) {
    y <- to_generic_projected(x, proj)
    if (wgs_datum) {
      datum <- sf::st_crs(4326)
    } else {
      datum <- NULL
    }
    plot(y,
      axes = TRUE, main = proj,
      graticule = sf::st_graticule(y, datum = datum),
      col_graticule = "#76c092"
    )
  })
}

if (requireNamespace("basetheme", quietly = TRUE)) {
  basetheme::basetheme("dark")
}
par(mfrow = c(2, 3))
plot_compare(x)
plot_compare(x, wgs_datum = TRUE)


sapply(c("laea", "utm", "aeqd", "pconic", "eqdc"), \(proj) {
  y <- to_generic_projected(
    uk,
    proj,
    return_as = "proj4"
  )
})

#                                                                                                                                  laea
#                                                   "+proj=laea +lon_0=-2.78488515753149 +lat_0=54.0378529685846 +ellps=WGS84 +no_defs"
#                                                                                                                                   utm
#                                                                                     "+proj=utm +zone=30 +north +ellps=WGS84 +no_defs"
#                                                                                                                                  aeqd
#                                                   "+proj=aeqd +lon_0=-2.78488515753149 +lat_0=54.0378529685846 +ellps=WGS84 +no_defs"
#                                                                                                                                pconic
# "+proj=pconic +lon_0=-2.78488515753149 +lat_0=54.0378529685846 +lat_1=60.8433810004329 +lat_2=49.8848260001199 +ellps=WGS84 +no_defs"
#                                                                                                                                  eqdc
#                           "+proj=eqdc +lon_0=-2.78488515753149 +lat_1=60.8433810004329 +lat_2=49.8848260001199 +ellps=WGS84 +no_defs"

Fwiw eqdc has no lat_0 param🙏

Ah yes, thanks for that! 🙂 Updated now.

Useful function indeed !
I suggest to use sf::st_coordinates(sf::st_centroid(sf::st_geometry(x))) here to avoid the st_centroid assumes attributes are constant over geometries warning.

Useful function indeed ! I suggest to use sf::st_coordinates(sf::st_centroid(sf::st_geometry(x))) here to avoid the st_centroid assumes attributes are constant over geometries warning.

Great suggestion - thankyou! I've added that in now :)