Open theroggy opened 4 months ago
@edzer are you sure the way this is handled in {sf} is actually correct? Meaning it correctly roundtrips and follows the spec?
@martinfleis which spec are you referring to?
@edzer I think I meant GeoPackage spec. From what I understand when @theroggy explained that, undefined SRS shall be saved as -1. But if we replicate the WKT from the sf issue linked above, it is not the case and it is saved as a custom SRS instead with ID 10000 or something like that (@theroggy please correct me here).
IIRC, the problem with the undefined SRS in GeoPackages is that they are assumed (and, after roundtrip read as) geographic coordinates with undefined datum. If you want to specify they are Cartesian but otherwise unknown, you can use the WKT
ENGCRS["Undefined Cartesian SRS with unknown unit",
EDATUM["Unknown engineering datum"],
CS[Cartesian, 2],
AXIS["X",unspecified,ORDER[1],LENGTHUNIT["unknown",0]],
AXIS["Y",unspecified,ORDER[2],LENGTHUNIT["unknown",0]]
]
which is what R package sf
writes to geopackages when the CRS is NA: https://github.com/r-spatial/sf/blob/main/R/read.R#L548
the undefined SRS in GeoPackages is that they are assumed (and, after roundtrip read as) geographic coordinates with undefined datum. If you want to specify they are Cartesian but otherwise unknown, you can use the WKT
Undefined in sense of SRS=0 is assumed geographic but SRS=-1 is assumed Cartesian. Citing Requirement 11 from the version 1.4 of the spec:
The record with an srs_id of -1 SHALL be used for undefined Cartesian coordinate reference systems. The record with an srs_id of 0 SHALL be used for undefined geographic coordinate reference systems.
That is where I think that the sf solution with custom WKT is not exactly reflecting the spec.
GDAL correctly understands -1 as Cartesian but, for unknown reasons, sets units to meters here. We discusses on the dev call earlier today that this should probably be fixed to resolve as Undefined Cartesian SRS with unknown units rather than with meters. Then setting SRS to -1 should do what we want it to without specifying custom WKT.
So, maybe raise an issue with GDAL?
As planned in the community meeting and suggested by @edzer , I opened an issue to discuss this further in GDAL: https://github.com/OSGeo/gdal/issues/9580
Just checking back on the status of this PR; it looks like this is potentially waiting on upstream changes to GDAL possibly in 3.9? But then what is the path forward on this for GDAL < 3.9?
Questions/remarks:
resolves #367