From 2c3b79f50da58fafe6d9c6d6e25eb43c6f06270f Mon Sep 17 00:00:00 2001 From: Jean-Romain Date: Fri, 3 Sep 2021 09:32:09 -0400 Subject: [PATCH] Remove most use of rgdal #466 --- NEWS.md | 1 + R/io_writeANY.R | 84 +------------ R/projection.R | 115 +++++------------- .../testthat/test-catalog-engine-automerge.R | 6 - tests/testthat/test-classify_noise.R | 2 - tests/testthat/test-grid_metrics.R | 2 +- tests/testthat/test-grid_terrain.R | 3 +- tests/testthat/test-las_check.R | 2 +- tests/testthat/test-projection.R | 70 ++++------- tests/testthat/test-spTransform.R | 16 ++- 10 files changed, 81 insertions(+), 220 deletions(-) diff --git a/NEWS.md b/NEWS.md index 0463e5cba..46b813f11 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,7 @@ If you are viewing this file on CRAN, please check [the latest news on GitHub](h - Fix: `quantize()` now preserves `NaN` values instead of converting them into minus infinite ([#460](https://github.com/Jean-Romain/lidR/issues/460)). - Fix: `stdmetrics_i()` now fails with an informative message when the sum of intensities is greater than `.Machine$integer.max` and become `double` ([#463](https://github.com/Jean-Romain/lidR/issues/463)) +- Change: `rgdal` will be retired in 2024. Code using `rgdal` internally has been removed. In many cases this will not change anything for users but in some cases it may fail at assigning an EPSG code to the LAS file. Also old version of `rgdal` built with old version of `gdal` and `proj` are no longer supported ([#466](https://github.com/Jean-Romain/lidR/issues/466)) ## lidR v3.1.4 (Release date: 2021-06-22) diff --git a/R/io_writeANY.R b/R/io_writeANY.R index 9b074b657..51d916dc6 100644 --- a/R/io_writeANY.R +++ b/R/io_writeANY.R @@ -54,87 +54,15 @@ writeANY = function(x, path, drivers) return(path) } -# Function from raster and adapted to support multiple output format writeSpatial = function(x, filename, overwrite = FALSE, ...) { filename <- normalizePath(filename, winslash = "/", mustWork = FALSE) - driver <- guess_driver_can_write(filename) - layer <- tools::file_path_sans_ext(basename(filename)) + x <- sf::st_as_sf(x) - if (file.exists(filename) & !overwrite) - stop('file exists, use overwrite=TRUE to overwrite it') - - if (!inherits(x, 'Spatial')) - stop('To write a shapefile you need to provide an object of class Spatial*') # nocov - - if (inherits(x, 'SpatialGrid')) - stop('These data cannot be written to a shapefile') # nocov - - if (inherits(x, 'SpatialPixels')) - { - if (.hasSlot(x, 'data')) - x <- as(x, 'SpatialPointsDataFrame') - else - x <- as(x, 'SpatialPoints') - - warning('Writing SpatialPixels to a shapefile. Writing to a raster file format might be more desirable') - } - - if (!.hasSlot(x, 'data')) - { - if (inherits(x, 'SpatialPolygons')) - x <- sp::SpatialPolygonsDataFrame(x, data.frame(ID = 1:length(x)), match.ID = FALSE) - else if (inherits(x, 'SpatialLines')) - x <- sp::SpatialLinesDataFrame( x, data.frame(ID = 1:length(x)), match.ID = FALSE) - else if (inherits(x, 'SpatialPoints')) - x <- sp::SpatialPointsDataFrame( x, data.frame(ID = 1:length(x)), match.ID = FALSE) - else - stop('These data cannot be written to a shapefile') # no cov - } - - rgdal::writeOGR(x, filename, layer, driver = driver, overwrite_layer = overwrite, ...) -} - -# Function non exported from sf and build from sf source code because CRAN does not accept to use operator ::: -# nocov start -guess_driver_can_write = function(dsn) -{ - assert_is_a_string(dsn) - - extension_map <- list( - bna = "BNA", csv = "CSV", e00 = "AVCE00", gdb = "OpenFileGDB", - geojson = "GeoJSON", gml = "GML", gmt = "GMT", gpkg = "GPKG", - gps = "GPSBabel", gtm = "GPSTrackMaker", gxt = "Geoconcept", - jml = "JML", kml = "KML", map = "WAsP", mdb = "Geomedia", - nc = "netCDF", ods = "ODS", osm = "OSM", pbf = "OSM", shp = "ESRI Shapefile", - sqlite = "SQLite", vdv = "VDV", xls = "xls", xlsx = "XLSX") - - prefix_map <- list( - couchdb = "CouchDB", db2odbc = "DB2ODBC", dods = "DODS", - gft = "GFT", mssql = "MSSQLSpatial", mysql = "MySQL", oci = "OCI", - odbc = "ODBC", pg = "PostgreSQL", sde = "SDE") - - drivers <- sf::st_drivers() - - drv <- extension_map[tolower(tools::file_ext(dsn))] - - if (any(grep(":", gsub(":[/\\]", "/", dsn)))) - drv <- prefix_map[tolower(strsplit(dsn, ":")[[1]][1])] - - drv <- unlist(drv) - - if (is.null(drv)) - stop("Could not guess driver for ", dsn, call. = FALSE) - - if (is.na(drv)) - stop("Could not guess driver for ", dsn, call. = FALSE) - - if (is.na(match(drv, drivers$name))) - stop(unlist(drv), " driver not available in supported drivers.'", call. = FALSE) - - if (!drivers[match(drv, drivers$name), "write"]) - stop("Driver ", drv, " cannot write.", call. = FALSE) + if (isTRUE(overwrite)) + append = FALSE + else + append = NA - return(drv) + sf::st_write(x, filename, append = append, quiet = TRUE) } -# nocov end diff --git a/R/projection.R b/R/projection.R index 926a752d2..0368019c9 100644 --- a/R/projection.R +++ b/R/projection.R @@ -241,50 +241,44 @@ setMethod("projection<-", "LAS", function(x, value) # header depending on the LAS format (1.4 or above) # In addition we need this to work both with old and new proj and gdal - proj6 <- rgdal::new_proj_and_gdal() - if (is(value, "CRS")) { proj4 <- value@projargs CRS <- value - wkt <- if (proj6) comment(value) else rgdal::showWKT(proj4) - if (is.null(wkt)) wkt = "" - epsg <- .find_epsg_code(CRS) + crs <- sf::st_crs(CRS) + wkt <- crs$wkt + epsg <- crs$epsg } else if (is(value, "crs")) { - if (proj6) - { - wkt <- value$wkt - CRS <- sp::CRS(SRS_string = wkt) - proj4 <- CRS@projargs - epsg <- .find_epsg_code(value) - } - else # nocov start - { - proj4 <- value$proj4string - CRS <- sp::CRS(proj4) - wkt <- rgdal::showWKT - epsg <- .find_epsg_code(value) - } # nocov end + wkt <- value$wkt + CRS <- sp::CRS(SRS_string = wkt) + crs <- value + proj4 <- CRS@projargs + epsg <- crs$epsg } else if (is.character(value)) { CRS <- sp::CRS(SRS_string = value) proj4 <- CRS@projargs - wkt <- if (proj6) comment(CRS) else rgdal::showWKT(proj4) - epsg <- .find_epsg_code(CRS) + crs <- sf::st_crs(CRS) + wkt <- crs$wkt + epsg <- crs$epsg } else if (is.numeric(value)) { epsg <- value CRS <- epsg2CRS(epsg) + crs <- sf::st_crs(CRS) proj4 <- CRS@projargs - wkt <- if (proj6) comment(CRS) else rgdal::showWKT(proj4) + wkt <- crs$wkt } else stop("'value' is not a CRS or a string or a number.") + if (is.null(epsg) | is.na(epsg)) + epsg = 0 + if (use_wktcs(x)) { if (is.null(wkt) || wkt == "") @@ -457,70 +451,24 @@ epsg2CRS <- function(epsg, fail = FALSE) wkt2CRS <- function(wkt, fail = FALSE) { - if (utils::packageVersion("rgdal") > "1.4.8" && rgdal::new_proj_and_gdal()) + if (!rgdal::new_proj_and_gdal()) { - crs <- tryCatch( - { - sp::CRS(SRS_string = wkt) - }, - error = function(e) - { - if (!fail) - return(sp::CRS(NA_character_)) - else - stop("Invalid WKT string", call. = FALSE) - }) + return(sp::CRS(NA_character_)) } - else # nocov start - { - proj <- wkt2proj(wkt, fail) - crs <- sp::CRS(proj, doCheckCRSArgs = FALSE) # doCheckCRSArgs = FALSE added in 2.2.4 after #323 - } # nocov end - return(crs) -} - -.find_epsg_code <- function(x) -{ - if (is(x, "CRS")) + crs <- tryCatch( { - if (is.na(x@projargs)) - return(0L) - - # Extract epsg code if any - epsg <- sub("\\+init=epsg:(\\d+).*", "\\1", x@projargs) - epsg <- suppressWarnings(as.integer(epsg)) - - if (!is.na(epsg)) return(epsg) - - ## We were not able to extract the epsg code earlier, we can retrieve it - ## with rgdal and proj4 string - epsg <- rgdal::showEPSG(x@projargs) - - # We are still unable to find an epsg code - if (epsg == "OGRERR_UNSUPPORTED_SRS") - return(0L) - else - return(epsg) - } - else if (is(x, "crs")) + sp::CRS(SRS_string = wkt) + }, + error = function(e) { - if (!is.null(x$epsg)) - return(x$epsg) - - if (is.na(x$input)) - return(0L) - - epsg <- sub("\\EPSG:(\\d+).*", "\\1", x$input) - epsg <- suppressWarnings(as.integer(epsg)) - - if (!is.na(epsg)) - return(epsg) + if (!fail) + return(sp::CRS(NA_character_)) else - return(0L) - } - else - stop("Internal error: x is not a CRS or a crs", call. = FALSE) + stop("Invalid WKT string", call. = FALSE) + }) + + return(crs) } # nocov start @@ -528,8 +476,8 @@ epsg2proj <- function(epsg, fail = FALSE) { tryCatch( { - wkt <- rgdal::showWKT(glue::glue("+init=epsg:{epsg}")) - rgdal::showP4(wkt) + crs <- sp::CRS(SRS_string = glue::glue("EPSG:{epsg}")) + crs@projargs }, error = function(e) { @@ -544,7 +492,8 @@ wkt2proj <- function(wkt, fail = FALSE) { tryCatch( { - rgdal::showP4(wkt) + crs <- sp::CRS(SRS_string = glue::glue("EPSG:{epsg}")) + crs@projargs }, error = function(e) { diff --git a/tests/testthat/test-catalog-engine-automerge.R b/tests/testthat/test-catalog-engine-automerge.R index fb89c029b..98c18eaf6 100644 --- a/tests/testthat/test-catalog-engine-automerge.R +++ b/tests/testthat/test-catalog-engine-automerge.R @@ -5,12 +5,6 @@ ctg@output_options$drivers$Raster$param$overwrite = TRUE ctg@output_options$drivers$Spatial$param$overwrite = TRUE ctg@output_options$drivers$SimpleFeature$param$delete_dsn = TRUE -if (utils::packageVersion("rgdal") > "1.4.8" && rgdal::new_proj_and_gdal()) { - expected_crs <- sp::CRS(SRS_string = "EPSG:26917") -} else { - expected_crs <- sp::CRS("+init=epsg:26917") -} - rtest <- function(cluster, layers = 1L) { las <- readLAS(cluster) if (is.empty(las)) return(NULL) diff --git a/tests/testthat/test-classify_noise.R b/tests/testthat/test-classify_noise.R index b9706352a..0340bedae 100644 --- a/tests/testthat/test-classify_noise.R +++ b/tests/testthat/test-classify_noise.R @@ -1,7 +1,5 @@ context("classify_noise") -rgdal::set_thin_PROJ6_warnings(TRUE) - las = clip_rectangle(topography, 273450, 5274350, 273550, 5274450) set.seed(314) diff --git a/tests/testthat/test-grid_metrics.R b/tests/testthat/test-grid_metrics.R index 666e1f398..ca92b89e1 100644 --- a/tests/testthat/test-grid_metrics.R +++ b/tests/testthat/test-grid_metrics.R @@ -1,7 +1,7 @@ context("grid_metrics") las <- lidR:::generate_las(2000) -projection(las) <- sp::CRS("+init=epsg:4326") +projection(las) <- 4326 test_that("grid_metrics returns a named RasterLayer", { diff --git a/tests/testthat/test-grid_terrain.R b/tests/testthat/test-grid_terrain.R index b22650f87..c9f9eff4e 100644 --- a/tests/testthat/test-grid_terrain.R +++ b/tests/testthat/test-grid_terrain.R @@ -1,8 +1,7 @@ context("grid_terrain") -rgdal::set_thin_PROJ6_warnings(TRUE) las <- lidR:::generate_las(5000) -projection(las) <- sp::CRS("+init=epsg:4326") +projection(las) <- 4326 las@data[, Z := round(Z + 0.1*X + 0.1*Y + sin(0.01*X) - sin(0.1*Y) + sin(0.003*X*Y),3)] tdtm <- lidR:::rOverlay(las, 1) diff --git a/tests/testthat/test-las_check.R b/tests/testthat/test-las_check.R index 154fe7583..f006bb0ad 100644 --- a/tests/testthat/test-las_check.R +++ b/tests/testthat/test-las_check.R @@ -24,7 +24,7 @@ las1@header@PHB[["Point Data Format ID"]] <- 25 LASfile <- system.file("extdata", "extra_byte.laz", package = "rlas") las2 <- readLAS(LASfile, select = "xyz") las2@header@PHB$`Global Encoding`$WKT = TRUE -wkt(las2) <- rgdal::showWKT(sp::CRS("+init=epsg:26917")@projargs) +wkt(las2) <- "PROJCS[\"unknown\",GEOGCS[\"unknown\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-81],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]" las2@proj4string <- sp::CRS() las3 <- las1 diff --git a/tests/testthat/test-projection.R b/tests/testthat/test-projection.R index 07ab98ce5..feb8f51fa 100644 --- a/tests/testthat/test-projection.R +++ b/tests/testthat/test-projection.R @@ -4,14 +4,11 @@ las = random_10_points test_that("Internal projection conversion works", { - wkt <- rgdal::showWKT("+init=epsg:2008") + skip_if(!rgdal::new_proj_and_gdal()) - if (rgdal::new_proj_and_gdal()) { - expected <- sp::CRS(SRS_string = "EPSG:2008") - } else { - expected <- sp::CRS("+init=epsg:2008") - expected@projargs <- sub("\\+init=epsg:\\d+\\s", "", expected@projargs) - } + wkt <- "PROJCS[\"NAD27(CGQ77) / SCoPQ zone 2\",GEOGCS[\"NAD27(CGQ77)\",DATUM[\"North_American_Datum_1927_CGQ77\",SPHEROID[\"Clarke 1866\",6378206.4,294.978698213898,AUTHORITY[\"EPSG\",\"7008\"]],AUTHORITY[\"EPSG\",\"6609\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4609\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-55.5],PARAMETER[\"scale_factor\",0.9999],PARAMETER[\"false_easting\",304800],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]" + + expected <- sp::CRS(SRS_string = "EPSG:2008") expect_equal(lidR:::epsg2CRS(2008), expected) expect_equal(lidR:::epsg2CRS(200800), sp::CRS()) @@ -36,68 +33,51 @@ test_that("projection<- with EPSG works", { test_that("projection<- with CRS works", { - projection(las) <- sp::CRS("+init=epsg:26917") + skip_if(!rgdal::new_proj_and_gdal()) + projection(las) <- sp::CRS(SRS_string = "EPSG:26917") expect_equal(epsg(las), 26917) - projection(las) <- sp::CRS("+init=epsg:26918") - + projection(las) <- sp::CRS(SRS_string = "EPSG:26918") expect_equal(epsg(las), 26918) }) + test_that("crs<- with epsg code works", { - crs(las) <- sp::CRS("+init=epsg:26917") + skip_if(!rgdal::new_proj_and_gdal()) + crs(las) <- sp::CRS(SRS_string = "EPSG:26917") expect_equal(epsg(las), 26917) - crs(las) <- sp::CRS("+init=epsg:26918") - + crs(las) <- sp::CRS(SRS_string = "EPSG:26918") expect_equal(epsg(las), 26918) }) test_that("projection<- with wkt code works", { + skip_if(!rgdal::new_proj_and_gdal()) + las@header@PHB[["Global Encoding"]][["WKT"]] <- TRUE - if (rgdal::new_proj_and_gdal()) - { - projection(las) <- sp::CRS(SRS_string = "EPSG:26919") - expect_match(wkt(las), "NAD83 / UTM zone 19N") - - projection(las) <- sp::CRS(SRS_string = "EPSG:26918") - expect_match(wkt(las), "NAD83 / UTM zone 18N") - } - else - { - projection(las) <- sp::CRS("+init=epsg:26919") - expect_match(wkt(las), "PROJCS") - - projection(las) <- sp::CRS("+init=epsg:26918") - expect_match(wkt(las), "PROJCS") - } + projection(las) <- sp::CRS(SRS_string = "EPSG:26919") + expect_match(wkt(las), "NAD83 / UTM zone 19N") + + projection(las) <- sp::CRS(SRS_string = "EPSG:26918") + expect_match(wkt(las), "NAD83 / UTM zone 18N") }) test_that("crs<- with wkt code works", { + skip_if(!rgdal::new_proj_and_gdal()) + las@header@PHB[["Global Encoding"]][["WKT"]] <- TRUE - if (rgdal::new_proj_and_gdal()) - { - crs(las) <- sp::CRS(SRS_string = "EPSG:26919") - expect_match(wkt(las), "NAD83 / UTM zone 19N") - - crs(las) <- sp::CRS(SRS_string = "EPSG:26918") - expect_match(wkt(las), "NAD83 / UTM zone 18N") - } - else - { - crs(las) <- sp::CRS("+init=epsg:26919") - expect_match(wkt(las), "PROJCS") - - crs(las) <- sp::CRS("+init=epsg:26918") - expect_match(wkt(las), "PROJCS") - } + crs(las) <- sp::CRS(SRS_string = "EPSG:26919") + expect_match(wkt(las), "NAD83 / UTM zone 19N") + + crs(las) <- sp::CRS(SRS_string = "EPSG:26918") + expect_match(wkt(las), "NAD83 / UTM zone 18N") }) test_that("epsg<- works", { diff --git a/tests/testthat/test-spTransform.R b/tests/testthat/test-spTransform.R index 9411fc97e..9b6d8ccdc 100644 --- a/tests/testthat/test-spTransform.R +++ b/tests/testthat/test-spTransform.R @@ -2,11 +2,17 @@ context("spTransform") LASfile <- system.file("extdata", "Megaplot.laz", package = "lidR") las <- readLAS(LASfile, select = "xyz") -crs <- sp::CRS("+init=epsg:26918") # Fails on Solaris: crs <- sp::CRS(SRS_string = "EPSG:26918") -las2 <- spTransform(las, crs) +if (rgdal::new_proj_and_gdal()) +{ + crs <- sp::CRS(SRS_string = "EPSG:26918") + las2 <- spTransform(las, crs) +} test_that("Datum tranformation works", { + + skip_if(!rgdal::new_proj_and_gdal()) + expect_equal(mean(las$X), 684879, tol = 1) expect_equal(mean(las$Y), 5017900, tol = 1) expect_equal(mean(las2$X), 214383, tol = 1) @@ -14,11 +20,17 @@ test_that("Datum tranformation works", { }) test_that("Bounding box is updated", { + + skip_if(!rgdal::new_proj_and_gdal()) + expect_equivalent(bbox(las), matrix(c(684766.39, 5017773.08, 684993.29, 5018007.25), ncol = 2)) expect_equivalent(bbox(las2), matrix(c(214261.66, 5021517.24, 214504.91, 5021767.46), ncol = 2), tol = 0.01) }) test_that("Bounding header is updated", { + + skip_if(!rgdal::new_proj_and_gdal()) + expect_equivalent(unlist(las@header@PHB[c("Min X", "Min Y", "Max X", "Max Y")]), c(684766.39, 5017773.08, 684993.29, 5018007.25)) expect_equivalent(unlist(las2@header@PHB[c("Min X", "Min Y", "Max X", "Max Y")]), c(214261.66, 5021517.24, 214504.91, 5021767.46), tol = 0.1) })