1# Copyright 2006-2021 Roger Bivand 2 3readOGR <- function(dsn, layer, verbose=TRUE, p4s=NULL, 4 stringsAsFactors=as.logical(NA), 5 drop_unsupported_fields=FALSE, 6 pointDropZ=FALSE, dropNULLGeometries=TRUE, useC=TRUE, 7 disambiguateFIDs=FALSE, addCommentsToPolygons=TRUE, encoding=NULL, 8 use_iconv=FALSE, swapAxisOrder=FALSE, require_geomType=NULL, 9 integer64="no.loss", GDAL1_integer64_policy=FALSE, 10 morphFromESRI=NULL, dumpSRS=FALSE, enforce_xy=NULL, 11 D3_if_2D3D_points=FALSE, missing_3D=0) { 12 if (missing(dsn)) stop("missing dsn") 13 stopifnot(is.character(dsn)) 14 stopifnot(length(dsn) == 1L) 15# copy sf::st_read.default usage 16 if (length(dsn) == 1 && file.exists(dsn)) 17 dsn <- enc2utf8(normalizePath(dsn)) 18 if (nchar(dsn) == 0) stop("empty name") 19 if (missing(layer)){ 20 layers <- ogrListLayers(dsn=dsn) 21 if (length(layers) == 0L) stop("missing layer") 22 if (length(layers) > 0L) layer <- c(layers[1]) 23 if (length(layers) > 1L) 24 warning("First layer ", layer, 25 " read; multiple layers present in\n", dsn, 26 ", check layers with ogrListLayers()") 27 } else layer <- enc2utf8(layer) 28# stop("missing dsn") 29# if (missing(layer)) stop("missing layer") 30 if (nchar(layer) == 0) stop("empty name") 31 integer64 <- match.arg(integer64, 32 c("allow.loss", "warn.loss", "no.loss")) 33 34 int64 <- switch(integer64, 35 "allow.loss"=1L, 36 "warn.loss"=2L, 37 "no.loss"=3L) 38 if (GDAL1_integer64_policy) int64 <- 4L 39# adding argument for SHAPE_ENCODING environment variable 121124 40 stopifnot(is.logical(use_iconv)) 41 stopifnot(length(use_iconv) == 1) 42 if (!is.null(encoding)) { 43 stopifnot(is.character(encoding)) 44 stopifnot(length(encoding) == 1) 45 } 46 WKB <- c("wkbPoint", "wkbLineString", "wkbPolygon", "wkbMultiPoint", 47 "wkbMultiLineString", "wkbMultiPolygon", "wkbGeometryCollection") 48 if (!is.null(require_geomType)) { 49 stopifnot(is.character(require_geomType) && 50 length(require_geomType)==1) 51 m_require_geomType <- match(require_geomType, WKB) 52 stopifnot(!is.na(m_require_geomType) || m_require_geomType <= 3) 53 } 54 if(is.na(stringsAsFactors)) { 55 stringsAsFactors <- 56 if(getRversion() < "4.1.0") 57 default.stringsAsFactors() 58 else 59 FALSE 60 } 61 62 suppressMessages(ogr_info <- ogrInfo(dsn=dsn, layer=layer, 63 encoding=encoding, use_iconv=use_iconv, 64 swapAxisOrder=swapAxisOrder, require_geomType=require_geomType, 65 morphFromESRI=morphFromESRI, dumpSRS=dumpSRS, 66 D3_if_2D3D_points=D3_if_2D3D_points)) 67 HAS_FEATURES <- TRUE 68 if (!ogr_info$have_features) { 69 if (dropNULLGeometries) { 70 stop("no features found") 71 } else { 72 warning("no features found; proceeding to atttributes only") 73 HAS_FEATURES <- FALSE 74 } 75 } 76 if (is.null(ogr_info$nListFields)) nListFields <- 0 77 else nListFields <- ogr_info$nListFields 78# 121130 RSB trap no field case (from PostGIS, Mathieu Basille) 79 if (ogr_info$nitems > 0) { 80 nodata_flag <- FALSE 81 if (strsplit(getGDALVersionInfo(), " ")[[1]][2] < "2") { 82 keep <- ogr_info$iteminfo$typeName %in% c("Integer", "Real", 83 "String", "Date", "Time", "DateTime", "IntegerList", 84 "RealList", "StringList") 85 } else { 86 keep <- ogr_info$iteminfo$typeName %in% c("Integer", "Real", 87 "String", "Date", "Time", "DateTime", "IntegerList", 88 "RealList", "StringList", "Integer64", "Integer64List") 89 } 90 if (nListFields > 0) 91 ListFields <- as.integer(ogr_info$iteminfo$maxListCount) 92 if (drop_unsupported_fields) { 93 iflds <- as.integer((1:ogr_info$nitems)-1) 94 iflds <- iflds[keep] 95 fldnms <- ogr_info$iteminfo$name[keep] 96 if (nListFields > 0) ListFields <- ListFields[keep] 97 if (any(!keep)) warning(paste("Fields dropped:", 98 paste(ogr_info$iteminfo$name[!keep], collapse=" "))) 99 } else { 100 if (any(!keep)) stop(paste("Unsupported field type:", 101 paste(ogr_info$iteminfo$typeName[!keep], collapse=" "))) 102 iflds <- as.integer((1:ogr_info$nitems)-1) 103 fldnms <- ogr_info$iteminfo$name 104 } 105 int64_found <- grep("Integer64", ogr_info$iteminfo$typeName) 106 } else { 107 int64_found <- integer(0L) 108 nodata_flag <- TRUE 109 iflds <- integer(0) 110 } 111 fids <- ogrFIDs(dsn=dsn, layer=layer) 112 if (attr(fids, "i") != attr(fids, "nf")) { 113 retain <- 1:attr(fids, "i") 114 afids <- 0:(attr(fids, "nf")-1) 115 deleted <- afids[!(afids %in% fids[retain])] 116 warning(paste("Deleted feature IDs: ", paste(deleted, 117 collapse=", "))) 118 fids <- fids[retain] 119 } else { 120 retain <- NULL 121 } 122 attributes(fids) <- NULL 123# suggestion by Dylan Beaudette 110620 124 non_unique_fids <- max(table(fids)) > 1 125 if (non_unique_fids) { 126 if (disambiguateFIDs) { 127 fids <- seq_along(fids) # if present, make new FIDs 128 } else { 129 stop("FIDs not unique") 130 } 131 } 132 if (verbose) { 133 cat("OGR data source with driver:", ogr_info$driver, "\n") 134 cat("Source: \"", dsn, '\", layer: \"', layer, "\"", '\n', 135 sep='') 136 cat("with", length(fids), "features") 137 if (!is.null(attr(ogr_info, "require_geomType"))) 138 cat(";\nSelected", attr(ogr_info, "require_geomType"), 139 "feature type, with", sum(attr(ogr_info, "keepGeoms")), 140 "rows") 141 cat("\n") 142 cat("It has", length(iflds), "fields") 143 if (nListFields > 0) 144 cat(", of which", nListFields, "list fields") 145 cat("\n") 146 if (length(int64_found > 0L)) { 147 if (GDAL1_integer64_policy) { 148 cat("Integer64 fields read as doubles: ", 149 paste(fldnms[int64_found], collapse=" "), "\n") 150 } else { 151 if (integer64 == "no.loss") { 152 cat("Integer64 fields read as strings: ", 153 paste(fldnms[int64_found], collapse=" "), "\n") 154 } else { 155 cat("Integer64 fields read as signed 32-bit integers: ", 156 paste(fldnms[int64_found], collapse=" "), "\n") 157 } 158 } 159 } 160 161 } 162 163# adding argument for SHAPE_ENCODING environment variable 121124 164 if (!use_iconv && !is.null(encoding) && 165 ogr_info$driver == "ESRI Shapefile") { 166 oSE <- getCPLConfigOption("SHAPE_ENCODING") 167 tull <- setCPLConfigOption("SHAPE_ENCODING", encoding) 168 } 169 if (nodata_flag) { 170 dlist <- list(FID=as.integer(fids)) 171 } else { 172 attr(iflds, "nListFields") <- as.integer(nListFields) 173 nflds <- length(iflds) 174 if (nListFields > 0) { 175 attr(iflds, "ListFields") <- ListFields 176 nflds <- nflds + sum(ListFields) - nListFields 177 fldnms1 <- NULL 178 for (i in seq(along=ListFields)) { 179 if (ListFields[i] == 0) fldnms1 <- c(fldnms1, fldnms[i]) 180 else fldnms1 <- c(fldnms1, 181 paste(fldnms[i], 1:ListFields[i], sep="")) 182 } 183 stopifnot(length(fldnms1) == nflds) 184 fldnms <- fldnms1 185 } 186 attr(iflds, "nflds") <- as.integer(nflds) 187 attr(iflds, "int64") <- as.integer(int64) 188 dlist <- .Call("ogrDataFrame", as.character(dsn), 189 enc2utf8(as.character(layer)), as.integer(fids), iflds, PACKAGE="rgdal") 190 names(dlist) <- make.names(fldnms ,unique=TRUE) 191 192 if (use_iconv && !is.null(encoding)) { 193 for (i in seq(along=dlist)) { 194 if (is.character(dlist[[i]])) { 195 dlist[[i]] <- iconv(dlist[[i]], from=encoding) 196 } 197 } 198 } 199 } 200 if (!use_iconv && !is.null(encoding) && 201 ogr_info$driver == "ESRI Shapefile") { 202 tull <- setCPLConfigOption("SHAPE_ENCODING", oSE) 203 } 204 data <- data.frame(dlist, row.names=fids, 205 stringsAsFactors=stringsAsFactors) 206 rm(dlist) 207 gc(verbose = FALSE) 208 209 if (!HAS_FEATURES) { 210 return(data) 211 } 212 213# suggestion by Paul Hiemstra 070817 214# if (is.null(morphFromESRI)) { 215# if (ogr_info$driver == "ESRI Shapefile") morphFromESRI <- TRUE 216# else morphFromESRI <- FALSE 217# } 218# if (ogr_info$driver != "ESRI Shapefile" && morphFromESRI) 219# morphFromESRI <- FALSE 220# stopifnot(is.logical(morphFromESRI)) 221# stopifnot(length(morphFromESRI) == 1) 222# stopifnot(is.logical(dumpSRS)) 223# stopifnot(length(dumpSRS) == 1) 224# prj <- .Call("ogrP4S", as.character(dsn), enc2utf8(as.character(layer)), as.logical(morphFromESRI), as.logical(dumpSRS), PACKAGE="rgdal") 225 226# prj <- OGRSpatialRef(dsn=dsn, layer=layer, morphFromESRI=morphFromESRI, 227# dumpSRS=dumpSRS, driver=ogr_info$driver, enforce_xy=enforce_xy) 228 229 prj <- ogr_info$p4s 230 231 if (!is.null(p4s)) { 232 if (!is.na(prj)) { 233 warning("p4s= argument given as: ", p4s, "\n and read as: ", 234 c(prj), "\n read string overridden by given p4s= argument value") 235 } 236 } else { 237 p4s <- prj 238 } 239 240 if (!is.na(p4s) && nchar(p4s) == 0) p4s <- as.character(NA) 241# if (new_proj_and_gdal()) wkt2 <- ogr_info$wkt2 242 oCRS <- new("CRS", projargs=p4s) 243 if (new_proj_and_gdal()) comment(oCRS) <- ogr_info$wkt2 244 245 geometry <- .Call("R_OGR_CAPI_features", as.character(dsn), 246 enc2utf8(as.character(layer)), comments=addCommentsToPolygons, 247 PACKAGE="rgdal") 248 if (is.null(retain)) { 249 eType <- geometry[[4]] 250 with_z <- geometry[[6]] 251 isNULL <- as.logical(geometry[[7]]) 252 gFeatures <- geometry[[5]] 253 gComments <- geometry[[8]] 254 } else { 255 eType <- geometry[[4]][retain] 256 with_z <- geometry[[6]][retain] 257 isNULL <- as.logical(geometry[[7]])[retain] 258 gFeatures <- geometry[[5]][retain] 259 gComments <- geometry[[8]][retain] 260 } 261 rm(geometry); 262 gc(verbose = FALSE) 263 if (any(isNULL)) { 264 eType <- eType[!isNULL] 265 with_z <- with_z[!isNULL] 266 } 267 elevated <- attr(ogr_info$with_z, "elevated") 268 u_with_z <- unique(sort(with_z)) 269 if (length(u_with_z) != 1L) { 270 if (elevated) { 271 if (pointDropZ) { 272 warning("Dropping mixed third dimension") 273 u_with_z <- 0L 274 } else { 275 warning("Setting any missing third dimension to ", missing_3D) 276 u_with_z <- 1L 277 } 278 } else { 279 stop(paste("Multiple # dimensions:", paste((u_with_z + 2), 280 collapse=":"))) 281 } 282 } 283 if (u_with_z < 0 || u_with_z > 1) stop( 284 paste("Invalid # dimensions:", (u_with_z + 2))) 285 286 eType[eType == 5L] <- 2L 287 eType[eType == 6L] <- 3L 288 289 u_eType <- unique(sort(eType)) 290 291 t_eType <- table(eType) 292 if (is.null(require_geomType)) { 293 keepGeoms <- NULL 294 if (length(u_eType) > 1L) stop( 295 paste("Multiple incompatible geometries:", 296 paste(paste(WKB[as.integer(names(t_eType))], t_eType, sep=": "), 297 collapse="; "))) 298# if (length(u_eType) == 2L) { 299# if (u_eType[1] == 2 && u_eType[2] == 5) u_eType = 2 300# else if (u_eType[1] == 3 && u_eType[2] == 6) u_eType = 3 301# else stop(paste("Multiple incompatible geometries:", 302# paste(paste(WKB[as.integer(names(t_eType))], t_eType, sep=": "), 303# collapse="; "))) 304# } 305 } else { 306 if (!require_geomType %in% WKB[as.integer(names(t_eType))]) 307 stop(require_geomType, "not in", WKB[as.integer(names(t_eType))]) 308 u_eType <- match(require_geomType, WKB) 309 keepGeoms <- WKB[eType] == require_geomType 310 message("NOTE: keeping only ", sum(keepGeoms), " ", require_geomType, 311 " of ", length(keepGeoms), " features\n") 312 } 313 314 if (length(gFeatures) != length(fids)) stop("Feature mismatch") 315 316 if (any(isNULL)) { 317 if (dropNULLGeometries) { 318 warning(paste("Dropping null geometries:", paste(which(isNULL), 319 collapse=", "))) 320 gFeatures <- gFeatures[!isNULL] 321 data <- data[!isNULL, , drop=FALSE] 322 fids <- fids[!isNULL] 323 gComments <- gComments[!isNULL] 324 } else { 325 warning(paste("Null geometries found:", paste(which(isNULL), 326 collapse=", "))) 327 warning("dropNULLGeometries FALSE, returning only data for null-geometry features") 328 return(data[isNULL, , drop=FALSE]) 329 } 330 } 331 332 if (!is.null(require_geomType)) { 333 gFeatures <- gFeatures[keepGeoms] 334 data <- data[keepGeoms, , drop=FALSE] 335 fids <- fids[keepGeoms] 336 gComments <- gComments[keepGeoms] 337 } 338 339 if (u_eType == 1) { # points 340 if (u_with_z == 0 || pointDropZ) { 341 if (swapAxisOrder) { 342 coords <- do.call("rbind", lapply(gFeatures, 343 function(x) c(x[[1]][[2]], x[[1]][[1]]))) 344 } else { 345 coords <- do.call("rbind", lapply(gFeatures, 346 function(x) c(x[[1]][[1]], x[[1]][[2]]))) 347 } 348 } else { 349 if (elevated) { 350 if (swapAxisOrder) { 351 coords <- do.call("rbind", lapply(gFeatures, 352 function(x) c(x[[1]][[2]], x[[1]][[1]], 353 ifelse(length(x[[1]]) == 2L, missing_3D, x[[1]][[3]])))) 354 } else { 355 coords <- do.call("rbind", lapply(gFeatures, 356 function(x) c(x[[1]][[1]], x[[1]][[2]], 357 ifelse(length(x[[1]]) == 2L, missing_3D, x[[1]][[3]])))) 358 } 359 } else { 360 if (swapAxisOrder) { 361 coords <- do.call("rbind", lapply(gFeatures, 362 function(x) c(x[[1]][[2]], x[[1]][[1]], 363 ifelse(length(x[[1]]) == 2L, missing_3D, x[[1]][[3]])))) 364 } else { 365 coords <- do.call("rbind", lapply(gFeatures, 366 function(x) c(x[[1]][[1]], x[[1]][[2]], 367 ifelse(length(x[[1]]) == 2L, missing_3D, x[[1]][[3]])))) 368 } 369 } 370 } 371# data <- data.frame(dlist) 372 row.names(data) <- NULL 373 res <- SpatialPointsDataFrame(coords=coords, data=data, 374 proj4string=oCRS) 375 } else if (u_eType == 2) { # lines 376 if (u_with_z != 0) warning("Z-dimension discarded") 377 n <- length(gFeatures) 378 lnList <- vector(mode="list", length=n) 379 for (i in 1:n) { 380 iG <- gFeatures[[i]] 381 m <- length(iG) 382 lnlist <- vector(mode="list", length=m) 383 for (j in 1:m) { 384 jG <- iG[[j]] 385 if (swapAxisOrder) { 386 lnlist[[j]] <- Line(cbind(jG[[2]], jG[[1]])) 387 } else { 388 lnlist[[j]] <- Line(cbind(jG[[1]], jG[[2]])) 389 } 390 } 391 lnList[[i]] <- Lines(lnlist, ID=as.character(fids[i])) 392 } 393 SL <- SpatialLines(lnList, proj4string=oCRS) 394# data <- data.frame(dlist, row.names=fids) 395 res <- SpatialLinesDataFrame(SL, data) 396 } else if (u_eType == 3) { # polygons 397 if (u_with_z != 0) warning("Z-dimension discarded") 398 if (useC) { 399# plList <- .Call("make_polygonslist", gFeatures, 400# as.list(as.character(fids)), PACKAGE="rgdal") 401 n <- length(gFeatures) 402 plList <- vector(mode="list", length=n) 403 for (i in 1:n) { 404 iG <- gFeatures[[i]] 405 if (swapAxisOrder) { 406 iG <- lapply(iG, function(x) { 407 tmp <- x[[1]]; x[[1]] <- x[[2]]; x[[2]] <- tmp 408 }) 409 } 410 if (addCommentsToPolygons) { 411 thisPL <- Polygons(.Call("make_Polygonlist", 412 iG, gComments[[i]], PACKAGE="rgdal"), 413 ID=as.character(fids[i])) 414 comment(thisPL) <- paste(gComments[[i]], 415 collapse=" ") 416 } else { 417 thisPL <- Polygons(.Call("make_Polygonlist", 418 iG, NULL, PACKAGE="rgdal"), 419 ID=as.character(fids[i])) 420 } 421 plList[[i]] <- thisPL 422 } 423 } else { 424 n <- length(gFeatures) 425 plList <- vector(mode="list", length=n) 426 for (i in 1:n) { 427 iG <- gFeatures[[i]] 428 if (swapAxisOrder) { 429 iG <- lapply(iG, function(x) { 430 tmp <- x[[1]]; x[[1]] <- x[[2]]; x[[2]] <- tmp 431 }) 432 } 433 m <- length(iG) 434 pllist <- vector(mode="list", length=m) 435 for (j in 1:m) { 436 jG <- iG[[j]] 437 cmat <- cbind(jG[[1]], jG[[2]]) 438 if (!identical(cmat[1,], cmat[nrow(cmat),])) { 439 cmat <- rbind(cmat, cmat[1,]) 440 warning(paste("Ring closed in Polygons", 441 i, "Polygon", j)) 442 } 443 t0 <- try(pllist[[j]] <- Polygon(cmat), 444 silent=TRUE) 445 if (inherits(t0, "try-error")) { 446 print(cmat) 447 print(t0) 448 stop("i: ", i, ", j: ", j, 449 ", Polygon error exit") 450 } 451 } 452 thisPL <- Polygons(pllist, ID=as.character(fids[i])) 453 if (addCommentsToPolygons) { 454 comment(thisPL) <- paste(gComments[[i]], 455 collapse=" ") 456 if (!isTRUE(all.equal(as.logical(gComments[[i]]), 457 sapply(slot(thisPL, "Polygons"), slot, "hole")))) 458 warning("comment/hole mismatch, geometry:", i) 459 } 460 plList[[i]] <- thisPL 461 } 462 } 463 rm(gFeatures) 464 gc(verbose = FALSE) 465 SP <- SpatialPolygons(plList, proj4string=oCRS) 466 rm(plList) 467 gc(verbose = FALSE) 468# data <- data.frame(dlist, row.names=fids) 469 res <- SpatialPolygonsDataFrame(SP, data, match.ID=FALSE) 470 } else stop(paste("Incompatible geometry:", u_eType)) 471 472 res 473} 474 475showWKT <- function(p4s, file=NULL, morphToESRI=FALSE, enforce_xy=NULL) { 476 477 if (!is.character(p4s)) stop("invalid p4s object") 478 stopifnot(length(p4s) == 1) 479 if (!is.logical(morphToESRI)) stop("invalid morphToESRI object") 480 stopifnot(length(morphToESRI) == 1) 481 stopifnot(!is.na(morphToESRI)) 482 if (!is.null(enforce_xy)) { 483 stopifnot(is.logical(enforce_xy)) 484 stopifnot(length(enforce_xy) == 1L) 485 stopifnot(!is.na(enforce_xy)) 486 } else { 487 enforce_xy <- get_enforce_xy() 488 } 489 attr(morphToESRI, "enforce_xy") <- enforce_xy 490 res <- .Call("p4s_to_wkt", as.character(p4s), morphToESRI, 491 PACKAGE="rgdal") 492 if (!is.null(file)) cat(res, "\n", sep="", file=file) 493 res 494} 495 496showP4 <- function(wkt, morphFromESRI=FALSE, enforce_xy=NULL) { 497 498 if (!is.character(wkt)) stop("invalid wkt object") 499 stopifnot(length(wkt) == 1) 500 if (!is.logical(morphFromESRI)) stop("invalid morphFromESRI object") 501 stopifnot(length(morphFromESRI) == 1) 502 stopifnot(!is.na(morphFromESRI)) 503 if (!is.null(enforce_xy)) { 504 stopifnot(is.logical(enforce_xy)) 505 stopifnot(length(enforce_xy) == 1L) 506 stopifnot(!is.na(enforce_xy)) 507 } else { 508 enforce_xy <- get_enforce_xy() 509 } 510 attr(morphFromESRI, "enforce_xy") <- enforce_xy 511 res <- .Call("wkt_to_p4s", as.character(wkt), 512 morphFromESRI, PACKAGE="rgdal") 513 res 514} 515 516 517showEPSG <- function(p4s, enforce_xy=NULL) { 518 519 if (!is.character(p4s)) stop("invalid p4s object") 520 stopifnot(length(p4s) == 1L) 521 stopifnot(!is.na(p4s)) 522 if (!is.null(enforce_xy)) { 523 stopifnot(is.logical(enforce_xy)) 524 stopifnot(length(enforce_xy) == 1L) 525 stopifnot(!is.na(enforce_xy)) 526 } else { 527 enforce_xy <- get_enforce_xy() 528 } 529 attr(p4s, "enforce_xy") <- enforce_xy 530 531 res <- .Call("ogrAutoIdentifyEPSG", p4s, PACKAGE="rgdal") 532 res 533} 534 535get_thin_PROJ6_warnings <- function() { 536 get("thin_PROJ6_warnings", envir=.RGDAL_CACHE) 537} 538 539set_thin_PROJ6_warnings <- function(value) { 540 stopifnot(is.logical(value)) 541 stopifnot(length(value) == 1L) 542 stopifnot(!is.na(value)) 543 assign("thin_PROJ6_warnings", value, envir=.RGDAL_CACHE) 544} 545 546get_rgdal_show_exportToProj4_warnings <- function() { 547 get("rgdal_show_exportToProj4_warnings", envir=.RGDAL_CACHE) 548} 549 550set_rgdal_show_exportToProj4_warnings <- function(value) { 551 stopifnot(is.logical(value)) 552 stopifnot(length(value) == 1L) 553 stopifnot(!is.na(value)) 554 assign("rgdal_show_exportToProj4_warnings", value, envir=.RGDAL_CACHE) 555} 556 557get_PROJ6_warnings_count <- function() { 558 get("PROJ6_warnings_count", envir=.RGDAL_CACHE) 559} 560 561get_enforce_xy <- function() { 562 get("enforce_xy", envir=.RGDAL_CACHE) 563} 564 565set_enforce_xy <- function(value) { 566 stopifnot(is.logical(value)) 567 stopifnot(length(value) == 1L) 568 stopifnot(!is.na(value)) 569 assign("enforce_xy", value, envir=.RGDAL_CACHE) 570} 571 572 573get_prefer_proj <- function() { 574 get("prefer_proj", envir=.RGDAL_CACHE) 575} 576 577set_prefer_proj <- function(value) { 578 stopifnot(is.logical(value)) 579 stopifnot(length(value) == 1L) 580 stopifnot(!is.na(value)) 581 assign("prefer_proj", value, envir=.RGDAL_CACHE) 582} 583 584 585showSRID <- function(inSRID, format="WKT2", multiline="NO", enforce_xy=NULL, EPSG_to_init=TRUE, prefer_proj=NULL) { 586 valid_WKT_formats <- c("SFSQL", "WKT1_SIMPLE", "WKT1", "WKT1_GDAL", 587 "WKT1_ESRI", "WKT2_2015", "WKT2_2018", "WKT2") # add WKT2_2019 ?? 588 valid_formats <- c("PROJ", valid_WKT_formats) 589 stopifnot(is.character(inSRID)) 590 stopifnot(length(inSRID) == 1L) 591 stopifnot(nzchar(inSRID)) 592 stopifnot(is.character(format)) 593 stopifnot(length(format) == 1L) 594 if (!(format %in% valid_formats)) stop("invalid format value") 595 out_format <- as.integer(NA) 596 if (format %in% valid_WKT_formats) out_format <- 1L 597 if (format == "PROJ") out_format <- 2L 598 if (!is.na(multiline) && is.logical(multiline)) { 599 multiline <- ifelse(multiline, "YES", "NO") 600 } 601 stopifnot(!is.na(multiline)) 602 stopifnot(is.character(multiline)) 603 stopifnot(length(multiline) == 1L) 604 if (!(multiline %in% c("YES", "NO"))) stop("invalid multiline value") 605 in_format <- as.integer(NA) 606 if (substring(inSRID, 1, 1) == " ") stop("string starts with space") 607 if (substring(inSRID, 1, 1) == "+") in_format = 1L 608 if (substring(inSRID, 1, 3) == "urn") in_format = 2L 609 if (substring(inSRID, 1, 2) == "PR") in_format = 3L 610 if (substring(inSRID, 1, 2) == "GE") in_format = 3L 611 if (substring(inSRID, 1, 2) == "BA") in_format = 3L 612 if (substring(inSRID, 1, 2) == "BO") in_format = 3L 613 if (substring(inSRID, 1, 1) == "S") in_format = 3L 614 if (substring(inSRID, 1, 2) == "CO") in_format = 3L 615 if (substring(inSRID, 1, 3) == "ENG") in_format = 3L 616 if (substring(inSRID, 1, 4) == "EPSG") in_format = 4L 617 if (substring(inSRID, 1, 4) == "ESRI") in_format = 5L 618 if (substring(inSRID, 1, 3) == "OGC") in_format = 5L 619 epsg <- as.integer(NA) 620 if (!is.null(prefer_proj)) { 621 stopifnot(is.logical(prefer_proj)) 622 stopifnot(length(prefer_proj) == 1L) 623 stopifnot(!is.na(prefer_proj)) 624 } else { 625 prefer_proj <- get_prefer_proj() 626 } 627 if (!is.na(in_format) && in_format == 4L) { 628 if (EPSG_to_init && !prefer_proj) { 629 in_format = 1L 630 inSRID <- paste0("+init=epsg:", substring(inSRID, 6, nchar(inSRID))) 631 } else { 632 epsg <- as.integer(substring(inSRID, 6, nchar(inSRID))) 633 } 634 } 635 format <- paste0("FORMAT=", format) 636 multiline <- paste0("MULTILINE=", multiline) 637 if (!is.null(enforce_xy)) { 638 stopifnot(is.logical(enforce_xy)) 639 stopifnot(length(enforce_xy) == 1L) 640 stopifnot(!is.na(enforce_xy)) 641 } else { 642 enforce_xy <- get_enforce_xy() 643 } 644 645 if (new_proj_and_gdal()) { 646 if (!is.na(in_format)) { 647 attr(in_format, "enforce_xy") <- enforce_xy 648 if (prefer_proj) { 649 if (in_format == 1L && !grepl("\\+type\\=crs", inSRID)) 650 inSRID <- paste0(inSRID, " +type=crs") 651 res <- try(.Call("P6_SRID_proj", as.character(inSRID), 652 as.character(format), as.character(multiline), 653 in_format, as.integer(epsg), 654 as.integer(out_format), PACKAGE="rgdal"), silent=TRUE) 655 } else { 656 res <- try(.Call("P6_SRID_show", as.character(inSRID), 657 as.character(format), as.character(multiline), 658 in_format, as.integer(epsg), 659 as.integer(out_format), PACKAGE="rgdal"), silent=TRUE) 660 } 661 if (inherits(res, "try-error")) 662 stop(unclass(attr(res, "condition"))$message, "\n", inSRID) 663 no_towgs84 <- ((is.null(attr(res, "towgs84"))) || 664 (all(nchar(attr(res, "towgs84")) == 0))) 665 if ((length(grep("towgs84|TOWGS84|Position Vector|Geocentric translations", c(res))) == 0L) 666 && !no_towgs84) warning("TOWGS84 discarded") 667 if ((!is.null(attr(res, "ellps"))) && 668 (nchar(attr(res, "ellps")) > 0L) && 669 (length(grep("ellps|ELLIPSOID", c(res))) == 0L)) { 670 if (length(grep("datum|DATUM", c(res))) == 0L) { 671 msg <- paste0("Discarded ellps ", attr(res, "ellps"), 672 " in Proj4 definition: ", c(res)) 673 } else { 674 msg <- "" 675 } 676 if (get_rgdal_show_exportToProj4_warnings()) { 677 if (!get_thin_PROJ6_warnings()) { 678 if (nchar(msg) > 0) warning(msg) 679 } else { 680 if (nchar(msg) > 0 && get("PROJ6_warnings_count", 681 envir=.RGDAL_CACHE) == 0L) { 682 warning(paste0("PROJ/GDAL PROJ string degradation in workflow\n repeated warnings suppressed\n ", msg)) 683 } 684 assign("PROJ6_warnings_count", 685 get("PROJ6_warnings_count", 686 envir=.RGDAL_CACHE) + 1L, envir=.RGDAL_CACHE) 687 } 688 } 689 } 690#warning("Discarded ellps ", attr(res, "ellps"), 691# " in Proj4 definition") 692 if ((!is.null(attr(res, "datum"))) 693 && (nchar(attr(res, "datum")) > 0L) 694# spatsoc tripped PROJ 8.0.0 ENSEMBLE 2021-02-22 695 && (length(grep("datum|DATUM|ENSEMBLE", c(res))) == 0L)) { 696 msg <- paste0("Discarded datum ", attr(res, "datum"), 697 " in Proj4 definition") 698 if (!no_towgs84 && (length(grep("towgs84", c(res))) > 0L)) 699 msg <- paste0(msg, ",\n but +towgs84= values preserved") 700 if (get_P6_datum_hard_fail()) stop(msg) 701 else { 702 if (get_rgdal_show_exportToProj4_warnings()) { 703 if (!get_thin_PROJ6_warnings()) { 704 warning(msg) 705 } else { 706 if (get("PROJ6_warnings_count", 707 envir=.RGDAL_CACHE) == 0L) { 708 warning(paste0("PROJ/GDAL PROJ string degradation in workflow\n repeated warnings suppressed\n ", msg)) 709 } 710 assign("PROJ6_warnings_count", 711 get("PROJ6_warnings_count", 712 envir=.RGDAL_CACHE) + 1L, envir=.RGDAL_CACHE) 713 } 714 } 715 } 716 } 717 res <- c(res) 718 } else { 719 stop("unknown input format") 720 } 721 } else { 722 if (!is.na(in_format) && in_format == 1L) { 723 res <- showWKT(inSRID) 724 } else { 725 stop("unknown input format") 726 } 727 } 728 res 729} 730 731