Database: register NZGD2000 -> ITRF96 transformation for NZGD2000 database

[rouault]

Related to OSGeo/PROJ-data#22

An entry is added in the other_transformation table, using a raw PROJ
string. Later, once the deformation model has been registered in EPSG, we'll
have to add code to map the EPSG transformation method and parameters to
+proj=defmodel

In the meantime, this works pretty well:

$ src/projinfo -s EPSG:4959 -t EPSG:7907
Candidate operations found: 1
-------------------------------------
Operation No. 1:

PROJ:NZGD2000-20180701, NZGD2000 to ITRF96 deformation model, unknown accuracy, New Zealand

PROJ string:
+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=rad +xy_out=deg

WKT2:2019 string:
COORDINATEOPERATION["NZGD2000 to ITRF96 deformation model",
    VERSION["20180701"],
    SOURCECRS[
        GEOGCRS["NZGD2000",
            DATUM["New Zealand Geodetic Datum 2000",
                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",4959]]],
    TARGETCRS[
        GEOGCRS["ITRF96",
            DATUM["International Terrestrial Reference Frame 1996",
                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",7907]]],
    METHOD["PROJ-based operation method: +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=rad +xy_out=deg"],
    USAGE[
        SCOPE["unknown"],
        AREA["New Zealand"],
        BBOX[-55.95,160.6,-25.88,-171.2]],
    ID["PROJ","NZGD2000-20180701"],
    REMARK["New Zealand Deformation Model. Defines the secular model (National Deformation Model) and patches for significant deformation events since 2000"]]

$ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cct -d 8 +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=axisswap +order=2,1 +step +proj=defmodel +model=nz_linz_nzgd2000-20180701.json +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=rad +xy_out=deg
 -40.99999402   172.99999938    0.00130333     2016.5000

$ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cs2cs -f "%.8f" EPSG:4959 EPSG:7907
-40.99999402	172.99999938 0.00130333 2016.5

[rouault]

CC @ccrook

[rouault]

Works well also using the CRS names instead of their codes

echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cs2cs -f "%.8f" NZGD2000 ITRF96
-40.99999402	172.99999938 0.00130333 2016.5

[ccrook]

@rouault Thanks - that is very cool. The final step in most NZ transformations will be to convert to/from current ITRF realisations, which will need to use the ITRF96-ITRFxxxx transformations used in NZ (which for historical reasons are the same as the NGA transformations but different from IERS). These are defined explicitly in EPSG, eg coordinate transformation EPSG:9082. Is there a way of encouraging PROJ to use them.

[rouault]

Is there a way of encouraging PROJ to use them.

Normally, PROJ has logic to infer the concatenations of existing steps, but this didn't work here for a reason I wasn't courageous enough to explore. And if it had worked automatically, there would have been a risk for it to use the worldwide ITRF96 <--> ITRFxx transformations rather than the NZ ones. So I've created manually records in the concatenated_operation table to be explicit. When submitting to EPSG, you'll probably want to submit them too.

So now the following works:

$ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cs2cs -f "%.8f" NZGD2000 ITRF96
-40.99999402	172.99999938 0.00130333 2016.5
$ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cs2cs -f "%.8f" NZGD2000 ITRF97
-40.99999439	172.99999925 0.07020264 2016.5
$ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cs2cs -f "%.8f" NZGD2000 ITRF2000
-40.99999404	172.99999908 0.03457527 2016.5
$ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cs2cs -f "%.8f" NZGD2000 ITRF2005
-40.99999382	172.99999909 -0.00214465 2016.5
$ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cs2cs -f "%.8f" NZGD2000 ITRF2008
-40.99999377	172.99999908 -0.00892236 2016.5
$ echo "-41 173 0 2016.5" | PROJ_NETWORK=ON src/cs2cs -f "%.8f" NZGD2000 ITRF2014
-40.99999377	172.99999911 -0.00786504 2016.5

[rouault]

@ccrook Do the definitions & results of the concatenated operations look good to you ?

[rouault]

@ccrook Gentle ping w.r.t above question

[ccrook]

Hi Even. This looks great - really cool. I have checked all the test results and they exactly match the LINZ coordinate conversion software.

I am trying to add this to a .gie for the NZGD2000 ITRF conversions, I tried

use_proj4_init_rules true
operation proj=pipeline step init=epsg:4959 inv step init=epsg:7911

which was read, but it failed to convert coordinates, ie

 PROJ_NETWORK=ON gie nzgd2000-20180701-itrf.gie 
-------------------------------------------------------------------------------
Reading file 'nzgd2000-20180701-itrf.gie'
     -----
     FAILURE in nzgd2000-20180701-itrf.gie(14):
     expected: 172.999999083 -40.999993769 -0.0089
     got:      173.000000000000   -41.000000000000
     deviation:  696.321525 mm,  expected:  1.000000 mm
-------------------------------------------------------------------------------
total:  1 tests succeeded,  0 tests skipped,  1 tests FAILED!
-------------------------------------------------------------------------------

But I haven't been able to make this work yet. Can you see what is wrong with this?

I hadn't seen the pipline step init inv step init before until I started grepping the gie files in the PROJ test files - very straightforward. Can I use the codes NZGD2000 etc in gie using a format like this. It didn't seem to work when I just replaced eg epsg:4959 with NZGD2000.

nzgd2000-20180701-itrf.gie.txt

[rouault]

operation proj=pipeline step init=epsg:4959 inv step init=epsg:7911

This won't work for your purposes. This uses the legacy init=epsg:XXXX syntax that expands CRS to proj.4 strings. This can only work for simple things, like projected CRS that have one single towgs84 transformation. Not for NZGD2000 transformation. We would perhaps need something emulating cs2cs syntax, but that doesn't exist yet in gie.

Your best option currently would probably be to add a new test in test/cli/testvarious , which uses cs2cs, and update test/cli/tv_out.dist. But just thinking this would require the grids to be available, so would be more complicated than that, since we don't want them to depend on the network, so you'd have to add those which would be needed in data/tests/