[PROJ] Obtaining possible transformation pipelines

[Even Rouault]

> > 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?

PROJ strings can accurately represent coordinate operations (even between CRS 
that can't be themselves losslessly represented as PROJ strings).

> 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?

Because you don't use 

proj_operation_factory_context_set_spatial_criterion(pjContext, 
operationContext,
PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION)

hint: if you had tried with projinfo, it would have advised you about that ;-)

The reason is that for EPSG:3111 to EPSG:7899, the area of use is "small", so
there are 3 possible datum shift operations whose area of use covers the area 
of use of the CRSs.
However, when using GDA94 -> GDA2020, the area of use of those CRS is much 
larger, and there's only the Helmert-based transformation that covers the 
whole area. If you specify PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION (or "--
spatial-test intersects" with projinfo), you ask PROJ to return coordinate 
operations whose area of use intersects at least partially with the one of the 
CRS, but do not necessarily cover them completely.

> 
> 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?

Directly from the EPSG dataset. Strangely enough, the accuracy advertized for 
those grids is 0.05m , whereas the one of the Helmert base transformation is 
0.01m. If you don't believe me, check EPSG:8048 and EPSG:8446 in
https://www.epsg-registry.org/


-- 
Spatialys - Geospatial professional services
http://www.spatialys.com