Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

raster crs representation #85

Closed
florisvdh opened this issue Sep 28, 2020 · 3 comments
Closed

raster crs representation #85

florisvdh opened this issue Sep 28, 2020 · 3 comments
Assignees
Labels
enhancement New feature or request solved_in_0.4.0

Comments

@florisvdh
Copy link
Member

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 florisvdh added the enhancement New feature or request label Sep 28, 2020
@florisvdh florisvdh self-assigned this Sep 28, 2020
@florisvdh
Copy link
Member Author

florisvdh commented Dec 9, 2020

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
Copy link
Member Author

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
Copy link
Member Author

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request solved_in_0.4.0
Projects
None yet
Development

No branches or pull requests

1 participant