1# Copyright 2003 (c) Barry Rowlingson 2# Modified 2006-12 Roger Bivand 3### 4### 5### Routines for ogr layer data source 6### 7### 8 9# 10ogrInfo <- function(dsn, layer, encoding=NULL, 11 use_iconv=FALSE, swapAxisOrder=FALSE, require_geomType=NULL, 12 morphFromESRI=NULL, dumpSRS=FALSE, enforce_xy=NULL, D3_if_2D3D_points=FALSE) { 13 if (missing(dsn)) stop("missing dsn") 14 stopifnot(is.character(dsn)) 15 stopifnot(length(dsn) == 1L) 16# copy sf::st_read.default usage 17 if (length(dsn) == 1 && file.exists(dsn)) 18 dsn <- enc2utf8(normalizePath(dsn)) 19 if (nchar(dsn) == 0) stop("empty dsn") 20 if (missing(layer)) { 21 layers <- ogrListLayers(dsn=dsn) 22 if (length(layers) == 0L) stop("missing layer") 23 if (length(layers) > 0L) layer <- enc2utf8(c(layers[1])) 24 if (length(layers) > 1L) 25 warning("First layer ", layer, 26 " read; multiple layers present in\n", dsn, 27 ", check layers with ogrListLayers()") 28 29 } else layer <- enc2utf8(layer) 30 if (nchar(layer) == 0) stop("empty name") 31 WKB <- c("wkbPoint", "wkbLineString", "wkbPolygon", "wkbMultiPoint", 32 "wkbMultiLineString", "wkbMultiPolygon", "wkbGeometryCollection") 33 if (!is.null(require_geomType)) { 34 stopifnot(is.character(require_geomType) && length(require_geomType)==1) 35 m_require_geomType <- match(require_geomType, WKB) 36 stopifnot(!is.na(m_require_geomType) || m_require_geomType <= 3) 37 } 38# a list with various ogr data source information 39 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 if (!use_iconv && !is.null(encoding)) { 47 oSE <- getCPLConfigOption("SHAPE_ENCODING") 48 tull <- setCPLConfigOption("SHAPE_ENCODING", encoding) 49 } 50 ogrinfo <- .Call("RGDAL_ogrInfo",as.character(dsn), as.character(layer), 51 PACKAGE = "rgdal") 52 if (!use_iconv && !is.null(encoding)) { 53 tull <- setCPLConfigOption("SHAPE_ENCODING", oSE) 54 } 55 56 if (swapAxisOrder) ogrinfo[[5]] <- ogrinfo[[5]][c(2,1,4,3)] 57 58 u_eType <- u_with_z <- null_geometries <- NULL 59 deleted_geometries <- NULL 60 retain <- NULL 61 have_features <- NULL 62 all_NULL <- FALSE 63 keepGeoms <- NULL 64 65 if (!is.na(ogrinfo[[1]])) { 66 67 fids <- ogrFIDs(dsn=dsn, layer=layer) 68 nrows_i <- attr(fids, "i") 69 have_features <- nrows_i > 0 70 if (have_features && (attr(fids, "i") != attr(fids, "nf"))) { 71 retain <- 1:attr(fids, "i") 72 afids <- 0:(attr(fids, "nf")-1) 73 deleted <- afids[!(afids %in% fids[retain])] 74 deleted_geometries <- paste("Deleted feature IDs:", paste(deleted, 75 collapse=", ")) 76 fids <- fids[retain] 77 } else { 78 deleted_geometries <- NULL 79 retain <- NULL 80 } 81# attributes(fids) <- NULL 82 if (have_features) { 83 eTypes <- .Call("R_OGR_types",as.character(dsn), as.character(layer), 84 PACKAGE = "rgdal") 85 if (is.null(retain)) { 86 eType <- eTypes[[4]] 87 with_z <- eTypes[[5]] 88 isNULL <- as.logical(eTypes[[6]]) 89 } else { 90 eType <- eTypes[[4]][retain] 91 with_z <- eTypes[[5]][retain] 92 isNULL <- as.logical(eTypes[[6]])[retain] 93 } 94 null_geometries <- NULL 95 if (any(isNULL)) { 96 all_NULL <- (sum(isNULL) == length(eType)) 97 eType <- eType[!isNULL] 98 with_z <- with_z[!isNULL] 99 null_geometries <- paste("Null geometry IDs:", 100 paste(which(isNULL), collapse=", ")) 101 } 102 if (!all_NULL) { 103 eType[eType == 5L] <- 2L 104 eType[eType == 6L] <- 3L 105 106 u_eType <- unique(sort(eType)) 107 u_with_z <- unique(sort(with_z)) 108 attr(u_with_z, "elevated") <- FALSE 109 if (length(u_with_z) != 1L) { 110 if (D3_if_2D3D_points) { 111 warning("Accepting mixed 2D and 3D as 3D") 112 u_with_z <- 1L 113 attr(u_with_z, "elevated") <- TRUE 114 } else { 115 stop(paste("Multiple # dimensions:", paste((u_with_z + 2), 116 collapse=":"))) 117 } 118 } 119 if (u_with_z < 0 || u_with_z > 1) stop( 120 paste("Invalid # dimensions:", (u_with_z + 2))) 121 122 t_eType <- table(eType) 123 if (is.null(require_geomType)) { 124 keepGeoms <- NULL 125 if (length(u_eType) > 1L) stop( 126 paste("Multiple incompatible geometries:", 127 paste(paste(WKB[as.integer(names(t_eType))], t_eType, sep=": "), 128 collapse="; "))) 129# if (length(u_eType) == 2L) { 130# if (u_eType[1] == 2 && u_eType[2] == 5) u_eType = 2 131# else if (u_eType[1] == 3 && u_eType[2] == 6) u_eType = 3 132# else stop(paste("Multiple incompatible geometries:", 133# paste(paste(WKB[as.integer(names(t_eType))], t_eType, sep=": "), 134# collapse="; "))) 135# } 136 } else { 137 if (!require_geomType %in% WKB[as.integer(names(t_eType))]) 138 stop(require_geomType, "not in", WKB[as.integer(names(t_eType))]) 139 u_eType <- match(require_geomType, WKB) 140 keepGeoms <- WKB[eType] == require_geomType 141 message("NOTE: keeping only ", sum(keepGeoms), " ", require_geomType, 142 " of ", length(keepGeoms), " features\n", 143 " note that extent applies to all features") 144 } 145 } else { 146 have_features <- FALSE 147 warning("ogrInfo: all features NULL") 148 } 149 } 150 } else { 151 ogrinfo[[1]] <- attr(ogrinfo[[1]], "dFIDs") 152 warning("ogrInfo: feature count overflow") 153 } 154 names(ogrinfo) <- c("nrows", "nitems", "iteminfo", "driver", "extent", 155 "nListFields") 156 if (ogrinfo$driver == "ESRI Shapefile") { 157 DSN <- dsn 158 if (!file.info(DSN)$isdir) DSN <- dirname(enc2utf8(normalizePath(dsn))) 159 DBF_fn <- paste(DSN, .Platform$file.sep, layer, ".dbf", sep = "") 160 if (file.exists(DBF_fn)) { 161 con <- file(DBF_fn, "rb") 162 vr <- readBin(con, "raw", n=32L) 163 ldid <- as.integer(vr[30]) 164 attr(ogrinfo, "LDID") <- ldid 165 close(con) 166 } else { 167 warning("ogrInfo: ", DBF_fn, " not found", sep="") 168 } 169 } 170 names(ogrinfo$iteminfo) <- c("name","type","length","typeName","maxListCount") 171 if (use_iconv && !is.null(encoding)) 172 ogrinfo$iteminfo$name <- iconv(ogrinfo$iteminfo$name, from=encoding) 173 ogrinfo$have_features <- have_features 174 ogrinfo$eType <- u_eType 175 ogrinfo$with_z <- u_with_z 176 ogrinfo$null_geometries <- null_geometries 177 ogrinfo$deleted_geometries <- deleted_geometries 178 ogrinfo$dsn <- dsn 179 ogrinfo$layer <- layer 180 if (is.null(morphFromESRI)) { 181 if (ogrinfo$driver == "ESRI Shapefile") morphFromESRI <- TRUE 182 else morphFromESRI <- FALSE 183 } 184 if (ogrinfo$driver != "ESRI Shapefile" && morphFromESRI) 185 morphFromESRI <- FALSE 186 p4s0 <- OGRSpatialRef(dsn, layer, morphFromESRI=morphFromESRI, 187 dumpSRS=dumpSRS, enforce_xy=enforce_xy) 188 if (new_proj_and_gdal()) wkt2 <- comment(p4s0) 189 ogrinfo$p4s <- c(p4s0) 190 if (new_proj_and_gdal()) ogrinfo$wkt2 <- wkt2 191 if (!is.null(require_geomType)) 192 attr(ogrinfo, "require_geomType") <- require_geomType 193 if (!is.null(keepGeoms)) attr(ogrinfo, "keepGeoms") <- keepGeoms 194 class(ogrinfo) <- "ogrinfo" 195 ogrinfo 196} 197 198print.ogrinfo <- function(x, ...) { 199 cat("Source: \"", x$dsn, '\", layer: \"', x$layer, "\"", '\n', sep='') 200 cat("Driver:", x$driver) 201 if (x$have_features) cat("; number of rows:", x$nrows, "\n") 202 WKB <- c("wkbPoint", "wkbLineString", "wkbPolygon", "wkbMultiPoint", 203 "wkbMultiLineString", "wkbMultiPolygon", "wkbGeometryCollection") 204 if (!is.null(attr(x, "require_geomType"))) { 205 cat(" selected geometry type:", attr(x, "require_geomType"), "with", 206 sum(attr(x, "keepGeoms")), "rows\n") 207 } 208 if (!x$have_features) cat(", no features found\n") 209 if (x$have_features) cat("Feature type:", paste(WKB[x$eType], 210 collapse=", "), "with", x$with_z+2, "dimensions\n") 211 if (!is.null(x$extent)) cat("Extent: (", x$extent[1], " ", 212 x$extent[2], ") - (", x$extent[3], " ", x$extent[4], ")\n", sep="") 213 if (!is.null(x$null_geometries)) cat(x$null_geometries, "\n") 214 if (!is.null(x$deleted_geometries)) cat(x$deleted_geometries, "\n") 215 if ((nchar(x$p4s) > 1) && !is.na(x$p4s)) cat("CRS:", x$p4s, "\n") 216 if (!is.null(attr(x, "LDID"))) cat("LDID:", attr(x, "LDID"), "\n") 217 cat("Number of fields:", x$nitems, "\n") 218 if (is.null(x$nListFields)) x$nListFields <- 0 219 if (x$nListFields > 0) cat("Number of list fields:", x$nListFields, "\n") 220 if (x$nitems > 0 && x$nListFields == 0) print(as.data.frame(x$iteminfo)[,1:4]) 221 if (x$nitems > 0 && x$nListFields > 0 && x$have_features) 222 print(as.data.frame(x$iteminfo)) 223 if (x$nitems > 0 && x$nListFields > 0 && !x$have_features) 224 print(as.data.frame(x$iteminfo)[,1:4]) 225 invisible(x) 226} 227 228 229ogrFIDs <- function(dsn, layer){ 230 if (missing(dsn)) stop("missing dsn") 231 stopifnot(is.character(dsn)) 232 stopifnot(length(dsn) == 1L) 233# copy sf::st_read.default usage 234 if (length(dsn) == 1 && file.exists(dsn)) 235 dsn <- enc2utf8(normalizePath(dsn)) 236 if (nchar(dsn) == 0) stop("empty name") 237 if (missing(layer)) stop("missing layer") 238 layer <- enc2utf8(layer) 239 if (nchar(layer) == 0) stop("empty name") 240 fids <- .Call("RGDAL_ogrFIDs",as.character(dsn),as.character(layer), PACKAGE = "rgdal") 241 if (attr(fids, "i") == 0L) warning("no features found") 242 fids 243} 244 245ogrDrivers <- function() { 246 if (strsplit(getGDALVersionInfo(), " ")[[1]][2] < "2") { 247 res <- .Call("ogr_GetDriverNames", PACKAGE="rgdal") 248 res <- as.data.frame(res, stringsAsFactors=FALSE) 249 } else { 250 res <- .Call('RGDAL_GetDriverNames', PACKAGE="rgdal") 251 if (!is.null(attr(res, "isVector"))) res$isVector <- attr(res, "isVector") 252 res <- as.data.frame(res, stringsAsFactors=FALSE) 253 res <- res[res$isVector,] 254 names(res)[3] <- "write" 255 } 256 res <- res[order(res$name),] 257 row.names(res) <- NULL 258 res 259} 260 261OGRSpatialRef <- function(dsn, layer, morphFromESRI=NULL, dumpSRS=FALSE, 262 driver=NULL, enforce_xy=NULL) { 263 stopifnot(is.character(dsn)) 264 stopifnot(length(dsn) == 1L) 265 if (is.null(driver)) driver <- attr(ogrListLayers(dsn), "driver") 266 if (is.null(morphFromESRI)) { 267 if (driver == "ESRI Shapefile") morphFromESRI <- TRUE 268 else morphFromESRI <- FALSE 269 } 270 if (driver != "ESRI Shapefile" && morphFromESRI) 271 morphFromESRI <- FALSE 272 stopifnot(is.logical(morphFromESRI)) 273 stopifnot(length(morphFromESRI) == 1L) 274# copy sf::st_read.default usage 275 if (length(dsn) == 1 && file.exists(dsn)) 276 dsn <- enc2utf8(normalizePath(dsn)) 277 layer <- enc2utf8(layer) 278 stopifnot(is.logical(dumpSRS)) 279 stopifnot(length(dumpSRS) == 1) 280 stopifnot(!is.na(dumpSRS)) 281 if (!is.null(enforce_xy)) { 282 stopifnot(is.logical(enforce_xy)) 283 stopifnot(length(enforce_xy) == 1L) 284 stopifnot(!is.na(enforce_xy)) 285 } else { 286 enforce_xy <- get_enforce_xy() 287 } 288 attr(dumpSRS, "enforce_xy") <- enforce_xy 289 wkt2 <- NULL 290 no_ellps <- FALSE 291 res <- .Call("ogrP4S", as.character(dsn), as.character(layer), 292 as.logical(morphFromESRI), dumpSRS, PACKAGE="rgdal") 293 if (!is.na(res) && new_proj_and_gdal()) { 294 no_towgs84 <- ((is.null(attr(res, "towgs84"))) || 295 (all(nchar(attr(res, "towgs84")) == 0))) 296 if ((length(grep("towgs84", c(res))) == 0L) && !no_towgs84) 297 warning("TOWGS84 discarded") 298 no_ellps <- (!is.null(attr(res, "ellps"))) && 299 (nchar(attr(res, "ellps")) > 0L) && 300 (length(grep("ellps", c(res))) == 0L) 301 no_ellps <- no_ellps && length(grep("datum", c(res))) == 0L 302 if (no_ellps) { 303 msg <- paste0("Discarded ellps ", attr(res, "ellps"), 304 " in Proj4 definition: ", c(res)) 305 if (get_rgdal_show_exportToProj4_warnings()) { 306 if (!get_thin_PROJ6_warnings()) { 307 warning(msg) 308 } else { 309 if (get("PROJ6_warnings_count", 310 envir=.RGDAL_CACHE) == 0L) { 311 warning(paste0("PROJ/GDAL PROJ string degradation in workflow\n repeated warnings suppressed\n ", msg)) 312 } 313 assign("PROJ6_warnings_count", 314 get("PROJ6_warnings_count", 315 envir=.RGDAL_CACHE) + 1L, envir=.RGDAL_CACHE) 316 } 317 } 318 } 319 if ((!is.null(attr(res, "datum"))) && (nchar(attr(res, "datum")) > 0L) 320 && (length(grep("datum", c(res))) == 0L)) { 321 msg <- paste0("Discarded datum ", attr(res, "datum"), 322 " in Proj4 definition: ", c(res)) 323 if (!no_towgs84 && (length(grep("towgs84", c(res))) > 0L)) 324 msg <- paste0(msg, ",\n but +towgs84= values preserved") 325 if (get_P6_datum_hard_fail()) stop(msg) 326 else { 327 if (get_rgdal_show_exportToProj4_warnings()) { 328 if (!get_thin_PROJ6_warnings()) { 329 warning(msg) 330 } else { 331 if (get("PROJ6_warnings_count", 332 envir=.RGDAL_CACHE) == 0L) { 333 warning(paste0("PROJ/GDAL PROJ string degradation in workflow\n repeated warnings suppressed\n ", msg)) 334 } 335 assign("PROJ6_warnings_count", 336 get("PROJ6_warnings_count", 337 envir=.RGDAL_CACHE) + 1L, envir=.RGDAL_CACHE) 338 } 339 } 340 } 341#warning(msg) 342 } 343 if (new_proj_and_gdal()) wkt2 <- attr(res, "WKT2_2018") 344 } 345 res <- c(res) 346 if (new_proj_and_gdal()) { 347 if (no_ellps) res <- showSRID(wkt2, "PROJ") 348 comment(res) <- wkt2 349 } 350 res 351} 352 353ogrListLayers <- function(dsn) { 354 if (missing(dsn)) stop("missing dsn") 355 stopifnot(is.character(dsn)) 356 stopifnot(length(dsn) == 1) 357# copy sf::st_read.default usage 358 if (length(dsn) == 1 && file.exists(dsn)) 359 dsn <- enc2utf8(normalizePath(dsn)) 360 if (nchar(dsn) == 0) stop("empty name") 361 if (!is.null(attr(dsn, "debug"))) { 362 stopifnot(is.logical(attr(dsn, "debug"))) 363 stopifnot(length(attr(dsn, "debug")) == 1) 364 } else { 365 attr(dsn, "debug") <- FALSE 366 } 367 layers <- .Call("RGDAL_ogrListLayers", dsn, PACKAGE = "rgdal") 368 n <- length(layers) 369 tmp <- layers[n] 370 layers <- layers[-n] 371 Encoding(layers) <- "UTF-8" 372 layers <- enc2native(layers) 373 attr(layers, "driver") <- tmp 374 attr(layers, "nlayers") <- (n-1) 375 layers 376} 377