GDAL 3.3.2 change

[rsbivand]

See OSGeo/gdal#4047 and https://lists.osgeo.org/pipermail/gdal-dev/2021-August/054600.html, probably enforcing 4326 axis order. See also advice to use OSRFindMatches() rather than AutoIdentifyEPSG().

[edzer]

Looks like I was wrong in #489?

[edzer]

sf checks cleanly against this RC, but when writing & reading a gpkg with a missing CRS, I got (from ogrinfo) this:

Layer SRS WKT:
(unknown)

before, and now I get

Layer SRS WKT:
GEOGCRS["Undefined geographic SRS",
    DATUM["unknown",
        ELLIPSOID["unknown",6378137,298.257223563,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 2,1

which, when reading back in R gives

Geodetic CRS:  Undefined geographic SRS

rather than what was before

CRS:           NA

That's quite a change, right?

[edzer]

FindMatches btw writes a cache to ~/.gdal/$(GDAL_VERSION_MAJOR).$(GDAL_VERSION_MINOR)/srs_cache

[rsbivand]

After revdeps, I think no packages roundtrip writing/reading files, so don't test them. The modal failure is that "+proj=utm" without a zone number now fails, and a number of packages (not just NetLogoR) have used it as a way of declaring an unspecified planar CRS.

[edzer]

Thanks! Good to know; we now know that that should be st_crs("LOCAL_CS[\"Undefined Cartesian SRS\"]"). We could even make this the default in sf, rather than st_crs(NA), since we treat NA CRS objects as having Cartesian coordinates with E/N axis order.

[rsbivand]

> library(sf)
Linking to GEOS 3.10.0dev, GDAL 3.3.2, PROJ 8.1.1
> as(st_crs("LOCAL_CS[\"Undefined Cartesian SRS\"]"), "CRS")
Coordinate Reference System:
Deprecated Proj.4 representation: NA 
WKT2 2019 representation:
ENGCRS["Undefined Cartesian SRS",
    EDATUM[""],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]] 
> sp::is.projected(as(st_crs("LOCAL_CS[\"Undefined Cartesian SRS\"]"), "CRS"))
[1] FALSE

[rsbivand]

With an as-yet uncommitted change in rgdal:

> library(sf)
Linking to GEOS 3.10.0dev, GDAL 3.3.2, PROJ 8.1.1
> obj <- st_crs("LOCAL_CS[\"Undefined Cartesian SRS\"]")
> sp::CRS(obj$wkt)
proj_get_ellipsoid: CRS has no geodetic CRS
proj_get_ellipsoid: Object is not a CRS or GeodeticReferenceFrame
proj_as_proj_string: Object type not exportable to PROJ
proj_get_ellipsoid: CRS has no geodetic CRS
proj_get_ellipsoid: Object is not a CRS or GeodeticReferenceFrame
Coordinate Reference System:
Deprecated Proj.4 representation: NA 
WKT2 2019 representation:
ENGCRS["Undefined Cartesian SRS",
    EDATUM[""],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]] 
Warning message:
In doTryCatch(return(expr), name, parentenv, handler) :
  export to PROJ failed: Unknown error (code 4096)

For sp, could we add an rda in data/ with this CRS object, and direct users to that rather than rolling their own?

[edzer]

Yes, but maybe that is taking things too fast?

[rsbivand]

Perhaps, let's see. A resolution in sf is more important.

[rsbivand]

There are problems with using $epsg to get the EPSG code: r-lidar/lidR#466 (comment).

[rsbivand]

The +proj=utm" problem is caused by OSGeo/PROJ#2672 from OSGeo/PROJ#2671 and applies in PROJ >= 8, probably backported to 8.0.1 and applies in 8.1.*. Previously a non-existing zone 0 was generated.