[PROJ] Obtaining possible transformation pipelines

[Nyall Dawson]

On Fri, 24 May 2019 at 17:44, Even Rouault <even.rouault at spatialys.com> wrote:
> >
> > - I'll use proj_coordoperation_get_accuracy as a way of showing users
> > the accuracy of the operations, and proj_is_deprecated to highlight
> > deprecated ones.
>
> You'll never get deprecated coordinate operations when using
> proj_create_operations(). They are automatically removed (deprecated == broken
> in EPSG terms. Which is different from superseded). You can only get
> deprecated objects by querying them explicitly by code.

Ok, makes sense. I'll remove that filter and checkbox for proj 6 builds.


> > - After users have selected their desired pipeline, I'll store the
> > WKT2 2018 representation of this. For future transformations between
> > that pair of coordinate systems (inside that project), I'll use
> > proj_create with the stored wkt2 string to create a transformation
> > using the desired pathway.
>
> Storing the PROJ string would also be reasonable (PROJ strings for coordinate
> operations are lossless for the purpose of transforming coordinates), but the
> WKT2 2018 representation will give you a more traceable way. That said, I will
> not guarantee at 100% that the WKT2 2018 representation of the whole pipeline
> can be roundtripped in all situations, so I'd suggest perhaps:
>
> 1) export the coord op as PROJ string
> 2) export the coord op as WKT2_2018 and instanciate it as a PJ with
> proj_create()
> 3) export the PJ of 2) as PROJ string, and compare (string equality) with 1).
> If they match you can use the WKT2_2018 just fine. Otherwise use the PROJ
> string of 1) (and report such cases so that we can check if that can be
> improved)

Reading this sounds like I'm best off just to use the proj strings and
be done with it. But I was under the impression that not all CRSes
could be represented as a proj string? (e.g. EPSG:2218). Is it still
possible to represent operations from crses like EPS:2218 to proj
strings with loss?



> > - From my tests it looks fine, but I'd love double-confirmation that
> > this approach will correctly handle transformations which require
> > pivots, such as Vicgrid GDA94 -> Vicgrid GDA2020.

I'm seeing something I can't explain here.

Using this code to list operations from EPSG:3111 to EPSG:7899, I see
(as expected) 3 results, all using GDA94->GDA2020 as a pivot:

  PJ *crs1 = proj_create_from_database( pjContext, "EPSG", "3111",
PJ_CATEGORY_CRS, false, nullptr );
  PJ *crs2 = proj_create_from_database( pjContext, "EPSG", "7899",
PJ_CATEGORY_CRS, false, nullptr );
  PJ_OPERATION_FACTORY_CONTEXT *operationContext =
proj_create_operation_factory_context( pjContext, nullptr );
  proj_operation_factory_context_set_grid_availability_use( pjContext,
operationContext, PROJ_GRID_AVAILABILITY_IGNORED );
  PJ_OBJ_LIST *ops = proj_create_operations( pjContext, crs1, crs2,
operationContext );

Results in:

Inverse of Vicgrid + GDA94 to GDA2020 (1) + Vicgrid
+proj=pipeline +step +inv +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36
+lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +step +proj=push
+v_3 +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=0.06155
+y=-0.01087 +z=-0.04019 +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979
+s=-0.009994 +convention=coordinate_frame +step +inv +proj=cart
+ellps=GRS80 +step +proj=pop +v_3 +step +proj=lcc +lat_0=-37
+lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000
+ellps=GRS80

Inverse of Vicgrid + GDA94 to GDA2020 (2) + Vicgrid
+proj=pipeline +step +inv +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36
+lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +step
+proj=hgridshift +grids=GDA94_GDA2020_conformal_and_distortion.gsb
+step +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38
+x_0=2500000 +y_0=2500000 +ellps=GRS80

Inverse of Vicgrid + GDA94 to GDA2020 (3) + Vicgrid
+proj=pipeline +step +inv +proj=lcc +lat_0=-37 +lon_0=145 +lat_1=-36
+lat_2=-38 +x_0=2500000 +y_0=2500000 +ellps=GRS80 +step
+proj=hgridshift +grids=GDA94_GDA2020_conformal.gsb +step +proj=lcc
+lat_0=-37 +lon_0=145 +lat_1=-36 +lat_2=-38 +x_0=2500000 +y_0=2500000
+ellps=GRS80

BUT... listing operations directly from GDA94 (EPSG:4283) to GDA2020
(EPSG:7844) only gives one result:

GDA94 to GDA2020 (1)
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert
+xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart
+ellps=GRS80 +step +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019
+rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994
+convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step
+proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step
+proj=axisswap +order=2,1

Why aren't the options using using the
GDA94_GDA2020_conformal_and_distoration.gsb and
GDA94_GDA2020_conformal.gsb listed here?

And as a follow up question -- when checking the accuracy of the 3
pipelines given from EPSG:3111 to EPSG:7899,
proj_coordoperation_get_accuracy reports that the transforms using the
grid shift files only have an accuracy of 0.05m, vs 0.01m for the
non-grid pipeline (counter to my expectations). How is this accuracy
derived?

Nyall


>
> Sounds good.
>
> Even
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com