Change order of horizontal and vertical operations when dealing with WKT1 / PROJ4 compound CRS #3127
Conversation
|
I'm not sure if I like this suggested change. A switch back to PROJ4 behaviour for PROJ4-style strings is less confusing for some users, but this is mainly an advantage for users that still use PROJ4-based code. For users of PROJ4-style strings in new PROJ versions it has a major disadvantage! Countries using a CRS with a horizontal grid will need to keep publishing two geoid grids forever: one geoid grid for an international CRS (e.g. ETRS89, ITRF) and one geoid grid with an interpolation CRS for the grid-shifted CRS. This is a burden for the authorities. For users this will be a source of confusion on what geoid grid to use in what situation (provided that the authorities will provide both geoid grids in the first place). See examples below for the transformation between ETRS89 (assuming a null-transformation from WGS84 for simplicity) and the national horizontal (RD) and vertical (NAP) CRS of the Netherlands. With the suggested change:
Without the suggested change:
(If you really want this change, I would strongly recommend to simultaneously fix the specification of an interpolation CRS in a Geodetic TIFF Grid, as we discussed here: https://lists.osgeo.org/pipermail/proj/2022-February/010532.html) |
What was the interpolation CRS of naptrans2018 : ETRS89 or Amersfoort (I see the CRS of nl_nsgi_naptrans2018.tif was ETRS89, but perhaps this was wrong) ? Was the PROJ.4 string and grids you indicate in 1) designed to work with PROJ >= 6, or earlier versions (if earlier versions, then it means the interpolation CRS of naptrans2018 was Amersfoort). From OSGeo/PROJ-data#80 where you removed naptrans2018 from PROJ-data, I understand this was for earlier versions, so this PR would make that use case work correctly again.
This PR doesn't change behavior for this transform (you need to add --3d to projinfo and cs2cs due to #3119 that was merged in 9.1dev). You'll get the same pipeline: $ bin/projinfo -s "+proj=lonlat +towgs84=0,0,0 +ellps=GRS80 +type=crs" -t "+proj=lonlat +towgs84=0,0,0 +ellps=GRS80 +geoidgrids=nlgeo2018.gtx +type=crs" -o PROJ --3d -q
That wouldn't help here given that currently, PROJ defers the opening of grids until it needs to transform a point (mostly for performance reasons when a lot of pipelines/grids are considered, like the tens of US state grids for NAD83 to WGS84), whereas the pipeline chaining the various steps (+proj=hgridshift / +proj=vgridshift) is established much before. Well, maybe something special could be done for grids that are specific in the CRS definition, but that would work only for GeoTIFF grids, not .gtx/.gsb that lack metadata. |
|
Too bad we couldn't get this and the related changes in |
|
The naptrans2018 grid file has Amersfoort as interpolation CRS. The PROJ4 strings (1) and (2) that I mentioned are designed to work with PROJ >=9.1 if this PR would be implemented (as well as with PROJ <6). I'll edit my previous message to make this clearer. The geoid grid file naptrans2018 needs a confusing interpolation CRS. Only to make it slightly less confusing, I recommend to fix the specification of an interpolation CRS in a Geodetic TIFF Grid (https://lists.osgeo.org/pipermail/proj/2022-February/010532.html). This indeed doesn't fix the behaviour of PROJ, but would only help to make human reading of the metadata easier (and to make conventional GeoTIFF software to use the interpolation CRS as horizontal CRS and thus showing the grid as an image on the correct location). My main objection against this PR is that we would need two different geoid grid files for the two use cases (1) and (2). While without this PR we need just one geoid grid file for both use cases (1) and (2) in PROJ >=6. |
Intuitively this is what I would expect to happen for the described transformation based on the PROJ.4-style input. I am not super happy about changing behaviour again thouvh. It will make it more or less impossible to predict the outcome of a transformation. Of course it can be argued that you shouldn't use PROJ.4-style strings to initialize transformations and as such you'll just have to live with the software being a bit flaky when using legacy features. At this point I don't have a strong preference for one over the other.
Perhaps the grid CRS could be tracked in the database alongside other grid information? That should allow PROJ to do the right thing without opening the grid first. |
That wouldn't be completely unprecedented. I know the Canadians have several variants of geoid models that differ by the interpolation CRS: cf https://github.com/OSGeo/PROJ-data/tree/master/ca_nrc#canada-cgvd2013cgg2013a-height----nad83csrs vs https://github.com/OSGeo/PROJ-data/tree/master/ca_nrc#canada-cgvd2013cgg2013a-height----itrf2008 |
To me that would only make sense if the input would be in a different CRS. It is weird that when transforming 3D from ETRS89 to RDNAP one would need a different geoid than when transforming ETRS89 to NAP height only. In the past we did have a geoid model using Amerfoort, but that used 3D Amersfoort, thus as interpolation CRS as well as as source CRS for the height. The Amersfoort CRS is now not being used anymore, except as an intermediate result before projecting to RD: "Use of geographic CRS Amersfoort (code 4289) for spatial referencing is discouraged".[1] The naptrans2018 grid would therefore only be needed for PROJ4 strings and no other software would use it. |
That could be figured out from what is available in the database : use the grid_alternatives table to find the official EPSG name of the grid, and then browse the grid_transformation table to find a transformation that use that grid, and cry if there are several records that use the same grid with different CRS :-) which I'm pretty sure happens for historical CRS / low accuracy grids that can be indifferently used against WGS84/ETRS89/NAD83. Anyway, given the lack of consensus on this, I guess I'll close this PR and #3126. |
Switzerland is another example of a country with two grid files for the same geoid model, depending on the interpolation CRS: The datum of CH1903+ tries to "fit the sea level" in Switzerland. That produces geoid undulations on the geoid file that uses LV95 on the order of centimeters or few meters, compared to the 50 m in the geoid file that uses ETRS89. So far those transformations are not in EPSG, but added "manually" in PROJ. |
|
What is the plan here? |
I'm not sure... No way to please everyone here |
Nice idea, but a more descriptive name of the flag would be clearer, something like: |
|
@Jochem-L I like the idea of a flag like that, explicitly choosing the behavior in either way. Not really how proj handled things in the past, but it really makes things explicit that weren't well defined before. |
|
ok, I've added a +geoid_crs=WGS84/horizontal_crs which behaves as specified in the below extract of the updated transformation.rst: The ``+geoid_crs=WGS84/horizontal_crs`` parameter was added in PROJ 9.1 to solve
an ambiguity that has been rampant for a long time. Before PROJ 6, the interpolation
CRS of the geoid transformation was assumed to be the one of the horizontal CRS.
Starting with PROJ 6, this was inadvertently changed to be WGS84.
When transforming from a CRS that specifies a horizontal datum shift
with ``+towgs84`` or ``+nadgrids``, and a vertical one with ``+geoidgrids`` with
``+geoid_crs=horizontal_crs``, to WGS84, the vertical transformation is done before
the horizontal one.
::
$ projinfo -s "+proj=longlat +ellps=GRS80 +nadgrids=@foo +geoidgrids=@bar +geoid_crs=horizontal_crs +type=crs" -t EPSG:4979 -o PROJ -q
+proj=pipeline
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +proj=vgridshift +grids=@bar +multiplier=1
+step +proj=hgridshift +grids=@foo
+step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m
+step +proj=axisswap +order=2,1
Conversely, when ``+geoid_crs=WGS84`` (or omitting it with PROJ >= 6),
the horizontal transformation is done before the vertical one.
::
$ projinfo -s "+proj=longlat +ellps=GRS80 +nadgrids=@foo +geoidgrids=@bar +geoid_crs=WGS84 +type=crs" -t EPSG:4979 -o PROJ -q
+proj=pipeline
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +proj=hgridshift +grids=@foo
+step +proj=vgridshift +grids=@bar +multiplier=1
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1
|
|
That's an elegant solution, although I don't really like the use of the name WGS84 in |
…ION PROJ4_GRIDS
Previously, when importing a WKT1 like
COMPD_CS["NAD83 + NAVD88 height",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4269"]],
VERT_CS["NAVD88 height",
VERT_DATUM["North American Vertical Datum 1988",2005,
EXTENSION["PROJ4_GRIDS","@foo.gtx"],
AUTHORITY["EPSG","5103"]],
UNIT["metre",1.0,
AUTHORITY["EPSG","9001"]],
AXIS["Gravity-related height",UP],
AUTHORITY["EPSG","5703"]]]
we modeled the vertical part as a BoundCRS of a vertical CRS with a hub
CRS being EPSG:4979 and a transformation from the vertical CRS to
EPSG:4979 using the specified grid.
This is quite questionable, and doesn't follow the semantics of the
GDAL <= 2.4 / PROJ < 6 era, where the grid was referenced to the
horizontal datum, rather than WGS 84.
So change this to use the promoted-to-3D geographic CRS of the compound CRS
as the hub CRS / geographic 3D of the vertical transformation.
Said otherwise, the WKT2 export of the above WKT is now:
COMPOUNDCRS["NAD83 + NAVD88 height",
GEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",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]],
ID["EPSG",4269]],
BOUNDCRS[
SOURCECRS[
VERTCRS["NAVD88 height",
VDATUM["North American Vertical Datum 1988"],
CS[vertical,1],
AXIS["gravity-related height",up,
LENGTHUNIT["metre",1]],
ID["EPSG",5703]]],
TARGETCRS[
GEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]],
ID["EPSG",6269]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,3],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["ellipsoidal height (h)",up,
ORDER[3],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
REMARK["Promoted to 3D from EPSG:4269"]]],
ABRIDGEDTRANSFORMATION["NAVD88 height to NAD83 ellipsoidal height",
METHOD["GravityRelatedHeight to Geographic3D"],
PARAMETERFILE["Geoid (height correction) model file","@foo.gtx",
ID["EPSG",8666]]]]]
…oidgrids, and modify order of horizontal and vertical steps Previously when importing "+proj=something +datum/+towgs84/+nadgrids +geoidgrids", the geoid transformation was referenced to WGS 84, instead of the horizontal datum, which was a breakage with respect to PROJ < 6 behaviour. Modify also the createOperations() logic in such situation to apply first the vertical transformation and then the horizontal one when going from the compound CRS to the final or intermediate geographic 3D CRS.
…al_crs is specified
rouault
force-pushed
the
wkt1_and_proj4_geoidgrid_import
branch
from
fa23235
to
906fc57
Compare
|
@kbevers opinion on latest developments ( #3127 (comment) ) ? |
I think I'm on the same page as @Jochem-L here. It solves the problem in nice way. I am not particularly happy about the implicit extension of the PROJ.4-syntax but I guess it's hard to avoid. |
This PR builds on top of #3125 and #3126, and specifically adds the following commit:
Previously when importing "+proj=something +datum/+towgs84/+nadgrids +geoidgrids",
the geoid transformation was referenced to WGS 84, instead of the
horizontal datum, which was a breakage with respect to PROJ < 6
behaviour.
Modify also the createOperations() logic in such situation to apply
first the vertical transformation and then the horizontal one when going
from the compound CRS to the final or intermediate geographic 3D CRS.
With all those changes, PROJ should behave very close to the PROJ < 6 era when doing transformations between CRS created from WKT1 and PROJ.4 and that have both horizontal & vertical transformations. The main effect is that the vertical grid is assumed to be registered against the geographic datum of the horizontal CRS, and no longer WGS 84 as was hardcoded in PROJ >= 6 and < 9.1. So a transformation doing "+nadgrids=nad_src +geoidgrids=geoid_src" to "+nadgrids=naddst +geoidgrids=geoid_dst" will now do "+proj=vgridshift +grids=geoid_src" -> "+proj=hgridshift +grids=nad_src" -> "+proj=hgridshift +grids=nad_dst" -> "+proj=vgridshift +grids=geoid_dst" (with appropriate direction of the transformation omitted here...)
@kbevers / @hobu : not sure if the above makes sense to you... Beyond being compatible with the legacy behaviour, it seems also more logical, although that can be argued and people could certainly have use cases where the vertical grid is referenced to the datum after the horizontal transformation (but if using WKT2 with explicit hub CRS for the vertical transformation, that should be correctly taken into account). The main question is "it appropriate at that time to bring back PROJ < 6 behaviour ?"
CC @Jochem-L as it might impact Dutch transformations using rdtrans and naptrans (but I believe there are considered legacy compared to a more modern approach mixing Helmert based and grid based transformation), as can been seen in the change of the expected result of one of the unit test. We had a discussion on the mailing list some time ago regarding the order of operations