Remove most use of rgdal #466

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)
 

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

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,79 +451,33 @@ 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
 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)
   {

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)

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)

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", {
 

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)

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