GDAL 3.3.2 changes/deprecates AutoIdentifyEPSG()

[rsbivand]

You use rgdal::showEPSG() in .find_epsg_code() in the "CRS" branch. This fails check with PROJ 8.1.1 and GDAL 3.3.2 (see also r-spatial/sf#1776).
testthat.Rout.fail.zip
You should not use the rgdal function; I do not know when or if these very old interfaces may be updated. Note that rgdal, rgeos and maptools will not be developed further, and will be withdrawn by 2024-01-01 at the latest.

[Jean-Romain]

Thank you for reporting. Let me explain some context so you will be able to guide me toward better solutions.

Objects from lidR inherit Spatial objects from sp. Thus they contain a CRS object from sp. I understand that the new standard is sf but lidR was created before sf. I now have to deal with backward compatibility and it is impossible to fully move to sf without breaking everything. I'm using more an more sf functions internally but the output and classes in lidR are still sp and raster based.

You use rgdal::showEPSG() in .find_epsg_code() in the "CRS" branch.

A LAS file can only store the EPSG code to record the CRS. At R level when a user want to assign a CRS to a LAS object they will likely use another spatial object this way projection(las) <- projection(polygon). If polygon contains a proj4string I must retrieve the EPSG code to be able to write this information to the LAS file. Similar problem arises with crs object from sf if polygon is an sf. Which non-rgdal function should I use to do that?

You should not use the rgdal function; I do not know when or if these very old interfaces may be updated. Note that rgdal, rgeos and maptools will not be developed further, and will be withdrawn by 2024-01-01 at the latest.

I the next release I already removed dependency to rgeos using sf only. I'm not using maptools, at least not at direct dependency. But rgdal yes. I enumerated 12 occurrences of 4 rgdalfunctions that I'm still using despite migrating to sf internally:

• rgdal::writeORG() to write sp object
• rgdal::new_proj_and_gdal() that one gave me a lot of trouble to make my code working on all platform on CRAN and for all users
• rgdal::showWKT() to convert a proj4 string into a WKT string because actually some LAS file format use WKT. But this is wrapped in a if (rgdal::new_proj_and_gdal()) and I'm using comment to retrieve the WKT of CRS object if the rgdal version is recent
• rgdal::showEPSG() as explained above

[rsbivand]

Good to hear of your progress with sf!

rgdal::writeOGR() should be replaced with coercion to "sf" and st_write(). rgdal::showWKT() should use coercion to "crs" and the $wkt method, rgdal::showEPSG() ditto and the $epsg method (which may fail for objects instantiated by +init=epsg: with current PROJ/GDAL). GDAL/PROJ versions can be found from sf_extSoftVersion().

> st_crs(sp::CRS("+init=epsg:26918"))$epsg
[1] NA
> st_crs(sp::CRS("EPSG:26918"))$epsg
[1] 26918
> st_crs("EPSG:26918")$epsg
[1] 26918
> st_crs("EPSG:26918")$wkt
[1] "PROJCRS[\"NAD83 / UTM zone 18N\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"UTM zone 18N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-75,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"North America - between 78°W and 72°W - onshore and offshore. Canada - Nunavut; Ontario; Quebec. United States (USA) - Connecticut; Delaware; Maryland; Massachusetts; New Hampshire; New Jersey; New York; North Carolina; Pennsylvania; Virginia; Vermont.\"],\n        BBOX[28.28,-78,84,-72]],\n    ID[\"EPSG\",26918]]"

[Jean-Romain]

Sure your options are valid but what about this case:

r = sp::CRS("+proj=utm +zone=18 +datum=NAD83 +units=m +no_defs ")
crs = sf::st_crs(r)
crs$epsg
#> NA
rgdal::showEPSG(r@projargs)
#> "26918"

I also made few more tests to see if it works in most probable use cases. It seems to me that showEPSG() works in more cases (see reprex at the end). Anyway, being able to set a CRS to a LAS object from other spatial object is not the main feature and I can safely move to sf even if some cases are no longer handled. However I'm more concerned by long time support and compatibility and I'd like your insight on this.

I said I removed rgeos but actually it is still in Suggest. If I remove rgeos the package check fails because raster no longer works without rgeos. If rgeos is withdrawn then raster too. You said that rgdal will be withdrawn too. This will impact sp and raster as well. I know that we now have sf, stars and terra but as I said I made lidR before sf and terra. Or to be exact I made lidR when sf was very young and I asked if I should base lidR on sf or sp and Edzer Pebesma himself advised me to keep going with sp (in an issue that no longer exist).

Considering the numerous previous issues you opened, I feel that one day or another lidR will no longer work along with raster and sp. I don't know when but this will happen and this day I will be forced to redo everything from strach with sf and terra/stars and break lidR.

So at this stage do you think I should start the work as soon as possible? I feel more and more that maintaining backward compatibility is a dead end issue and one day or another sp or raster (or both) with be removed from CRAN. So I must anticipate these days. And if I must move to more modern packages should I use terra or stars classes to produce raster output?

I'm comfortable with sf that I now use daily instead of sp when doing stuff outside of lidR but I'm still a user of raster and do not well know the differences between stars and terra.

---

# ===============
  
f = system.file("external/lux.shp", package="raster")

p = raster::shapefile(f)
crs = sf::st_crs(p)
crs$epsg
#> [1] 4326

p = sf::st_read(f)
crs = sf::st_crs(p)
crs$epsg
#> [1] 4326

lidR:::.find_epsg_code(crs)
#> [1] 4326

# ===============

f <- system.file("extdata", "lake_polygons_UTM17.shp", package = "lidR")

p = raster::shapefile(f)
crs = sf::st_crs(p)
crs$epsg
#> [1] 26917

p = sf::st_read(f)
crs = sf::st_crs(p)
crs$epsg
#> [1] 26917

lidR:::.find_epsg_code(crs)
#> [1] 26917

# ===========

tif = system.file("tif/L7_ETMs.tif", package = "stars")
r = stars::read_stars(tif)
crs = sf::st_crs(r)
crs$epsg
#> [1] NA

r = terra::rast(tif)
crs = sf::st_crs(r)
crs$epsg
#> [1] NA

r = raster::raster(tif)
#> Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum unknown in Proj4 definition
crs = sf::st_crs(r)
crs$epsg
#> [1] NA

rgdal::showEPSG(crs$input)
#> Error in rgdal::showEPSG(crs$input): Can't parse PROJ.4-style parameter string

# ===============

tif = system.file("ex/meuse.tif", package = "terra")
r = stars::read_stars(tif)
crs = sf::st_crs(r)
crs$epsg
#> [1] NA

r = terra::rast(tif)
crs = sf::st_crs(r)
crs$epsg
#> [1] NA

r = raster::raster(tif)
crs = sf::st_crs(r)
crs$epsg
#> [1] NA

rgdal::showEPSG(crs$input)
#> [1] "OGRERR_UNSUPPORTED_SRS"

# ==============

load(url("http://github.com/mgimond/Spatial/raw/master/Data/raster.RData"))
crs = sf::st_crs(bath)
crs$epsg
#> [1] NA
rgdal::showEPSG(crs$input)
#> [1] "4326"

[Jean-Romain]

crs = sf::st_crs(26917)
crs$epsg
# NULL

[rsbivand]

I've linked your findings to the relevant sf issue, as you see.

Our plan is to reduce risk by retiring rgdal, rgeos and maptools which I maintain, but which might not receive the attention they require. I retired in April, am 70, and now in good health, but the risk is there. For vector work, sf covers effectively everything (maybe stepping carefully around st_transform() for points, which may rather be sf_project() with many points to get to the speed of rgdal::spTransform() for points). For raster work, it is more complex, so sp classes will remain, but beyond that will depend on terra or sf for reading and writing raster files. Edzer announced this policy in his useR keynote.

[edzer]

I think the signals from upstream (GDAL community) are clear: we should stop assuming proj.4 strings can give us EPSG IDs for cases where the proj4 string implicitly assumes axis order long/lat and the EPSG definition has order lat/long.

[rsbivand]

@Jean-Romain , for me:

> library(sf)
Linking to GEOS 3.10.0dev, GDAL 3.3.2, PROJ 8.1.1
> crs = sf::st_crs(26917)
> crs$epsg
[1] 26917

(all most recent). What are your package and extSoftVersion details for #466 (comment) ?

[Jean-Romain]

I watched Edzer's useR keynote on Youtube. I now better understand the state of rgdal and rgeos. So much stuff has been built on top of these packages and we own you a lot. Thank you for everything and happy retirement.

library(sf)
Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

My packages are updated from the ubuntugis-unstable PPA so I guess I need to wait few more weeks/months before to reach your recent versions.

It is now clear that I have to fully move to sf and terra/stars even if it breaks backward compatibility. I can easily create new functions with new names that return sf/terra/stars objects instead of sp/raster (all my functions already accept sf inputs). However my classes are sp based and inherit Spatial.

setClass(
  Class = "LAS", contains = "Spatial",
  representation(data = "data.table", header = "LASheader", index = "list")
)

The goal of such inheritance was to get a slot for CRS and a compatibility by inheritance with other packages. For example

sf::st_crs(las) # works
raster::projection(las) # works

The biggest problem so far is that the CRS is a sp::CRS and I cannot move to a sf::crs. If I rebuild lidR based on sf I'd like to inherit an sf class to get a crs and a compatibility by inheritance. For obvious reasons of memory and speed a LiDAR point cloud cannot be stored into POINT. Here an example with only 80000 points

library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las = readLAS(LASfile)
sflas = sf::st_as_sf(as.spatial(las))
format(object.size(las), "Mb")
#> "5.3 Mb"
format(object.size(sflas), "Mb")
#> "37.7 Mb"

So my question is: which low level sf class, if any, can I inherit to be compatible with sf while not being an sf object.

[Jean-Romain]

That one fails too (at least with my versions of gdal and proj)

crs = sf::st_crs(sp::CRS("+init=epsg:26917"))
crs$epsg
# NA

crs = sf::st_crs(sp::CRS(SRS_string = "EPSG:26917"))
crs$epsg
# 26917

[rsbivand]

> crs = sf::st_crs(sp::CRS("+init=epsg:26917"))
> crs$epsg
[1] NA
> crs = sf::st_crs(26917)
> crs$epsg
[1] 26917
> crs = sf::st_crs("EPSG:26917")
> crs$epsg
[1] 26917
> crs = sf::st_crs("+init=epsg:26917")
> crs$epsg
[1] NA

When user input is a PROJ4 string, the connection is lost, in the WKT2 representation, the ID["EPSG",] node is not present.

[Jean-Romain]

I removed almost all occurrences of rgdal:: including showP4, showWKT, showEPSG, writeOGR. Honestly I'm not sure of what I'm doing. A lot of stuff was actually workarounds for cases where rgdal::new_proj_and_gdal() = FALSE. This include not up-to-date users but also the CRAN-Solaris. For example in my unit tests I have this to pass on Solaris.

crs <- sp::CRS("+init=epsg:26918") # Fails on Solaris: crs <- sp::CRS(SRS_string = "EPSG:26918")

But since I can no longer retrieve the EPSG code from sp::CRS("+init=epsg:26918") I have to use sp::CRS(SRS_string = "EPSG:26918"). That is an example among others. As a consequence I now skip a lot of tests if rgdal::new_proj_and_gdal() == FALSE and I removed workarounds to support old versions. Everything has been replaced by sf based code.

In lidR the occurence of rgdal:: are now

• rgdal::set_thin_PROJ6_warnings(TRUE) in unit tests setup
• Several rgdal::new_proj_and_gdal() to skip some unit tests + in the core code to fail nicely.

I'm expecting to get things cleaner when my classes will no longer inherit Spatial and I will be able to store crs instead of CRS

[rsbivand]

When you are closer, pleae let me know, and I'll try to test on a system with the same GEOS, PROJ and GDAL versions as used by CRAN Solaris. They are similar to those found on LTS CentOS and Debian typically used in large labs.

[Jean-Romain]

Well the current code passes CRAN check and the 1345 unit tests on github actions for R-devel and R-release. So you can test it right now to check if it passes with other versions of GEOS, PROJ, GDAL. Thanks

[rsbivand]

This is the check directory without the source and installed package:
lidR.Rcheck_GEOS_3.6.4_GDAL_2.2.4_PROJ_5.2.0.zip

There is one test error in test-projection.R:91:2, which when run live says:

> library(lidR)
Loading required package: raster
Loading required package: sp
lidR 3.1.5 using 2 threads. Help on <gis.stackexchange.com>. Bug report on <github.com/Jean-Romain/lidR>.
> random_10_points <- lidR:::generate_las(10)
> las = random_10_points
> epsg(las) <- 200800
Error: Invalid epsg code
> wkt(las) <- "INVALID"
Error: This object is not in a format that supports WKT CRS. Use an EPSG code.

so I think the error messages differ slightly from those in the test. I added skip(!rgdal::new_proj_and_gdal()) after the test_that( on line 89, and with that edit, it now passes. This is so small a change that it isn't worth a PR. With the revision, the check directory is:
lidR_after_edit.Rcheck.zip

[Jean-Romain]

Thank you for reporting. I identified the problem.

wkt(las) <- "INVALID"

Was expected to trigger an error saying that "INVALID" is not a valid WKT string but because there is no longer any support of "WKT" string conversion with old software instrumentation it returns sp::CRS() without error and fall to the next error saying that it is not possible to store a WKT string in this point cloud that is in a format that expect EPSG code.

It can only be a problem with old GDAL/PROJ + LAS format with WTK support. In this case "INVALID" will be actually written in the LAS file. Consequently I added an error instead.

test_that("Set an invalid code or WKT fails", {
  expect_error({ epsg(las) <- 200800 }, "Invalid epsg")

  if (rgdal::new_proj_and_gdal())
    expect_error({ wkt(las) <- "INVALID" }, "Invalid WKT")
  else
    expect_error({ wkt(las) <- "INVALID" }, "WKT strings are no longer supported in lidR with old versions of GDAL/PROJ")
})

[rsbivand]

Good, check gets past that error now, but fails at the line 100 on wkt2CRS():
lidR _test-projection.R_100_2.Rcheck.zip

[Jean-Romain]

Ho crap ! Of course it has cascading consequences... It's such a pitty it is so complex to test myself... anyway. I fixed the 3 errors again by skipping stuff for old GDAL/PROJ

By the way, when I will move to full sf to handle CRS stuff I guess I'll have the same troubles. What is the sf way to check old/new GDAL/PROJ ?

[rsbivand]

Next failure in test-las_check:83:
lidR_test-las_check.R_83_3.Rcheck.zip
For version checking, maybe see r-spatial/sf#1653, where we got to:

compareVersion(sf::sf_extSoftVersion()["PROJ"], "6.0.0") < 0 && compareVersion(sf::sf_extSoftVersion()["GDAL"], "3.0.0")

or something similar - more details in issue, which was about GEOS, but the same principles apply.

[Jean-Romain]

Always the same, trying to put a WKT string with old GDAL/PROJ fails. I greped occurrences of wkt()<- function call it should be ok now. Thank you for your time.

[rsbivand]

And line 104 too, I think it is wkt2CRS() that is where the message change bites, it being called in las_check.R once and projection.R twice as far as I can grep.

lidR_test-las_check.R_104.Rcheck.zip

[Jean-Romain]

Well, cascading effects go more in depth than expected. As I said earlier I was not sure of what I was doing. I could fix it with a skip but at this stage I want some feature to work even with old GDAL/PROJ

[Jean-Romain]

I tried to remove gdal and proj to install older versions but the oldest I can get in repos are gdal 3.0.1 and proj 6.3.0. So even making efforts it is impossible to reproduce 😞

[rsbivand]

Yes, I know, I can carry on merging from your upstream, generating the revised package tarball and running checks on my laptop, which has the old external software versions, at least for another 10 days. I think sf can generate WKT (not WKT2) on older PROJ/GDAL, and can turn WKT to a "CRS" object:

o <- st_crs(27700)
as(st_crs(o$wkt), "CRS")

might work for wkt2CRS(). in both pre-GDAL3/PROJ6 and GDAL3/PROJ6 settings.

Tomorrow I'm looking after grandchildren, so not so responsive ...

[Jean-Romain]

That seems the very best options if it works. I'm now using as(st_crs(wkt), "CRS") + st_crs(CRS)$epsg. This simplified a lot the code and I no longer have any occurrence of rgdal::new_gdal_and_proj() in the core code. I hope this work on older PROJ/GDAL.

However I still skip some unit tests because many of the test are written with sp::CRS(SRS_string = "EPSG:xxxxx") and I won't rewrite all tests.

Don't worry with the timing I also have kids to take care of 😉. Thanks for you help.

[rsbivand]

Now looks as though check is the same for PROJ 5.2.0 and GDAL 2.2.4 as PROJ 8.1.1 and GDAL 3.3.2 (latest). Using sf anyway bars GDAL < 2, but that is acceptable, only rgdal still supports GDAL 1, because it was started such a long time ago, and we haven't heard any complaints with regard to sf.

lidR_2109040930.Rcheck.zip

[Jean-Romain]

Very good! So we no longer have any occurrence of rgdal:: in the core code. Only set_thin_PROJ6_warnings and new_gdal_and_proj in unit tests. I moved rgdal to Suggest and updated the doc and error messages in few places so there is no longer any occurrence at all of the word rgdal in the entire core code.

I'll start a project of rebuilding my S4 classes without inheriting Spatial and I will try to maintain maximum backward compatibility. I will stop using sp internally and plan to find only some occurrence of sf::as_Spatial() for backward compatibility (this project is already started actually). I will eventually stop depending on raster and sp by 2024.

I'm closing this issue that is resolved. Thank you for your time

[annamuji]

I think I may have a similar problem. From a database df I want to transform the lon-lat data to utm. I am quite new to R and this line originally worked:

sp <- SpatialPoints(df[,c("longitude","latitude")],proj4string=CRS("+init=epsg:4326 +proj=longlat 
+ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 "))

however, when I transformed it into utm data using:

sp <- SpatialPoints(df[,c("longitude","latitude")],proj4string=CRS("+init=epsg:4326 +proj=longlat 
+ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 "))

an error appeared saying that I needed rgdal package.

After struggling with it I installed the libraries using the terminal (UBUNTU 18.04) but now I have the following output when I run the first order:

NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
Error in CRS("+init=epsg:4326 +proj=longlat \n+ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 ") : 
  no system list, errno: 2

Can anyone help me please?

Thank you

[rsbivand]

@annamuji this has nothing at all to do with this issue. Please do not hijack development issues to ask user questions. Please do use the R-sig-geo list after subscribing https://stat.ethz.ch/mailman/listinfo/R-SIG-Geo (this is where a motivated answer can be expected), or StackOverflow.

Please understand that you need to do much more yourself if you choose not to use CRAN binaries for Windows or MacOS; on other platforms, you have to install packages from source, satisfying their external requirements listed on their CRAN landing pages, for rgdal https://cran.r-project.org/package=rgdal. Your installation of rgdal is completely broken, because your system does not adequately satisfy the requirements. Very likely you have installed multiple versions of GDAL and/or PROJ, and only you or your local support can sort this out. See also https://github.com/r-spatial/sf/#multiple-gdal-geos-andor-proj-versions-on-your-system. Preferably use sf not *sp.

[annamuji]

ok. Thank you! sorry for hijacking.

I'll try to follow the instructions and if I don't find the solution I'll post it on the correct forum.