search

OSGeo / PROJ

PROJ - Cartographic Projections and Coordinate Transformations Library
https://proj.org
Other
1.75k stars 791 forks source link

Export error from Geographic 3D CRS in WKT1_ESRI for Non-Earth Parameters #2757

Closed dpmayerUSGS closed 3 years ago

dpmayerUSGS commented 3 years ago

Example of problem

I am using ogr2ogr to convert a vector VRT that describes 3D point data from a CSV to ESRI Shapefile format, using a 3D CRS, but I keep getting an error from PROJ that it can't export the CRS as WKT1_ESRI. The CRS I'm trying to use is for Mars.

Contents of VRT file:

<OGRVRTDataSource>
<OGRVRTLayer name="mini_llt">
<SrcDataSource>mini_llt.csv</SrcDataSource>
<LayerSRS>mini_llt.prj</LayerSRS>
<GeometryType>wkbPoint25D</GeometryType>
<GeometryField encoding="PointFromColumns" x="LONG_EAST" y="LAT_NORTH" z="TOPOGRAPHY" />
<Field name="LONG_EAST" src="LONG_EAST" type="Real" />
<Field name="LAT_NORTH" src="LAT_NORTH" type="Real" />
<Field name="TOPOGRAPHY" src="TOPOGRAPHY" type="Real" />
</OGRVRTLayer>
</OGRVRTDataSource>

Contents of mini_llt.prj:

GEOGCS["GCS_Mars_Sphere_2000",
       DATUM["D_Mars_Sphere_2000",
         SPHEROID["Mars_Sphere_2000",3396190.0,0.0]],
       PRIMEM["Reference Meridian",0.0],
       UNIT["Degree",0.0174532925199433]],
VERTCS["Mars_Sphere_2000",
       DATUM["D_Mars_Sphere_2000",
         SPHEROID["Mars_Sphere_2000",3396190.0,0.0]],
       PARAMETER["Vertical_Shift",0.0],
       PARAMETER["Direction",1.0],UNIT["Meter",1.0]]

projinfo command:

projinfo "$(cat mini_llt.prj)" -o WKT1:ESRI
> WKT1:ESRI string:
> Error when exporting to WKT1:ESRI: Cannot export this Geographic 3D CRS in WKT1_ESRI

Problem description

I know that the WKT I am using is structurally correct. If I substitute in Earth-relevant values, projinfo no longer complains:

cat wgs84_esri_3d.prj 
GEOGCS["GCS_WGS_1984",
       DATUM["D_WGS_1984",
         SPHEROID["WGS_1984",6378137.0,298.257223563]],
       PRIMEM["Greenwich",0.0],
       UNIT["Degree",0.0174532925199433]],
VERTCS["WGS_1984",
       DATUM["D_WGS_1984",
         SPHEROID["WGS_1984",6378137.0,298.257223563]],
       PARAMETER["Vertical_Shift",0.0],
       PARAMETER["Direction",1.0],UNIT["Meter",1.0]]

projinfo "$(cat wgs84_esri_3d.prj)" -o WKT1:ESRI
> WKT1:ESRI string:
> GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],VERTCS["WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Meter",1.0]]

However, if create a "hybrid" WKT filled with Mars-relevant values, but /change the string name of the datum/ to "D_WGS_1984", then PROJ no longer complains.

cat hybrid_wkt1.prj
GEOGCS["GCS_Mars_Sphere_2000",
       DATUM["D_WGS_1984",
         SPHEROID["Mars_Sphere_2000",3396190.0,0.0]],
       PRIMEM["Reference Meridian",0.0],
       UNIT["Degree",0.0174532925199433]],
VERTCS["Mars_Sphere_2000",
       DATUM["D_Mars_Sphere_2000",
         SPHEROID["Mars_Sphere_2000",3396190.0,0.0]],
       PARAMETER["Vertical_Shift",0.0],
       PARAMETER["Direction",1.0],UNIT["Meter",1.0]]

projinfo "$(cat hybrid_wkt1.prj)" -o WKT1:ESRI
> WKT1:ESRI string:
> GEOGCS["GCS_Mars_Sphere_2000",DATUM["D_WGS_1984",SPHEROID["Mars_Sphere_2000",3396190.0,0.0]],PRIMEM["Reference_Meridian",0.0],UNIT["Degree",0.0174532925199433]],VERTCS["WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Meter",1.0]]

All of this behavior is mirrored when I run ogr2ogr -f "ESRI Shapefile" out.shp mini_llt.vrt (updating the VRT to point to the different .prj files, of course). I get that export error message for the "pure Mars" case, but not for the WGS84 and "hybrid" cases.

Expected Output

I would expect projinfo to export WKT1:ESRI to the screen without error when given the "pure Mars" WKT posted above. The fact that the hybrid WKT works makes me think that PROJ is scanning for the string "D_WGS_1984" in order to validate the input.

Additional Context

I started out by trying to use WKT2 that described a 3D CRS, and which validated according to gdalsrsinfo -V, but ogr2ogr still spit out the export error about WKT1_ESRI, so I switched to trying to build WKT1 in the ESRI flavor.

Environment Information

Installation method

rouault commented 3 years ago

Well, this is a bit tricky. Fundamentally, a CompoundCRS whose vertical CRS is an ellipsoidal height is not considered by ISO 19111 as a valid representation. This should rather be represented as a Geographic 3D CRS. And if you look at the WKT:2019 output of your initial CRS, you'll see that PROJ ingests it as a Geographic 3D CRS

WKT2:2019 string:
GEOGCRS["GCS_Mars_Sphere_2000",
    DATUM["D_Mars_Sphere_2000",
        ELLIPSOID["Mars_Sphere_2000",3396190,0,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]],
    PRIMEM["Reference Meridian",0,
        ANGLEUNIT["Degree",0.0174532925199433]],
    CS[ellipsoidal,3],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["Degree",0.0174532925199433]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["Degree",0.0174532925199433]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

But on export, ESRI WKT doesn't support Geographic 3D CRS, but this somewhat hacky CompoundCRS representation. Another thing to have in mind is that when PROJ ingests WKT1:ESRI, it normalizes the ESRI datum names to EPSG datum names, using the aliases available in proj.db. When exporting this Geographic 3D CRS as a CompoundCRS, it does the reverse operation and tries to find back the ESRI datum names from the EPSG official names. This works for WGS 84 because there's such alias in the database between D_WGS_1984 and "World Geodetic System 1984", but not for D_Mars_Sphere_2000 (I see there's a "Mars2000(Sphere)" however in https://github.com/Esri/projection-engine-db-doc/blob/master/csv/pe_list_datum.csv#L892, which is imported in proj.db). Hence this fails here: https://github.com/OSGeo/PROJ/blob/d41c10c9e1e787a521a3841605a2e6024649eb3f/src/iso19111/crs.cpp#L1665

Do you have access to ESRI software to check that your minillt.prj is properly supported by it ? If so, we could perhaps amend the exportAsESRIWktCompoundCRSWithEllipsoidalHeight() function so that it goes on the export when the datum name starts with D , even when it is not known as an official ESRI datum name.

rouault commented 3 years ago

Could you also try the following in ESRI software:

GEOGCS["GCS_Mars_Sphere_2000",
       DATUM["Mars_2000_(Sphere)",
         SPHEROID["Mars_Sphere_2000",3396190.0,0.0]],
       PRIMEM["Reference Meridian",0.0],
       UNIT["Degree",0.0174532925199433]],
VERTCS["Mars_Sphere_2000",
       DATUM["Mars_2000_(Sphere)",
         SPHEROID["Mars_Sphere_2000",3396190.0,0.0]],
       PARAMETER["Vertical_Shift",0.0],
       PARAMETER["Direction",1.0],UNIT["Meter",1.0]]

or possibly (spheroid name also changed to "Mars2000(Sphere)" which seems to be the now official name in Esri world)

GEOGCS["GCS_Mars_Sphere_2000",
       DATUM["Mars_2000_(Sphere)",
         SPHEROID["Mars_2000_(Sphere)",3396190.0,0.0]],
       PRIMEM["Reference Meridian",0.0],
       UNIT["Degree",0.0174532925199433]],
VERTCS["Mars_Sphere_2000",
       DATUM["Mars_2000_(Sphere)",
         SPHEROID["Mars_2000_(Sphere)",3396190.0,0.0]],
       PARAMETER["Vertical_Shift",0.0],
       PARAMETER["Direction",1.0],UNIT["Meter",1.0]]

PROJ also rejects exporting it to ESRI:WKT, because technically there's no alias between a ESRI name and EPSG name for that one. But here clearly we should also check if the datum name isn't already an official ESRI datum name.

dpmayerUSGS commented 3 years ago

I don't have immediate access to ESRI software, but I can check in the next day or so, and I will reply with my results.

Ideally, I'd like to specify some flavor of WKT2 in the .prj that I pass to ogr2ogr, but the current version doesn't seem to like that. When converting to ESRI Shapefile, will OGR only write ESRI flavor WKT? Is there a way to override this?

rouault commented 3 years ago

but the current version doesn't seem to like that. When converting to ESRI Shapefile, will OGR only write ESRI flavor WKT?

yes, this is the only flavor. Technically we could export to something else, but the aim is to be compatible with Esri software when writing shapefiles, hence we restrict to Esri WKT. If you use GeoPackage, WKT:2015 will be written when a Geogrpahic 3D CRS is used.

dpmayerUSGS commented 3 years ago

@thareUSGS was kind enough to test shapefiles using all 3 WKT descriptions (my example and the 2 suggested by @rouault ) in ArcMap and ArcPro. Unfortunately, Arc managed to load the data but complain about "Unknown Coordinate System."

However, ogr2ogr converted a VRT using the following 2D Geographic CRS into ESRI-style WKT1, and Arc opened it up just fine:

GEOGCRS["GCS_Mars_2000",
    DATUM["D_Mars_2000",
        ELLIPSOID["Mars_2000_IAU",3396190.0,169.89444722361179,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Reference_Meridian",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]]]
rouault commented 3 years ago

Unfortunately, Arc managed to load the data but complain about "Unknown Coordinate System."

And it can't make sense at all of it ? (QGIS might for example also say that it is a unknown CRS but still be able to process it)

However, ogr2ogr converted a VRT using the following 2D Geographic CRS into ESRI-style WKT1, and Arc opened it up just fine:

ok, but that's a different beast. Just a Geographic 2D CRS. It would be abusive to silently convert a Geographic 3D CRS to a 2D one.

dpmayerUSGS commented 3 years ago

And it can't make sense at all of it ?

Arc manages to render the points. @thareUSGS will have to provide details. I test all 3 shapefiles in QGIS 3.16 and found that it renders the points, and recognizes them as 3D point geometry ("PointZ").

(I only mentioned the additional test using the Geographic 2D CRS to illustrate that there isn't something else, perhaps more fundamental, wrong with the data.)

rouault commented 3 years ago

Does the below GEOGCS["GCS_WGS_1984",....],VERTCS["WGS_1984",....] work in ArcXXXX ? (I strongly suspect it does, since it appears on Internet searches)

GEOGCS["GCS_WGS_1984",
       DATUM["D_WGS_1984",
         SPHEROID["WGS_1984",6378137.0,298.257223563]],
       PRIMEM["Greenwich",0.0],
       UNIT["Degree",0.0174532925199433]],
VERTCS["WGS_1984",
       DATUM["D_WGS_1984",
         SPHEROID["WGS_1984",6378137.0,298.257223563]],
       PARAMETER["Vertical_Shift",0.0],
       PARAMETER["Direction",1.0],UNIT["Meter",1.0]]

Could you also try the following which is the variation of your intial proposal but with 2000 and Sphere switched, consistently with what is found in the ESRI database ?

GEOGCS["GCS_Mars_2000_Sphere",
       DATUM["D_Mars_2000_Sphere",
         SPHEROID["Mars_2000_Sphere",3396190.0,0.0]],
       PRIMEM["Reference Meridian",0.0],
       UNIT["Degree",0.0174532925199433]],
VERTCS["Mars_2000_Sphere",
       DATUM["D_Mars_2000_Sphere",
         SPHEROID["Mars_2000_Sphere",3396190.0,0.0]],
       PARAMETER["Vertical_Shift",0.0],
       PARAMETER["Direction",1.0],UNIT["Meter",1.0]]
dpmayerUSGS commented 3 years ago

@rouault I know for fact the WGS84 example works without complaint in Arc. I will have to get back to you on the Mars example with 2000 and Sphere switched.

Just to clarify my earlier comment: Arc is loading the points from shapefiles using those other 3 sets of WKT and seems to place them in the correct location, but it complains that it doesn't recognize the CRS. This is consistent with ESRI's long history of tolerating "not quite perfect" CRSes. It seems to recognize that it's a geographic CRS, with coordinates expressed as latitude and longitude in degrees, but just didn't recognize the various strings in the WKT.

rouault commented 3 years ago

Fix in https://github.com/OSGeo/PROJ/pull/2902 such that the following WKT can be round-trip as WKT1:ESRI by PROJ

GEOGCS["Mars2000(Sphere)",DATUM["Mars2000(Sphere)",SPHEROID["Mars2000(Sphere)",3396190.0,0.0]],PRIMEM["Reference_Meridian",0.0],UNIT["Degree",0.0174532925199433]],VERTCS["Mars2000(Sphere)",DATUM["Mars2000(Sphere)",SPHEROID["Mars2000(Sphere)",3396190.0,0.0]],PARAMETER["Vertical_Shift",0.0],PARAMETER["Direction",1.0],UNIT["Meter",1.0]]