search

inbo / n2khab

R package with preprocessing functions and standard reference data for Flemish Natura 2000 (N2K) habitat (HAB) analyses
https://inbo.github.io/n2khab
GNU General Public License v3.0
2 stars 1 forks source link

raster crs representation #85

Closed florisvdh closed 3 years ago

florisvdh commented 4 years ago

Follow-up to #67 and #84.

CRS representation is to be revisited once raster from CRAN accepts simple EPSG-code as input for crs<- (as on GitHub), i.e. with raster>=3.3-16. This would allow dropping the Suggests: sp part.

Maybe, explicit CRS setting can then be just dropped as well, insisting on PROJ 6+/GDAL 3, on condition this omission causes no further clashes with other object's CRSes. The same may be considered for sf data sources at that time.

florisvdh commented 3 years ago

Maybe, explicit CRS setting can then be just dropped as well, insisting on PROJ 6+/GDAL 3

Nope, it will always clash for datasets constructed by ESRI software. The ESRI CRS specification (which is abundant in Geopunt data sources) of its 'Belge 1972 / Belgian Lambert 72' CRS, which it says to be 'Well known ID 31370', is syntactically different from the official EPSG specification of EPSG:31370. Notably, it specifies 'false easting' and 'false northing' with a precision of 10-5 m and 10-4 m respectively (EPSG: precision 10-3 m), and it switches 1st and 2nd standard parallels.

This results in negligable actual coordinate differences (supporting the approach of overriding the CRS with official EPSG:31370), but as the ESRI 'Belge 1972 / Belgian Lambert 72' CRS is technically different from official EPSG:31370 'Belge 1972 / Belgian Lambert 72' CRS, despite ESRI's reference to it, the ESRI variant results in clashes when confronted with actual EPSG:31370 (this can be observed e.g. in R and in GRASS). Overriding and setting as EPSG:31370 seems the best option.

Also the Geopunt metadata state that the CRS is EPSG:31370 (example) while technically it is different from EPSG:31370.

UPDATE: the ESRI version of 'Belge 1972 / Belgian Lambert 72' is known by PROJ as ESRI:102499 = "Belge_Lambert_1972_bad_FE_FN" (referring to bad false easting & northing). Compare below WKT2 strings.

projinfo EPSG:31370 -o WKT2:2019 ``` WKT2:2019 string: PROJCRS["Belge 1972 / Belgian Lambert 72", BASEGEOGCRS["Belge 1972", DATUM["Reseau National Belge 1972", ELLIPSOID["International 1924",6378388,297, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4313]], CONVERSION["Belgian Lambert 72", METHOD["Lambert Conic Conformal (2SP)", ID["EPSG",9802]], PARAMETER["Latitude of false origin",90, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8821]], PARAMETER["Longitude of false origin",4.36748666666667, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8822]], PARAMETER["Latitude of 1st standard parallel",51.1666672333333, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8823]], PARAMETER["Latitude of 2nd standard parallel",49.8333339, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8824]], PARAMETER["Easting at false origin",150000.013, LENGTHUNIT["metre",1], ID["EPSG",8826]], PARAMETER["Northing at false origin",5400088.438, LENGTHUNIT["metre",1], ID["EPSG",8827]]], CS[Cartesian,2], AXIS["easting (X)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["northing (Y)",north, ORDER[2], LENGTHUNIT["metre",1]], USAGE[ SCOPE["unknown"], AREA["Belgium - onshore"], BBOX[49.5,2.5,51.51,6.4]], ID["EPSG",31370]] ```
projinfo ESRI:102499 -o WKT2:2019 ``` WKT2:2019 string: PROJCRS["Belge_Lambert_1972_bad_FE_FN", BASEGEOGCRS["Belge 1972", DATUM["Reseau National Belge 1972", ELLIPSOID["International 1924",6378388,297, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4313]], CONVERSION["Belge_Lambert_1972_bad_FE_FN", METHOD["Lambert Conic Conformal (2SP)", ID["EPSG",9802]], PARAMETER["Latitude of false origin",90, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8821]], PARAMETER["Longitude of false origin",4.36748666666667, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8822]], PARAMETER["Latitude of 1st standard parallel",49.8333339, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8823]], PARAMETER["Latitude of 2nd standard parallel",51.1666672333333, ANGLEUNIT["degree",0.0174532925199433], ID["EPSG",8824]], PARAMETER["Easting at false origin",150000.01256, LENGTHUNIT["metre",1], ID["EPSG",8826]], PARAMETER["Northing at false origin",5400088.4378, LENGTHUNIT["metre",1], ID["EPSG",8827]]], CS[Cartesian,2], AXIS["(E)",east, ORDER[1], LENGTHUNIT["metre",1]], AXIS["(N)",north, ORDER[2], LENGTHUNIT["metre",1]], USAGE[ SCOPE["unknown"], AREA["Belgium - onshore"], BBOX[49.5,2.5,51.51,6.4]], ID["ESRI",102499]] ```
florisvdh commented 3 years ago

This would allow dropping the Suggests: sp part.

Meanwhile, it is advised to import or suggest rgdal for functions that use raster. While raster does not import it, raster does use rgdal (it's in its 'suggests') in functions used by n2khab. And this makes renv::init() and renv:snapshot() standard implementations incomplete, i.e. missing rgdal (example fix here).

florisvdh commented 3 years ago

this makes renv::init() and renv:snapshot() standard implementations incomplete, i.e. missing rgdal

From the documentation of renv::dependencies() and from tests, it can be concluded that rgdal will only be added to renv.lock by default when either a package (that is explicitly loaded in a project) imports rgdal (as defined in its DESCRIPTION file), or library(rgdal) explicitly appears in a project script.

Both sp and raster have rgdal in their Suggests section. So the challenge lies entirely with renv / the renv user; we should not try to handle this in n2khab.

So we won't follow:

Meanwhile, it is advised to import or suggest rgdal for functions that use raster.