R/st_crs2.R

In ggranga/fidolasen: Find, Download and Process Sentinel-2 Data

#' @title Retrieve coordinate reference system from sf or sfc object
#' @name st_crs2
#' @description This function is a wrapper for [sf::st_crs], unless
#'  threating numeric `character` strings as integers, and
#'  accepting also UTM timezones, paths of spatial files and paths of 
#'  text files containing WKT like .prj (see details) .
#' @param x numeric, character, or object of class \link{sf} or \link{sfc}, being:
#'  - EPSG code: numeric (e.g. `32632`) or character (in the form
#'      `"32632"` or `"EPSG:32632"`);
#'  - UTM zone: numeric (e.g. `32`, interpreted as 32 North) or character
#'      (e.g. `"32"` or `"32N"` for zone 32 North, `"32S"` for 32 South);
#'  - WKT test: passed as character string or as path of a text file 
#'      containing it (e.g. the path of a .prj file);
#'  - PROJ.4 string, passed as character (e.g. 
#'      `"+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"` 
#'      (**NOTE**: this representation is deprecated with PROJ >= 6
#'      -- see http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html --
#'      so a warning is returned using it, unless the string contains only
#'      the epsg code -- e.g. `"+init=epsg:32632"`, in which case the EPSG
#'      code is taken);
#'  - path of a spatial file (managed by [sf::st_read] or [stars::read_stars]),
#'      passed as character string of length 1;
#'  - spatial file of class \link{sf} or \link{sfc}.
#' @param ... other parameters passed to [sf::st_crs].
#' @return An object of class \link{crs} of length 2.
#' @details See [sf::st_crs] for details.
#' @importFrom sf gdal_crs st_crs
#' @export
#' @author Luigi Ranghetti, phD (2019)
#' @references L. Ranghetti, M. Boschetti, F. Nutini, L. Busetto (2020).
#'  "sen2r": An R toolbox for automatically downloading and preprocessing 
#'  Sentinel-2 satellite data. _Computers & Geosciences_, 139, 104473. 
#'  \doi{10.1016/j.cageo.2020.104473}, URL: \url{https://sen2r.ranghetti.info/}.
#' @note License: GPL 3.0
#' @examples 
#' ## CRS from EPSG
#' st_crs2(32609)
#' st_crs2("EPSG:32609")
#' 
#' ## CRS from UTM zone
#' st_crs2(9)
#' st_crs2("09")
#' st_crs2("9N")
#' st_crs2("09S")
#' 
#' ## CRS from WKT (string or path)
#' (wkt_32n <- sf::st_as_text(sf::st_crs(32609)))
#' st_crs2(wkt_32n)
#' writeLines(wkt_32n, wkt_32n_path <- tempfile())
#' st_crs2(wkt_32n_path)
#' 
#' \dontrun{
#' ## CRS from spatial file path
#' raster_path <- system.file(
#'   "extdata/out/S2A2A_20190723_022_Barbellino_BOA_10.tif", 
#'   package="sen2r"
#' )
#' vector_path <- system.file(
#'   "extdata/vector/barbellino.geojson", 
#'   package="sen2r"
#' )
#' try( st_crs2(raster_path) )
#' st_crs2(vector_path)
#' 
#' ## CRS from spatial files
#' st_crs2(stars::read_stars(raster_path))
#' st_crs2(sf::read_sf(vector_path))
#' 
#' ## CRS from PROJ.4 string
#' # (avoid using this with PROJ >= 6!)
#' st_crs2("+init=epsg:32609") # this makes use of the EPSG code
#' st_crs2("+proj=utm +zone=9 +datum=WGS84 +units=m +no_defs")
#' st_crs2(raster::raster(raster_path)) # st_crs(raster) uses the PROJ.4 as input
#' }


st_crs2 <- function(x, ...) UseMethod("st_crs2")

## character: several cases (see)
#' @export
st_crs2.character <- function(x, ...) {
  
  ## case 1: EPSG code / UTM zone
  x_epsg <- if (grepl("^(([0-5]?[0-9])|60)[Nn]?$", x)) {
    # x: UTM zone North -> integer EPSG
    as.integer(paste0(
      "326", 
      str_pad2(
        gsub("^(([0-5]?[0-9])|60)[Nn]?$", "\\1", x),
        2, "left", "0"
      )
    ))
  } else if (grepl("^(([0-5]?[0-9])|60)[Ss]$", x)) {
    # x: UTM zone South -> integer EPSG
    as.integer(paste0(
      "327", 
      str_pad2(
        gsub("^(([0-5]?[0-9])|60)[Ss]$", "\\1", x),
        2, "left", "0"
      )
    ))
  } else if (grepl("^[0-9]+$", x)) {
    # x: EPSG (integer, numeric or character) -> integer EPSG
    as.integer(x)
  } else if (grepl("^[Ee][Pp][Ss][Gg]\\:[0-9]+$", x)) {
    # x: EPSG (in the form "EPSG:xxx") -> integer EPSG
    as.integer(gsub("^[Ee][Pp][Ss][Gg]\\:([0-9]+)$", "\\1", x))
  } else if (grepl("^\\+init\\=epsg:[0-9]+$", tolower(x))) {
    # x: PROJ.4 with only EPSG -> integer EPSG
    as.integer(gsub("^\\+init\\=epsg:([0-9]+)$", "\\1", tolower(x)))
  } else {
    NULL
  }
  if (!is.null(x_epsg)) {
    return(sf::st_crs(x_epsg, ...))
  }
  
  ## case 2: PROJ.4
  if (grepl("^\\+[a-z]+\\=", x)) {
    # x: PROJ.4 -> character PROJ.4 with warning
    print_message(
      type = "warning",
      "Using PROJ.4 strings is deprecated with PROJ >= 6 ",
      "(see https://www.r-spatial.org/r/2020/03/17/wkt.html)."
    )
    return(sf::st_crs(x, ...))
  }
  
  ## case 3: file path
  if (file.exists(as.character(x))) {
    # x: file path -> spatial file or WKT
    x2 <- tryCatch(
      # x: path of a vector file -> sf
      st_read(x, quiet = TRUE),
      warning = function(w) {
        if (grepl("no simple feature geometries present\\: returning a data\\.frame", w)) {
          # x: path of a tabular file -> x (st_crs will return the proper error)
          x
        } else {st_read(x, quiet = TRUE)}
      },
      error = function(e) {tryCatch(
        # x: path of a text file with WKT -> crs
        suppressWarnings(sf::st_crs(paste(readLines(x), collapse=""))),
        error = function(e) {tryCatch(
          # x: path of a raster file -> stars proxy
          gdal_crs(x),
          error = function(e) {
            # x: path of a non supported file -> x (st_crs will return the proper error)
            x
          }
        )}
      )}
    )
    return(sf::st_crs(x2, ...))
  }
  
  ## case 4: WKT and other characters
  if (grepl("^((PROJCR?S)|(GEOGCR?S))\\[.+\\]$", x)) {
    # x: WKT string -> crs
    return(sf::st_crs(x, ...))
  }
  
  ## any other case: pass to st_crs as is
  sf::st_crs(x, ...)
  
}

## integer or numeric (EPGS / UTM zone): threat as character
#' @export
st_crs2.integer <- function(x, ...) {st_crs2.character(as.character(x), ...)}
#' @export
st_crs2.numeric <- function(x, ...) {st_crs2.character(as.character(x), ...)}

## classes already managed by st_crs()
#' @export
st_crs2.default <- function(x, ...) {
  if (missing(x)) {sf::st_crs(NA)} else {sf::st_crs(x, ...)}
}