1# alternative, but more limiting from sp/R/CRS-methods.R, https://github.com/edzer/sp/pull/31 @hughjonesd 2# (no longer used) 3#identicalCRS1 = function(x, y) { 4# args_x <- strsplit(x, " +")[[1]] 5# args_y <- strsplit(y, " +")[[1]] 6# setequal(args_x, args_y) 7#} 8 9# this function establishes whether two crs objects are semantically identical. This is 10# the case when: (1) they are completely identical (including NA), or (2) GDAL considers 11# them equivalent 12#' @export 13Ops.crs <- function(e1, e2) { 14 if (nargs() == 1) 15 stop(paste("unary", .Generic, "not defined for \"crs\" objects"), call. = FALSE) 16 17 cmp <- switch(.Generic, "==" =, "!=" = TRUE, FALSE) 18 if (!cmp) 19 stop(paste("operation", .Generic, "not supported for crs objects"), call. = FALSE) 20 if (.Generic == "!=") 21 !(e1 == e2) 22 else { # "==": check semantic equality 23 if (isTRUE(all.equal(e1, e2))) # includes both are NA_crs_ 24 TRUE 25 else if (is.na(e1) || is.na(e2)) # only one of them is NA_crs_ 26 FALSE 27 else 28 isTRUE(try(CPL_crs_equivalent(e1, e2), silent = TRUE)) # use GDAL's srs1->IsSame(srs2) 29 } 30} 31 32#' Retrieve coordinate reference system from object 33#' 34#' Retrieve coordinate reference system from sf or sfc object 35#' @name st_crs 36#' @param x numeric, character, or object of class \link{sf} or \link{sfc} 37#' @param ... ignored 38#' @export 39#' @return If \code{x} is numeric, return \code{crs} object for EPSG:\code{x}; 40#' if \code{x} is character, return \code{crs} object for \code{x}; 41#' if \code{x} is of class \code{sf} or \code{sfc}, return its \code{crs} object. 42#' @details The *crs functions create, get, set or replace the \code{crs} attribute 43#' of a simple feature geometry list-column. This attribute is of class \code{crs}, 44#' and is a list consisting of \code{input} (user input, e.g. "EPSG:4326" or "WGS84" 45#' or a proj4string), and \code{wkt}, an automatically generated wkt2 representation of the crs. 46#' If \code{x} is identical to the wkt2 representation, and the CRS has a name, this name 47#' is used for the \code{input} field. 48#' 49#' Comparison of two objects of class \code{crs} uses the GDAL function 50#' \code{OGRSpatialReference::IsSame}. 51#' @return Object of class \code{crs}, which is a list with elements \code{input} (length-1 character) 52#' and \code{wkt} (length-1 character). 53#' Elements may be \code{NA} valued; if all elements are \code{NA} the CRS is missing valued, and coordinates are 54#' assumed to relate to an arbitrary Cartesian coordinate system. 55st_crs = function(x, ...) UseMethod("st_crs") 56 57#' @name st_crs 58#' @export 59st_crs.sf = function(x, ...) st_crs(st_geometry(x), ...) 60 61#' @name st_crs 62#' @export 63st_crs.numeric = function(x, ...) { 64 make_crs(paste0("EPSG:", x)) 65} 66 67 68#' @name st_crs 69#' @export 70st_crs.character = function(x, ...) { 71 if (is.na(x)) 72 NA_crs_ 73 else { 74 crs = make_crs(x) 75 if (is.na(crs)) 76 stop(paste("invalid crs:", x)) 77 # if we input wkt2, and CRS has a name, use it: 78 if (identical(x, crs$wkt) && !identical(crs$Name, "unknown")) 79 crs$input = crs$Name 80 crs 81 } 82} 83 84fix_crs = function(x) { 85 if (all(c("epsg", "proj4string") %in% names(x))) { 86 message("old-style crs object detected; please recreate object with a recent sf::st_crs()") 87 x = unclass(x) 88 if (!is.na(x$epsg)) 89 st_crs(x$epsg) 90 else 91 st_crs(x$proj4string) 92 } else 93 x 94} 95 96 97#' @name st_crs 98#' @param parameters logical; \code{FALSE} by default; if \code{TRUE} return a list of coordinate reference system parameters, with named elements \code{SemiMajor}, \code{InvFlattening}, \code{units_gdal}, \code{IsVertical}, \code{WktPretty}, and \code{Wkt} 99#' @export 100st_crs.sfc = function(x, ..., parameters = FALSE) { 101 crs = fix_crs(attr(x, "crs")) 102 if (parameters) { 103 if (is.na(crs)) 104 list() 105 else 106 crs_parameters(crs) 107 } else 108 crs 109} 110 111#' @name st_crs 112#' @export 113st_crs.bbox = function(x, ...) { 114 crs = attr(x, "crs") 115 if (is.null(crs)) 116 NA_crs_ 117 else 118 crs 119} 120 121#' @name st_crs 122#' @export 123st_crs.CRS = function(x, ...) { 124 if (is.null(comment(x)) || CPL_proj_version() < "6.0.0" || 125 CPL_gdal_version() < "3.0.0") 126 st_crs(x@projargs) 127 else { 128 ret = st_crs(comment(x)) 129 name = ret$Name 130 ret$input = if (name == "unknown") 131 x@projargs 132 else 133 name 134 ret 135 } 136} 137 138#' @name st_crs 139#' @export 140st_crs.crs = function(x, ...) x 141 142#' @export 143st_crs.default = function(x, ...) NA_crs_ 144 145#' Set or replace coordinate reference system from object 146#' 147#' Set or replace retrieve coordinate reference system from object 148#' @name st_crs 149#' @param value one of (i) character: a string accepted by GDAL, (ii) integer, a valid EPSG value (numeric), or (iii) an object of class \code{crs}. 150#' @details In case a coordinate reference system is replaced, no transformation takes 151#' place and a warning is raised to stress this. 152#' 153#' @export 154`st_crs<-` = function(x, value) UseMethod("st_crs<-") 155 156#' @name st_crs 157#' @examples 158#' sfc = st_sfc(st_point(c(0,0)), st_point(c(1,1))) 159#' sf = st_sf(a = 1:2, geom = sfc) 160#' st_crs(sf) = 4326 161#' st_geometry(sf) 162#' @export 163`st_crs<-.sf` = function(x, value) { 164 st_crs(x[[ attr(x, "sf_column") ]]) = value 165 x 166} 167 168# return crs object from crs, integer, or character string 169make_crs = function(x) { 170 171 if (inherits(x, "CRS")) { 172 x = if (is.null(comment(x)) || (CPL_proj_version() < "6.0.0" || 173 CPL_gdal_version() < "3.0.0")) 174 175 x@projargs 176 else 177 comment(x) # WKT2 178 } 179 if (is.numeric(x) && !is.na(x)) 180 x = paste0("EPSG:", x) 181 # return: 182 if (is.na(x)) 183 NA_crs_ 184 else if (inherits(x, "crs")) 185 x 186 else if (is.character(x)) { 187 if (grepl("+init=epsg:", x) && 188 sf_extSoftVersion()[["proj.4"]] >= "6.0.0" && 189 sf_extSoftVersion()[["proj.4"]] < "6.3.1") { # nocov start FIXME: 190 x = strsplit(x, " ")[[1]] 191 if (length(x) > 1) 192 warning(paste("the following proj4string elements are ignored:", 193 paste(x[-1], collapse = " "), "; remove the +init=epsg:XXXX to undo this")) 194 x = paste0("EPSG:", as.integer(substr(x[1], 12, 20))) # nocov end 195 } 196 CPL_crs_from_input(x) 197 } else 198 stop(paste("cannot create a crs from an object of class", class(x)), call. = FALSE) 199} 200 201#' @name st_crs 202#' @examples 203#' sfc = st_sfc(st_point(c(0,0)), st_point(c(1,1))) 204#' st_crs(sfc) = 4326 205#' sfc 206#' @export 207`st_crs<-.sfc` = function(x, value) { 208 209 if (is.null(attr(x, "crs"))) 210 start_crs = NA_crs_ 211 else 212 start_crs = st_crs(x) 213 214 end_crs = make_crs(value) 215 216 if (!is.na(start_crs) && !is.na(end_crs) && start_crs != end_crs) 217 warning("st_crs<- : replacing crs does not reproject data; use st_transform for that", call. = FALSE) 218 219 structure(x, crs = end_crs) 220} 221 222#' @export 223`st_crs<-.bbox` = function(x, value) { 224 structure(x, crs = make_crs(value)) 225} 226 227 228#' @name st_crs 229#' @examples 230#' sfc = st_sfc(st_point(c(0,0)), st_point(c(1,1))) 231#' library(dplyr) 232#' x = sfc %>% st_set_crs(4326) %>% st_transform(3857) 233#' x 234#' @export 235st_set_crs = function(x, value) { 236 st_crs(x) = value 237 x 238} 239 240#' Assert whether simple feature coordinates are longlat degrees 241#' 242#' Assert whether simple feature coordinates are longlat degrees 243#' @param x object of class \link{sf} or \link{sfc}, or otherwise an object of a class that has an \link{st_crs} method returning a \code{crs} object 244#' @return TRUE if x has geographic coordinates, FALSE if it has projected coordinates, or NA if \code{is.na(st_crs(x))}. 245#' @export 246st_is_longlat = function(x) { 247 crs = st_crs(x) 248 if (is.na(crs)) 249 NA 250 else { 251 ret = crs_parameters(crs)$IsGeographic 252 if (ret && inherits(x, c("sf", "sfc", "stars"))) { 253 bb = st_bbox(x) 254 # check for potentially meaningless value range: 255 eps = sqrt(.Machine$double.eps) 256 if (all(!is.na(unclass(bb))) && 257 (bb["xmin"] < (-180-eps) || bb["xmax"] > (360+eps) || bb["ymin"] < (-90-eps) || bb["ymax"] > (90+eps))) 258 warning("bounding box has potentially an invalid value range for longlat data") 259 } 260 ret 261 } 262} 263 264# a = "b" => a is the proj.4 unit (try: cs2cs -lu); "b" is the udunits2 unit 265udunits_from_proj = list( 266# PROJ.4 UDUNITS 267 `km` = as_units("km"), 268 `m` = as_units("m"), 269 `dm` = as_units("dm"), 270 `cm` = as_units("cm"), 271 `mm` = as_units("mm"), 272 `kmi` = as_units("nautical_mile"), 273 `in` = as_units("in"), 274 `ft` = as_units("ft"), 275 `yd` = as_units("yd"), 276 `mi` = as_units("mi"), 277 `fath` = as_units("fathom"), 278 `ch` = as_units("chain"), 279 `link` = as_units("link", check_is_valid = FALSE), # not (yet) existing; set in .onLoad() 280 `us-in` = as_units("us_in", check_is_valid = FALSE), 281 `us-ft` = as_units("US_survey_foot"), 282 `us-yd` = as_units("US_survey_yard"), 283 `us-ch` = as_units("chain"), 284 `us-mi` = as_units("US_survey_mile"), 285 `ind-yd` = as_units("ind_yd", check_is_valid = FALSE), 286 `ind-ft` = as_units("ind_ft", check_is_valid = FALSE), 287 `ind-ch` = as_units("ind_ch", check_is_valid = FALSE) 288) 289 290crs_parameters = function(x, with_units = TRUE) { 291 stopifnot(inherits(x, "crs")) 292 if(is.na(x)) return(list(NA)) 293 294 ret = CPL_crs_parameters(x) 295 units(ret$SemiMajor) = as_units("m") 296 units(ret$SemiMinor) = as_units("m") 297 if (with_units) 298 ret$ud_unit = if (isTRUE(ret$IsGeographic)) 299 as_units("arc_degree") # FIXME: is this always true? 300 else if (is.null(x$units)) 301 as_units("m") 302 else if (is.character(udunits_from_proj[[x$units]])) 303 as_units(udunits_from_proj[[x$units]]) 304 else 305 udunits_from_proj[[x$units]] 306 ret 307} 308 309epsg = function(x) { 310 if (is.na(x)) 311 NA_integer_ 312 else if (grepl("^EPSG:", x[["input"]])) 313 as.integer(gsub("^EPSG:(\\d+)\\b.*$", "\\1", x[["input"]])) 314 else 315 crs_parameters(x, with_units = FALSE)[["epsg"]] 316} 317 318proj4string = function(x) { 319 if (is.na(x)) 320 NA_character_ 321 else 322 crs_parameters(x, with_units = FALSE)[["proj4string"]] 323} 324 325 326#' @name st_as_text 327#' @param projjson logical; if TRUE, return projjson form (requires GDAL 3.1 and PROJ 6.2), else return well-known-text form 328#' @param pretty logical; if TRUE, print human-readable well-known-text representation of a coordinate reference system 329#' @export 330st_as_text.crs = function(x, ..., projjson = FALSE, pretty = FALSE) { 331 if (is.na(x)) 332 NA_character_ 333 else if (projjson) { 334 if (sf_extSoftVersion()["GDAL"] < "3.1.0" || sf_extSoftVersion()["proj.4"] < "6.2.0") 335 stop("ProjJson requires GDAL >= 3.1.0 and PROJ >= 6.2.0") 336 crs_parameters(x)$ProjJson 337 } else { # wkt: 338 if (pretty) 339 crs_parameters(x)$WktPretty 340 else 341 crs_parameters(x)$Wkt 342 } 343} 344 345 346#' @name st_crs 347#' @details 348#' \code{NA_crs_} is the \code{crs} object with missing values for \code{input} and \code{wkt}. 349#' @export 350NA_crs_ = structure( 351 list(input = NA_character_, 352 wkt = NA_character_), 353 class = "crs") 354 355#' @name st_crs 356#' @export 357#' @method is.na crs 358is.na.crs = function(x) { 359 identical(x, NA_crs_) 360} 361 362#' @name st_crs 363#' @param name element name 364#' @export 365#' @examples 366#' st_crs("EPSG:3857")$input 367#' st_crs(3857)$proj4string 368#' st_crs(3857)$b # numeric 369#' st_crs(3857)$units # character 370#' @details the \code{$} method for \code{crs} objects retrieves named elements 371#' using the GDAL interface; named elements include 372#' \code{"SemiMajor"}, \code{"SemiMinor"}, \code{"InvFlattening"}, \code{"IsGeographic"}, 373#' \code{"units_gdal"}, \code{"IsVertical"}, \code{"WktPretty"}, \code{"Wkt"}, 374#' \code{"Name"}, \code{"proj4string"}, \code{"epsg"}, \code{"yx"} and 375#' \code{"ud_unit"} (this may be subject to changes in future GDAL versions). 376#' @export 377`$.crs` = function(x, name) { 378 379 if (!is.null(x[["proj4string"]])) { # old-style object: 380 warning("CRS uses proj4string, which is deprecated.") 381 x = st_crs(x[["proj4string"]]) # FIXME: should this be only for some transition period? Add test? 382 } 383 if (is.na(x)) 384 NA_character_ 385 else if (is.numeric(name) || name %in% names(x)) 386 x[[name]] 387 else { 388 p = crs_parameters(x, with_units = FALSE) 389 if (name %in% names(p)) 390 p[[name]] 391 else { 392 tryNum = function(x) { n = suppressWarnings(as.numeric(x)); if (is.na(n)) x else n } 393 p4s = strsplit(p$proj4string, " ")[[1]] 394 p4s2 = strsplit(p4s, "=") 395 vals = lapply(p4s2, function(x) if (length(x) == 1) TRUE else tryNum(x[2])) 396 names(vals) = substring(sapply(p4s2, function(x) x[1]), 2) 397 vals[[name]] 398 } 399 } 400} 401 402#' @export 403print.crs = function(x, ...) { 404 cat("Coordinate Reference System:") 405 if (is.na(x)) { 406 cat(" NA\n") 407 } else { 408 cat("\n") 409 if (is.na(x$input)) 410 cat(" No user input\n") 411 else 412 cat(" User input:", x$input, "\n") 413 414 # print wkt: 415 if (!is.na(x$wkt)) 416 cat(" wkt:\n", x$wkt, "\n", sep = "") 417 } 418} 419 420#' @name st_crs 421#' @export 422#' @details format.crs returns NA if the crs is missing valued, or else 423#' the name of a crs if it is different from "unknown", or 424#' else the user input if it was set, or else its "proj4string" representation; 425format.crs = function(x, ...) { 426 if (is.na(x)) 427 NA_character_ 428 else { 429 p = crs_parameters(x) 430 if (p$Name == "unknown") { 431 if (x$input == "unknown") 432 x$proj4string 433 else 434 x$input 435 } else 436 x$Name 437 } 438} 439 440 441#' @export 442st_crs.Raster = function(x, ...) { 443 crsobj <- raster::crs(x) 444 st_crs(crsobj) # nocov 445} 446 447#' @export 448st_crs.Spatial = function(x, ...) { 449 if (! requireNamespace("sp", quietly = TRUE)) 450 stop("package sp required, please install it first") 451 st_crs(x@proj4string) # nocov 452} 453 454#' @name st_crs 455#' @param authority_compliant logical; specify whether axis order should be 456#' handled compliant to the authority; if omitted, the current value is printed. 457#' @details 458#' \code{st_axis_order} can be used to get and set the axis order: \code{TRUE} 459#' indicates axes order according to the authority 460#' (e.g. EPSG:4326 defining coordinates to be latitude,longitude pairs), \code{FALSE} 461#' indicates the usual GIS (display) order (longitude,latitude). This can be useful 462#' when data are read, or have to be written, with coordinates in authority compliant order. 463#' The return value is the current state of this (\code{FALSE}, by default). 464#' @return \code{st_axis_order} returns the (logical) current value if called without 465#' argument, or (invisibly) the previous value if it is being set. 466#' @export 467#' @examples 468#' pt = st_sfc(st_point(c(0, 60)), crs = 4326) 469#' # st_axis_order() only has effect in GDAL >= 2.5.0: 470#' st_axis_order() # query default: FALSE means interpret pt as (longitude latitude) 471#' st_transform(pt, 3857)[[1]] 472#' old_value = FALSE 473#' if (sf_extSoftVersion()["GDAL"] >= "2.5.0") 474#' (old_value = st_axis_order(TRUE)) 475#' # now interpret pt as (latitude longitude), as EPSG:4326 prescribes: 476#' st_axis_order() # query current value 477#' st_transform(pt, 3857)[[1]] 478#' st_axis_order(old_value) # set back to old value 479st_axis_order = function(authority_compliant = logical(0)) { 480 ret = CPL_axis_order_authority_compliant(authority_compliant) 481 if (length(authority_compliant)) 482 invisible(ret) 483 else 484 ret 485} 486