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.

Sign up for GitHub

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

Jump to bottom

Forthcoming PROJ 9.6.0 breaks geofi #52

Open
rsbivand opened this issue Mar 14, 2025 · 6 comments
Open

Forthcoming PROJ 9.6.0 breaks geofi #52

rsbivand opened this issue Mar 14, 2025 · 6 comments

Comments

@rsbivand
Copy link

As forthcoming PROJ 9.6.0 announced here https://lists.osgeo.org/pipermail/proj/2025-March/011738.html adds NKG transformations: Add support for EUREF-FIN in Finish transformations (#4399) in the NEWS section, geofi now breaks at line 152 in geofi_spatial_analysis.Rmd, where sf::st_intersection fails with differing CRS:

Error in geos_op2_geom("intersection", x, y, ...) : 
  st_crs(x) == st_crs(y) is not TRUE
> st_crs(muni)
Coordinate Reference System:
  User input: EPSG:3067 
  wkt:
PROJCRS["EUREF-FIN / TM35FIN(E,N)",
    BASEGEOGCRS["EUREF-FIN",
        DATUM["EUREF-FIN",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ANCHOREPOCH[1997]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",10690]],
    CONVERSION["TM35FIN",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",27,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Finland - onshore and offshore."],
        BBOX[58.84,19.08,70.09,31.59]],
    USAGE[
        SCOPE["Pan-European conformal mapping at scales larger than 1:500,000."],
        AREA["Finland - between 24°E and 30°E, onshore and offshore."],
        BBOX[59.64,23.99,70.09,30]],
    ID["EPSG",3067]]
> st_crs(bounding_box_point)
Coordinate Reference System:
  User input: EPSG:3067 
  wkt:
PROJCRS["ETRS89 / TM35FIN(E,N)",
    BASEGEOGCRS["ETRS89",
        ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
            MEMBER["European Terrestrial Reference Frame 1989"],
            MEMBER["European Terrestrial Reference Frame 1990"],
            MEMBER["European Terrestrial Reference Frame 1991"],
            MEMBER["European Terrestrial Reference Frame 1992"],
            MEMBER["European Terrestrial Reference Frame 1993"],
            MEMBER["European Terrestrial Reference Frame 1994"],
            MEMBER["European Terrestrial Reference Frame 1996"],
            MEMBER["European Terrestrial Reference Frame 1997"],
            MEMBER["European Terrestrial Reference Frame 2000"],
            MEMBER["European Terrestrial Reference Frame 2005"],
            MEMBER["European Terrestrial Reference Frame 2014"],
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[0.1]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["TM35FIN",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",27,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Finland - onshore and offshore."],
        BBOX[58.84,19.08,70.09,31.59]],
    ID["EPSG",3067]]
> st_crs(muni) == st_crs(helsinki)
[1] TRUE
> st_crs(muni) == st_crs(bounding_box_point)
[1] FALSE

There is a further error at line 279.

I guess that conditioning on the PROJ version will be needed, using for example package_version(sf::sf_extSoftVersion()["PROJ"]).
@rsbivand
Copy link
Author

rsbivand commented Mar 14, 2025 •
edited
Loading

The problem arises from "canned" sf objects such as municipality_central_localities giving:

> st_crs(muni) == st_crs(point)
[1] FALSE

You may duplicate these for PROJ < 9.6.0 and PROJ >= 9.6.0 perhaps, or for the vignette now set the CRS of point (line 63) to that of muni, which is generated on-the-fly on reading. Note that EPSG definitions are not fixed across releases of the database, 9.6.0 EPSG v12.004, 9.5.1 EPSG v11.022. See also https://lists.osgeo.org/pipermail/proj/2025-March/011772.html.

@rsbivand
Copy link
Author

rsbivand commented Mar 14, 2025 •
edited
Loading

@muuankarski If CRAN updates to PROJ 9.6.0 when Debian or other check platform does, geofi will fail on CRAN, so please try to address this before that happens. When you have a draft fix but don't have 9.6.0, I can run the draft on a system with 9.6.0. And see: https://lists.osgeo.org/pipermail/proj/2025-March/011774.html.

@pitkant
Copy link
Member

pitkant commented Mar 31, 2025

@muuankarski If I understood the problem correctly I think this could be fixed by implementing the wording from the vignette "All the data you can obtain using geofi is transformed automatically into EPSG:3067." from line 56:

geofi/vignettes/geofi_spatial_analysis.Rmd

Line 56 in f7fa1d9

When using spatial data in R it necessary to have all data in same coordinate reference system (CRS). You can check the CRS of you `sf`-object with `sf::st_crs()`-function. All the data you can obtain using `geofi` is transformed automatically into `EPSG:3067`. Most of spatial data providers in Finland provide their data in the same CRS.

In get_municipalities.R EPSG:3067 is used only if the sf::st_crs(sf_obj) is NA.

geofi/R/get_municipalities.R

Lines 51 to 54 in f7fa1d9

if (is.na(sf::st_crs(sf_obj))) {
warning("Coercing CRS to epsg:3067 (ETRS89 / TM35FIN)", call. = FALSE)
sf::st_crs(sf_obj) <- 3067
}

Maybe it was interpreted as NA before but now it isn't when EUREF-FIN is supported? Anyway if it is changed so that 3067 is used in all cases as default then the problem would go away?

@rsbivand
Copy link
Author

@pitkant I do not think this is the case. Users of PROJ < 9.6.0 will see a different EPSG:3067 than users of PROJ >= 9.6.0 anyway, and the code error occurred because comparing the stored data values of the CRS created with PROJ < 9.6.0 are not equal to the CRS of the object created by get_municipalities (the CRS created on the fly from PROJ >= 9.6.0), so sf topological operations and predicates fail. https://lists.osgeo.org/pipermail/proj/2025-March/011787.html explains why the EPSG definition for 3067 was changed.

@pitkant
Copy link
Member

pitkant commented Mar 31, 2025

@rsbivand thank you for clarifying! It's just curious to me that bounding_box_point, which is created by subsetting the muni object and then turned into a bbox, would then have different CRS than the original muni object. From that perspective both should have a CRS created on the fly, or am I missing something?

@muuankarski
Copy link
Collaborator

muuankarski commented Mar 31, 2025 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@muuankarski @rsbivand @pitkant