lidR vulnerable to forthcoming changes in sp and rgdal #299
Comments
|
Thank you for reporting. It does not bring much trouble and it is easy to fix in my case. I just need to find how to install gdal3 and proj6. |
Jean-Romain
added
Bug
|
Major problems remain as of sp 1.4-1, rgdal 1.5-6, sf 0.9-0: Taking test-LAS, line 85: The internals of |
|
I'm able to reproduce and I fixed those issues that were not critical. However I also have 500+ warnings because the following now throw a warnings: wkt <- rgdal::showWKT("+init=epsg:2008")
p4 <- rgdal::showP4(wkt)
sp::CRS(p4)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
#> datum Unknown based on Clarke 1866 ellipsoid in CRS definition
#> CRS arguments:
#> +proj=tmerc +lat_0=0 +lon_0=-55.5 +k=0.9999 +x_0=304800 +y_0=0
#> +ellps=clrk66 +units=m +no_defs
wkt <- rgdal::showWKT("+init=epsg:2949")
p4 <- rgdal::showP4(wkt)
sp::CRS(p4)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
#> datum Unknown based on GRS80 ellipsoid in CRS definition
#> CRS arguments:
#> +proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 +x_0=304800 +y_0=0
#> +ellps=GRS80 +units=m +no_defsI can fixed ~250 of them with I need to be able to transform an epsg code into a proj4 that does not trigger those warnings. |
|
Thanks for moving on this. See if |
|
I tested rgdal::showSRID("EPSG:2949", "PROJ")
#> Warning in rgdal::showSRID("EPSG:2949", "PROJ"): Discarded datum
#> NAD83_Canadian_Spatial_Reference_System in CRS definition
#> [1] "+proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +units=m +no_defs"And rgdal::set_thin_PROJ6_warnings(TRUE)
rgdal::showSRID("EPSG:2949", "PROJ")
#> Warning in rgdal::showSRID("EPSG:2949", "PROJ"): PROJ6/GDAL3 PROJ string degradation in workflow
#> repeated warnings suppressed
#> Discarded datum NAD83_Canadian_Spatial_Reference_System in CRS definition
#> [1] "+proj=tmerc +lat_0=0 +lon_0=-70.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +units=m +no_defs"
Ok but how we do that? |
|
The warnings are just warnings, they are meant to alert the user that an important part of the Proj4 representation ( |
|
And moving away from sp will not help. terra uses simple CRS, but they will break soon, and need fixing, since terra uses GDAL and GDAL uses WKT2, and will make Proj4 strings defunct. If you have a simple test-bed script for your classes/workflow, I can look at how it might be developed to use WKT2. |
Yes for sure. But on CRAN a non fixable warning is likely to be an archived package...
So, now the errors are fixed can I install back rgdal 1.4-8 for not being annoyed with all these warnings and wait for the CRAN release safely? If I'm understanding well,
That is basically what I doing as we already discussed. I'm reading an epsg code or a WKT string (rarely) and I'm making a CRS compatible with the geospatial ecosystem (i.e. a proj4). Then, the CRS is propagated to the outputs usually |
There are 4 classes:
setClass(
Class = "LAS", contains = "Spatial",
representation(data = "data.table", header = "LASheader"))
setClass(
Class = "LAScatalog",
contains = "SpatialPolygonsDataFrame",
representation(
chunk_options = "list",
processing_options = "list",
output_options = "list",
input_options = "list" ))Below some simple tests related to CRS that can be performed. If library(lidR)
library(testthat)
# Using rgdal 1.5-6
# Make a LAS
LASfile <- system.file("extdata", "Megaplot.laz", package = "lidR")
las <- readLAS(LASfile)
# Make a RasterLayer
bbox <- extent(las)
r <- raster(bbox)
r[] <- sample(1:ncell(r))
# Make a SpatialPointsDataFrame
spdf <- as.spatial(lasfilterground(las))
# Tests set 1
expect_equal(raster::projection(las), "+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs")
expect_warning(raster::projection(r) <- raster::projection(las), NA)
expect_warning(raster::projection(spdf) <- raster::projection(las), NA)
expect_warning(sp::proj4string(spdf) <- sp::proj4string(las), NA)
# Test set 2 - change the espg code for a dummy one that is known to trigger a warning
epsg(las) <- 2008
expect_equal(raster::projection(las), "+proj=tmerc +lat_0=0 +lon_0=-55.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=clrk66 +units=m +no_defs")
expect_equal(las@header@VLR[["GeoKeyDirectoryTag"]][["tags"]][[2]][["value offset"]], 2008)
expect_warning(raster::projection(r) <- raster::projection(las), NA)
expect_warning(raster::projection(spdf) <- raster::projection(las), NA)
expect_warning(sp::proj4string(spdf) <- sp::proj4string(las), NA)
# Test set 3 - set a projection to a LAS object
projection(las) <- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"
expect_equal(epsg(las), 32617)
# Test set 4 - LAScatalog get attributed a CRS
ctg = readLAScatalog(LASfile)
expect_equal(sp::proj4string(ctg), "+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs") |
|
Warnings are not CRAN warnings. R CMD check just runs examples, which hard-fail only when they hit a Please (re-)read https://www.r-spatial.org/r/2020/03/17/wkt.html to see that is being done in sp. The useless legacy Proj4 string is retained in the would simply be Please re-make the example without raster - it has not adapted to PROJ changes yet, and without testthat, just put the expected value in a comment. Probably all testing of expected values of CRS objects should wait until the situation has stabilised, say by mid-2022 (sorry, but everything has been disorganised since 2017 already, so mid-2022 is a sober estimate. |
|
Thank you for the link, this clarifies a lot the question. I'll review my code to check if I need to change something. Especially I may have used bad assignment practices in the past. The current state of the code should be free of such error but I can't swear that unit tests are all good.
This is what I did so far but I changed that recently because invalid epsg codes now print a message in console ( crs = tryCatch(sp::CRS("+init=epsg:1234"), error = function(x) {sp::CRS()})
#> proj_create: crs not found
sp::CRS(SRS_string='EPSG:4326')
#> Error in sp::CRS(SRS_string = "EPSG:4326") :
#> no arguments in initialization listBelow a simplified set of tests. I understand I should not test the proj strings. Most of the time my tests actually aims to check is (1) the CRS was propagated to the output and (2) the CRS was properly read from files. I will see if I can simplify the tests to rely less on a strict string matching. library(lidR) # Using rgdal 1.5-6
# Make a LAS
LASfile <- system.file("extdata", "Megaplot.laz", package = "lidR")
las <- readLAS(LASfile)
sp::proj4string(las) # "+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs"
epsg(las) <- 2008
sp::proj4string(las) # "+proj=tmerc +lat_0=0 +lon_0=-55.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=clrk66 +units=m +no_defs"
epsg(las) # 2008
projection(las) <- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"
epsg(las) # 32617
ctg = readLAScatalog(LASfile)
sp::proj4string(ctg) # "+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs" |
|
Thanks, useful. I've edited and simplified error throwing and reporting in rgdal, in SVN revision 951 (R-forge) I see: So with R-Forge rgdal >= rev 951, you might do this: or maybe: to alert a user to the invalid string, here the EPSG code, but if you were passing a WKT, an invalid set of specifications (then via the new You will only see changes in behaviour with sp 1.4-1 (CRAN) and rgdal from R-forge: install.packages("rgdal", repos="http://R-Forge.R-project.org"), installed from source and build on PROJ >= 6 and GDAL >= 3. I'm usure which you have. I see: so the so the datum is dropped in GDAL from the EPSG definition. I think the logic in your |
|
When reading a las file the only available information is the epsg. So I trust the epsg. No choice. For newer formats we have WKT and in this case epsg codes are not involved. But both espg or WKT are transformed to proj to go to into a
I'm up-to-date with all geosptial stuff. I have sp 1.4.1, rgdal 1.5-6 from rforge, GDAL 3, PROJ 7 but yet I can't reproduce your warnings
Well setMethod("epsg<-", "LAS", function(object, value)
{
crs <- sp::CRS(paste0("+init=epsg:", value))
epsg(object@header) <- value # set the epsg code in the header so we can write a georeferenced files
raster::projection(object) <- crs # set a CRS compatible with R ecosytem.
return(object)
})
setMethod("wkt<-", "LAS", function(object, value)
{
proj <- rgdal::showP4(value) # can be replaced by sp::CRS(SRS_string = value)
crs <- sp::CRS(proj)
wkt(object@header) <- value # set the WKT in the header so we can write a georeferenced file
raster::projection(object) <- crs # set a CRS compatible with R ecosytem.
return(object)
})By the way you should not have been able to do Actually the more complex function is projection(las) <- projection(object)And projection is more complex (but yet easy to understand). Here a sightly simplified version. It looks at the proj, converts it to WKT if the format supports WKT (LAS 1.4) or it retrieves the EPSG if the file format supports only epsg. There is room for improvement with the new # x is a LAS
# value is a CRS
setMethod("projection<-", "LAS", function(x, value)
{
proj4 <- value@projargs
# the object is in a format that supports WKT (LAS 1.4)
if (use_wktcs(x))
{
wkt <- rgdal::showWKT(proj4) # can now be replaced by comment(value)
wkt(x@header) <- wkt
x@proj4string <- sp::CRS(proj4)
return(x)
}
# the object is in a format that supports only EPSG code
# (might be simplified by using the WKT in `comment(value)`)
else
{
# Extract epsg code (if any) using +init=epsg:xxxx part
epsg <- sub("\\+init=epsg:(\\d+).*", "\\1", proj4)
epsg <- suppressWarnings(as.integer(epsg))
# We were not able to extract the epsg code, we can retrieve it with rgdal
if (is.na(epsg)) epsg <- rgdal::showEPSG(proj4)
# We are still unable to find an epsg code
if (epsg == "OGRERR_UNSUPPORTED_SRS")
warning("EPSG code not found: header not updated. Try to use the function epsg() manually to ensure CRS will be written in file.", call. = FALSE)
# We found the epsg code
else
epsg(x@header) <- epsg
x@proj4string <- value
return(x)
}
}) |
Running revdep checks for current rgdal on R-Forge - see:
https://stat.ethz.ch/pipermail/r-sig-geo/2019-November/027801.html
shows the errors in the attached test failure log, related to use of PROJ&/GDAL3
and required changes to sp and rgdal. If useful find a regerence to a docker
image in this thread:
r-spatial/discuss#28
Changes will occur quite fast, and packages need to be prepared.
The text was updated successfully, but these errors were encountered: