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

concatenated operations #3819

Closed
barry-gallagher opened this issue Jul 11, 2023 · 4 comments
Closed

concatenated operations #3819

barry-gallagher opened this issue Jul 11, 2023 · 4 comments

Comments

@barry-gallagher
Copy link

Hello,

We are making a lot of vertical datums, crs and transforms which seem to work thanks to #3616 (and #3617). We wanted to add some concatenated operations and those seem to work unless I add a geographic 3d to geographic 3d step at the end. Is the problem in trying to add a 3d transform at the end of what was mostly vertical, related to our custom authority or missing a step or method code? Thanks in advance for any help or pointers on what I have done wrong, I've attached my prototype proj.db.zip.

For example, we have tidal datums of Mean Lower Low Water (MLLW) and Local Mean Sea Level (LMSL) and a transform between them. We also have a transform from LMSL to Geoid09. I've tried to create two concatenated operations, one to go from MLLW to EPSG:4893 NAD83(NSRS2007) and the other to go from MLLW to EPSG:6319 NAD83(2011). The first one works but the second fails.

Both transforms start with the same three steps.

  1. MLLW (NOAA:752) -> LMSL (NOAA:348) using grid_transformation NOAA:754
  2. LMSL (NOAA:348) -> Geoid09(Tried both EPSG:5500 and EPSG:5703) using grid_transformation NOAA:350
  3. Geoid09 -> NAD83(NSRS2007) (EPSG:4893) using grid_transform EPSG:9173

The second transform adds one more step.
4) NAD83(NSRS2007) (EPSG:4893) to NAD83(2011) (EPSG:6319) using grid_transform EPSG:8559

Running from our tidal EPSG:4759+NOAA:752 NAD83(NSRS2007) 2D+MLLW to EPSG:4893 NAD83(NSRS2007) works
Trying to extend to NAD83(2011) fails.

projinfo
Rel. 9.2.1, June 1st, 2023


cs2cs "EPSG:4759+NOAA:752" "EPSG:4893" < coords.txt
26d36'N 97dW -22.923


cs2cs "EPSG:4759+NOAA:752" "EPSG:6319" < coords.txt
Rel. 9.2.1, June 1st, 2023
<cs2cs>:
cannot initialize transformation
cause: (null)
program abnormally terminated

Projinfo doesn't return anything for for my attempt above but by leaving off the horizontal component it does find the concatenated_operation entry. It gives the message "Error when exporting to PROJ string: Invalid nature of source and/or targetCRS for Geographic/Geocentric conversion" which I don't understand or find reference to.

projinfo -s EPSG:4759+NOAA:752 -t EPSG:6319
Candidate operations found: 0

projinfo -s NOAA:752 -t EPSG:6319
Candidate operations found: 1
-------------------------------------
Operation No. 1:

NOAA:758, TXshlmat01_8301_MLLW_1983-2001_NAD83(2011), 0.0 m, Texas - Gulf of Mexico, offshore north Port Aransas, south to US/Mexico border

PROJ string:
Error when exporting to PROJ string: Invalid nature of source and/or targetCRS for Geographic/Geocentric conversion

WKT2:2019 string:
CONCATENATEDOPERATION["TXshlmat01_8301_MLLW_1983-2001_NAD83(2011)",
    SOURCECRS[
        VERTCRS["TXshlmat01_8301_MLLW_1983-2001",

Running projinfo on the MLLW to NSRS2007 seems to work whether I supply the horizontal or not.

projinfo -s NOAA:752 -t EPSG:4893
Candidate operations found: 1
-------------------------------------
Operation No. 1:

NOAA:756, TXshlmat01_8301_MLLW_1983-2001_NAD83(NSRS2007), 0.0 m, Texas - Gulf of Mexico, offshore north Port Aransas, south to US/Mexico border

PROJ string:
+proj=pipeline
  +step +proj=vgridshift +grids=TXshlmat01_8301/mllw.gtx +multiplier=1
  +step +proj=vgridshift +grids=TXshlmat01_8301/tss.gtx +multiplier=1
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_noaa_geoid09_conus.tif +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

WKT2:2019 string:
CONCATENATEDOPERATION["TXshlmat01_8301_MLLW_1983-2001_NAD83(NSRS2007)",
    SOURCECRS[
        VERTCRS["TXshlmat01_8301_MLLW_1983-2001",
            VDATUM["MLLW_1983-2001"],
            CS[vertical,1],
                AXIS["gravity-related height (H)",up,
                    LENGTHUNIT["metre",1]],
            ID["NOAA",752]]],
    TARGETCRS[
        GEOGCRS["NAD83(NSRS2007)",
            DATUM["NAD83 (National Spatial Reference System 2007)",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,3],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["ellipsoidal height (h)",up,
                    ORDER[3],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",4893]]],
    STEP[
        COORDINATEOPERATION["TXshlmat01_8301_MLLW_1983-2001_LMSL_1983-2001",
            SOURCECRS[
                VERTCRS["TXshlmat01_8301_MLLW_1983-2001",
                    VDATUM["MLLW_1983-2001"],
                    CS[vertical,1],
                        AXIS["gravity-related height (H)",up,
                            LENGTHUNIT["metre",1]],
                    ID["NOAA",752]]],
            TARGETCRS[
                VERTCRS["TXshlmat01_8301_LMSL_1983-2001",
                    VDATUM["LMSL_1983-2001"],
                    CS[vertical,1],
                        AXIS["gravity-related height (H)",up,
                            LENGTHUNIT["metre",1]],
                    ID["NOAA",348]]],
            METHOD["Vertical Offset by Grid Interpolation (gtx)",
                ID["EPSG",1084]],
            PARAMETERFILE["Vertical offset file","TXshlmat01_8301/mllw.gtx"],
            OPERATIONACCURACY[0.0],
            ID["NOAA",754],
            REMARK["Texas - Gulf of Mexico, offshore north Port Aransas, south to US/Mexico border to Mean Lower Low Water 1983-2001 National Tidal Datum Epoch to Local Mean Sea Level 1983-2001 National Tidal Datum Epoch"]]],
    STEP[
        COORDINATEOPERATION["TXshlmat01_8301_LMSL_1983-2001_NAVD88_GEOID09_CONUS",
            SOURCECRS[
                VERTCRS["TXshlmat01_8301_LMSL_1983-2001",
                    VDATUM["LMSL_1983-2001"],
                    CS[vertical,1],
                        AXIS["gravity-related height (H)",up,
                            LENGTHUNIT["metre",1]],
                    ID["NOAA",348]]],
            TARGETCRS[
                COMPOUNDCRS["NAD83(NSRS2007) + NAVD88 height",
                    GEOGCRS["NAD83(NSRS2007)",
                        DATUM["NAD83 (National Spatial Reference System 2007)",
                            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]]],
                    VERTCRS["NAVD88 height",
                        VDATUM["North American Vertical Datum 1988"],
                        CS[vertical,1],
                            AXIS["gravity-related height (H)",up,
                                LENGTHUNIT["metre",1]]],
                    ID["EPSG",5500]]],
            METHOD["Vertical Offset by Grid Interpolation (gtx)",
                ID["EPSG",1084]],
            PARAMETERFILE["Vertical offset file","TXshlmat01_8301/tss.gtx"],
            OPERATIONACCURACY[0.0],
            ID["NOAA",350],
            REMARK["Texas - Gulf of Mexico, offshore north Port Aransas, south to US/Mexico border to Local Mean Sea Level 1983-2001 National Tidal Datum Epoch to NAVD88_GEOID09_CONUS"]]],
    STEP[
        COORDINATEOPERATION["Inverse of NAD83(NSRS2007) to NAVD88 height (1)",
            SOURCECRS[
                VERTCRS["NAVD88 height",
                    VDATUM["North American Vertical Datum 1988"],
                    CS[vertical,1],
                        AXIS["gravity-related height (H)",up,
                            LENGTHUNIT["metre",1]],
                    ID["EPSG",5703]]],
            TARGETCRS[
                GEOGCRS["NAD83(NSRS2007)",
                    DATUM["NAD83 (National Spatial Reference System 2007)",
                        ELLIPSOID["GRS 1980",6378137,298.257222101,
                            LENGTHUNIT["metre",1]]],
                    PRIMEM["Greenwich",0,
                        ANGLEUNIT["degree",0.0174532925199433]],
                    CS[ellipsoidal,3],
                        AXIS["geodetic latitude (Lat)",north,
                            ORDER[1],
                            ANGLEUNIT["degree",0.0174532925199433]],
                        AXIS["geodetic longitude (Lon)",east,
                            ORDER[2],
                            ANGLEUNIT["degree",0.0174532925199433]],
                        AXIS["ellipsoidal height (h)",up,
                            ORDER[3],
                            LENGTHUNIT["metre",1]],
                    ID["EPSG",4893]]],
            METHOD["Inverse of Geographic3D to GravityRelatedHeight (gtx)",
                ID["INVERSE(EPSG)",9665]],
            PARAMETERFILE["Geoid (height correction) model file","us_noaa_geoid09_conus.tif"],
            OPERATIONACCURACY[0.05],
            ID["INVERSE(DERIVED_FROM(EPSG))",9173],
            REMARK["Uses Geoid09 hybrid model. Replaced by 2012 model (CT code 6326)."]]],
    OPERATIONACCURACY[0.0],
    USAGE[
        SCOPE["Coastal hydrography."],
        AREA["Texas - Gulf of Mexico, offshore north Port Aransas, south to US/Mexico border"],
        BBOX[-90,-180,90,180]],
    ID["NOAA",756],
    REMARK["TXshlmat01_8301_MLLW_1983-2001 to NAD83(NSRS2007)"]]


projinfo -s EPSG:4759+NOAA:752 -t EPSG:4893
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, TXshlmat01_8301_MLLW_1983-2001_NAD83(NSRS2007), 0.0 m, Texas - Gulf of Mexico, offshore north Port Aransas, south to US/Mexico border

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=TXshlmat01_8301/mllw.gtx +multiplier=1
  +step +proj=vgridshift +grids=TXshlmat01_8301/tss.gtx +multiplier=1
  +step +proj=vgridshift +grids=us_noaa_geoid09_conus.tif +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

WKT2:2019 string:
COORDINATEOPERATION["TXshlmat01_8301_MLLW_1983-2001_NAD83(NSRS2007)",
    SOURCECRS[
        COMPOUNDCRS["NAD83(NSRS2007) + TXshlmat01_8301_MLLW_1983-2001",
            GEOGCRS["NAD83(NSRS2007)",
                DATUM["NAD83 (National Spatial Reference System 2007)",
                    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",4759]],
            VERTCRS["TXshlmat01_8301_MLLW_1983-2001",
                VDATUM["MLLW_1983-2001"],
                CS[vertical,1],
                    AXIS["gravity-related height (H)",up,
                        LENGTHUNIT["metre",1]],
                ID["NOAA",752]]]],
    TARGETCRS[
        GEOGCRS["NAD83(NSRS2007)",
            DATUM["NAD83 (National Spatial Reference System 2007)",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,3],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["ellipsoidal height (h)",up,
                    ORDER[3],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",4893]]],
    METHOD["PROJ-based operation method (approximate): +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=TXshlmat01_8301/mllw.gtx +multiplier=1 +step +proj=vgridshift +grids=TXshlmat01_8301/tss.gtx +multiplier=1 +step +proj=vgridshift +grids=us_noaa_geoid09_conus.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1"],
    OPERATIONACCURACY[0.0],
    USAGE[
        SCOPE["unknown"],
        AREA["Texas - Gulf of Mexico, offshore north Port Aransas, south to US/Mexico border"],
        BBOX[-90,-180,90,180]],
    REMARK["TXshlmat01_8301_MLLW_1983-2001 to NAD83(NSRS2007)"]]

Thanks again,
Barry
rouault added a commit to rouault/PROJ that referenced this issue Jul 11, 2023
@rouault
fixStepsDirection(): enable to define a concatenated operationg endin…
9a99520 
…g with a NADCOND (3D) transformation (use case of refs OSGeo#3819)
@rouault
Copy link
Member

rouault commented Jul 11, 2023

I've submitted #3820 with a number of related fixes, but that doesn't solve your end-to-end issue.
#3820 now enables to instanciate NOAA:758 by automatically promoting the EPSG:8559 transformation to refer to 3D CRS instead of the 2D ones of EPSG.
When #3820 is applied and trying cs2cs "EPSG:4759+NOAA:752" "EPSG:6319", there are a number of inconsistent concatenated operations (missing steps) that make PROJ fail. I've added a PROJ_IGNORE_INSTANCIATION_ERRORS env variable that you can set to YES to skip those broken records, but ultimately the database should be fixed.

Basically you're pushing PROJ beyond its limits with the use of complex concatenated operations between a VertCRS and a GeogCRS, and then hoping that PROJ will apply that to a CompoundCRS
I'm not clear which pipeline PROJ should apply for NAD83(NSRS2007) + TXshlmat01_8301_MLLW_1983-2001 to NAD83(2011).
I haven't debugged that deeply, but I believe that when (if ?) PROJ finds the NOAA:758 operation (..), which is a transformation between TXshlmat01_8301_MLLW_1983-2001 and NAD83(2011), it is confused because NAD83(NSRS2007) != NAD83(2011). Currently it thus tries to do in an initial step NAD83(2011) -> NAD83(NSRS2007) and then applies the NOAA:756 operation, which is also fishy because the target CRS of that operation is NAD83(NSRS2007) and not NAD83(2011)
Your best hope is probably to define completely the pipeline from EPSG:4759+NOAA:752 to EPSG:6319 as a PROJ pipeline.

rouault added a commit to rouault/PROJ that referenced this issue Jul 11, 2023
@rouault
fixStepsDirection(): enable to define a concatenated operationg endin…
70f1393 
…g with a NADCOND (3D) transformation (use case of refs OSGeo#3819)
@barry-gallagher
Copy link
Author

Thanks for your quick response. Are you suggesting I register a custom compound crs for EPSG:4759+NOAA:752 and then make an other_transformation with the steps in a pipeline? If I add a compound_crs (say NOAA:3000) to the proj.db would PROJ realize that EPSG:4759+NOAA:752 was equivalent and use the other_transformation for NOAA:3000?

Maybe I should ask for more general advice. We have these vertical datums, which depending on when they were made, relate back to different geographic 3d crs NAD83 (CORS96, NSRS2007, 2011 etc) or IGS/ITRF (08, 14 etc). I was trying to make it easier for a user to get from a vertical datum and its corresponding horizontal CRS to a different vertical datum in another horizontal crs. Maybe someone needs the MLLW from a NAD83(NSRS2007) model to be transformed to MHHW in ITRF14. From your answer #3616 before about PROJ only making transform that makes two jumps (like a->b->c but not a->b->c->d) then I thought I would supply concatenated operations so everything could meet at NAD83(2011) and there would always be a two step a->NAD83(2011)->c path to be found.

Would you still suggest making pipelines (other_transformations), some other strategy or should I just stop at the 3d CRS applicable to a specific datum and tell people to learn to use PROJ?

rouault added a commit to rouault/PROJ that referenced this issue Jul 15, 2023
@rouault
fixStepsDirection(): enable to define a concatenated operationg endin…
f1d29da 
…g with a NADCOND (3D) transformation (use case of refs OSGeo#3819)
rouault added a commit that referenced this issue Jul 17, 2023
@rouault
Merge pull request #3820 from rouault/fix_3819
c5b74f1 
Various fixes related to #3819
@rouault
Copy link
Member

rouault commented Jul 17, 2023

Are you suggesting I register a custom compound crs for EPSG:4759+NOAA:752 and then make an other_transformation with the steps in a pipeline?

yes

If I add a compound_crs (say NOAA:3000) to the proj.db would PROJ realize that EPSG:4759+NOAA:752 was equivalent and use the other_transformation for NOAA:3000?

If the name of NOAA:3000 is exactly {name of EPSG:759} + {name of NOAA:752}, then yes I believe it should manage to identify it. If not, that could be relatively easily fixed I believe

then I thought I would supply concatenated operations so everything could meet at NAD83(2011) and there would always be a two step a->NAD83(2011)->c path to be found.

That makes sense. It is just that PROJ is confused when having just the vertical transformation "TXshlmat01_8301_MLLW_1983-2001 to NAD83(2011)" to infer how to do the more complex "NAD83(NSRS2007) + TXshlmat01_8301_MLLW_1983-2001 to NAD83(2011)"

Would you still suggest making pipelines (other_transformations), some other strategy or should I just stop at the 3d CRS applicable to a specific datum

I'm not sure which answer to give that. The more registered transformations you add, the easier & more reliable PROJ will behave, but the more complicated maintenance of the registry will be. Having a pivotal CRS like NAD83(2011) as you suggest can makes sense I guess

@barry-gallagher
Copy link
Author

Thanks very much. I will try making the compound crs entries with naming as you suggest.

@barry-gallagher barry-gallagher mentioned this issue Jul 24, 2023
Closed
@barry-gallagher barry-gallagher mentioned this issue Apr 17, 2024
Open
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
2 participants
@rouault @barry-gallagher