[PROJ] Obtaining possible transformation pipelines

classic Classic list List threaded Threaded
12 messages Options
Reply | Threaded
Open this post in threaded view
|

[PROJ] Obtaining possible transformation pipelines

Nyall Dawson
Hi list,

So, I think I'm finally onto the last piece of the QGIS proj 4 -> 6
puzzle. I think I've got the way forward here, but I'm after
confirmation that this is the right path to take.

The situation is:

- In current QGIS versions (proj 4 builds), a custom table of datum
transforms exists. This table maps src -> dest crs pairs via a
corresponding grid shift file. E.g. 4283 -> 7844 has entries for
GDA94_GDA2020_conformal_and_distortion.gsb and
GDA94_GDA2020_conformal.gsb.

- There's potentially multiple rows in this table per src/dest pair,
representing each of the different paths for transforming between
these systems

- Depending on users settings, either the "best" (an extremely loose
definition - basically the first non-deprecated path!) is selected
from this table when first transforming between a pair of systems, OR
users are shown a dialog listing all available pathways and they can
make a manual selection from these. (Settings are stored in projects,
with an option to save as the default choice on a particular install).

- This transformation is used for all future operations between those
two coordinate systems in that project

Now, obviously under proj 6 we want to dump all that, and use proj 6's
(far) superior intelligence and db to handle this. I would still like
to be able to (optionally) present users with a list of possible
transformation pathways which they can select from though.

Currently, I'm thinking of doing this by:

- Using proj_create_operations to obtain a list of the possible
transforms between a source/dest reference system

- Either using proj_get_name or proj_as_wkt as a way to represent
these transforms and present them to users (or is there a nicer
approach? The names I've encountered tend to be vague, e.g. "Inverse
of Vicgrid + GDA94 to GDA2020 (1) + Vicgrid"/"Inverse of Vicgrid +
GDA94 to GDA2020 (2) + Vicgrid", yet the wkt string is overly lengthy
and complex for most users).

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

- I'll use proj_operation_factory_context_set_grid_availability_use
with either PROJ_GRID_AVAILABILITY_IGNORED or
PROJ_GRID_AVAILABILITY_USED_FOR_SORTING so that users can see
transforms which they are missing the grids for, and highlight these
differently in the gui (and block their use)

- I'll confirm grid availability using
proj_coordoperation_is_instantiable, and
proj_coordoperation_get_grid_used in order to show the grid details
(and download url) to users

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

Questions:

- Does this approach seem reasonable? Is there anything I'm
misunderstanding, or better ways I could approach this?

- I never seem to obtain any deprecated results using
proj_create_operations -- are these filtered in the operation factory
context by default? If so, is there a way to enable them?

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

Many thanks!
Nyall
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Even Rouault-2
Nyall,

>
> - Either using proj_get_name or proj_as_wkt as a way to represent
> these transforms and present them to users (or is there a nicer
> approach? The names I've encountered tend to be vague, e.g. "Inverse
> of Vicgrid + GDA94 to GDA2020 (1) + Vicgrid"/"Inverse of Vicgrid +
> GDA94 to GDA2020 (2) + Vicgrid", yet the wkt string is overly lengthy
> and complex for most users).

No better suggestion.

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

>
> - I'll use proj_operation_factory_context_set_grid_availability_use
> with either PROJ_GRID_AVAILABILITY_IGNORED or
> PROJ_GRID_AVAILABILITY_USED_FOR_SORTING so that users can see
> transforms which they are missing the grids for, and highlight these
> differently in the gui (and block their use)

Sounds good

>
> - I'll confirm grid availability using
> proj_coordoperation_is_instantiable, and
> proj_coordoperation_get_grid_used in order to show the grid details
> (and download url) to users

Sounds good

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

> - I never seem to obtain any deprecated results using
> proj_create_operations -- are these filtered in the operation factory
> context by default? If so, is there a way to enable them?

See above. Removal of them is hardcoded.

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

Sounds good.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Nyall Dawson
On Fri, 24 May 2019 at 17:44, Even Rouault <[hidden email]> 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
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Even Rouault-2
> > 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
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Nyall Dawson
On Fri, 24 May 2019 at 20:48, Even Rouault <[hidden email]> wrote:

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

Thanks -- this was exactly what I was missing.

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

Odd. I suspect this was due to a misinterpretation of
http://www.icsm.gov.au/sites/default/files/DatumMattersT1FactSheet.pdf
. I'll find out what's required to get this fixed.

(hopefully) my last question:

If I have a proj string representing a coordinate operation, is there
any way to determine what grids are required for using it when those
grids are NOT available for use on the system?

Elsewhere I've used using proj_coordoperation_get_grid_used_count, but
this requires an existing coordinate operation object. Attempting to
create the (not available) coordinate operation using proj_create
fails with "Error -38: failed to load datum shift file". So as far as
I can see, there's no current way in API to determine the actual list
of missing grids which is preventing creation of a coordinate
operation. Or am I missing something?

Nyall
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Kristian Evers-2
Nyall,

I think you can use proj_coordoperation_get_grid_used() for that:

https://proj4.org/development/reference/functions.html#_CPPv433proj_coordoperation_get_grid_usedP10PJ_CONTEXTPK2PJiPPKcPPKcPPKcPPKcPiPiPi

/Kristian

-----Oprindelig meddelelse-----
Fra: PROJ <[hidden email]> På vegne af Nyall Dawson
Sendt: 27. maj 2019 03:38
Til: Even Rouault <[hidden email]>
Cc: PROJ <[hidden email]>
Emne: Re: [PROJ] Obtaining possible transformation pipelines

On Fri, 24 May 2019 at 20:48, Even Rouault <[hidden email]> wrote:

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

Thanks -- this was exactly what I was missing.

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

Odd. I suspect this was due to a misinterpretation of
http://www.icsm.gov.au/sites/default/files/DatumMattersT1FactSheet.pdf
. I'll find out what's required to get this fixed.

(hopefully) my last question:

If I have a proj string representing a coordinate operation, is there
any way to determine what grids are required for using it when those
grids are NOT available for use on the system?

Elsewhere I've used using proj_coordoperation_get_grid_used_count, but
this requires an existing coordinate operation object. Attempting to
create the (not available) coordinate operation using proj_create
fails with "Error -38: failed to load datum shift file". So as far as
I can see, there's no current way in API to determine the actual list
of missing grids which is preventing creation of a coordinate
operation. Or am I missing something?

Nyall
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Nyall Dawson
On Mon, 27 May 2019 at 17:42, Kristian Evers <[hidden email]> wrote:
>
> Nyall,
>
> I think you can use proj_coordoperation_get_grid_used() for that:
>
> https://proj4.org/development/reference/functions.html#_CPPv433proj_coordoperation_get_grid_usedP10PJ_CONTEXTPK2PJiPPKcPPKcPPKcPPKcPiPiPi

I can't use that one, because it requires an existing coordoperation
object, and I can't make the coordoperation when the grid doesn't
exist...

Nyall

>
> /Kristian
>
> -----Oprindelig meddelelse-----
> Fra: PROJ <[hidden email]> På vegne af Nyall Dawson
> Sendt: 27. maj 2019 03:38
> Til: Even Rouault <[hidden email]>
> Cc: PROJ <[hidden email]>
> Emne: Re: [PROJ] Obtaining possible transformation pipelines
>
> On Fri, 24 May 2019 at 20:48, Even Rouault <[hidden email]> wrote:
> >
> > > 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.
>
> Thanks -- this was exactly what I was missing.
>
> >
> > >
> > > 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/
>
> Odd. I suspect this was due to a misinterpretation of
> http://www.icsm.gov.au/sites/default/files/DatumMattersT1FactSheet.pdf
> . I'll find out what's required to get this fixed.
>
> (hopefully) my last question:
>
> If I have a proj string representing a coordinate operation, is there
> any way to determine what grids are required for using it when those
> grids are NOT available for use on the system?
>
> Elsewhere I've used using proj_coordoperation_get_grid_used_count, but
> this requires an existing coordinate operation object. Attempting to
> create the (not available) coordinate operation using proj_create
> fails with "Error -38: failed to load datum shift file". So as far as
> I can see, there's no current way in API to determine the actual list
> of missing grids which is preventing creation of a coordinate
> operation. Or am I missing something?
>
> Nyall
> _______________________________________________
> PROJ mailing list
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/proj
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Kristian Evers-2
I'm not sure this will work for your particular use case (since you are instantiating
a PJ from a proj-string), but there is proj_operation_factory_context_set_grid_availability_use
which allow you to ignore unavailable grids when creating PJ objects from using a
PJ_OPERATION_FACTORY_CONTEXT. Here's an example from the tests:

https://github.com/OSGeo/proj.4/blob/d67203a6f76a74f5ac029ff052dbcc72e3b59624/test/unit/test_c_api.cpp#L1316-L1351

which may guide you in the right direction.

/Kristian

-----Oprindelig meddelelse-----
Fra: Nyall Dawson <[hidden email]>
Sendt: 27. maj 2019 10:35
Til: Kristian Evers <[hidden email]>
Cc: Even Rouault <[hidden email]>; PROJ <[hidden email]>
Emne: Re: [PROJ] Obtaining possible transformation pipelines

On Mon, 27 May 2019 at 17:42, Kristian Evers <[hidden email]> wrote:
>
> Nyall,
>
> I think you can use proj_coordoperation_get_grid_used() for that:
>
> https://proj4.org/development/reference/functions.html#_CPPv433proj_coordoperation_get_grid_usedP10PJ_CONTEXTPK2PJiPPKcPPKcPPKcPPKcPiPiPi

I can't use that one, because it requires an existing coordoperation
object, and I can't make the coordoperation when the grid doesn't
exist...

Nyall

>
> /Kristian
>
> -----Oprindelig meddelelse-----
> Fra: PROJ <[hidden email]> På vegne af Nyall Dawson
> Sendt: 27. maj 2019 03:38
> Til: Even Rouault <[hidden email]>
> Cc: PROJ <[hidden email]>
> Emne: Re: [PROJ] Obtaining possible transformation pipelines
>
> On Fri, 24 May 2019 at 20:48, Even Rouault <[hidden email]> wrote:
> >
> > > 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.
>
> Thanks -- this was exactly what I was missing.
>
> >
> > >
> > > 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/
>
> Odd. I suspect this was due to a misinterpretation of
> http://www.icsm.gov.au/sites/default/files/DatumMattersT1FactSheet.pdf
> . I'll find out what's required to get this fixed.
>
> (hopefully) my last question:
>
> If I have a proj string representing a coordinate operation, is there
> any way to determine what grids are required for using it when those
> grids are NOT available for use on the system?
>
> Elsewhere I've used using proj_coordoperation_get_grid_used_count, but
> this requires an existing coordinate operation object. Attempting to
> create the (not available) coordinate operation using proj_create
> fails with "Error -38: failed to load datum shift file". So as far as
> I can see, there's no current way in API to determine the actual list
> of missing grids which is preventing creation of a coordinate
> operation. Or am I missing something?
>
> Nyall
> _______________________________________________
> PROJ mailing list
> [hidden email]
> https://lists.osgeo.org/mailman/listinfo/proj
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Even Rouault-2
In reply to this post by Nyall Dawson
> If I have a proj string representing a coordinate operation, is there
> any way to determine what grids are required for using it when those
> grids are NOT available for use on the system?

I don't think so. As you find out, PROJ will refuse to instanciate a
coordinate operation from a PROJ string that has a grids not present on the
system.

Asthere are not so many ways of expressing grids in PROJ grids, basically you
could look for +grids=xxxxx

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Nyall Dawson
On Tue, 28 May 2019 at 03:23, Even Rouault <[hidden email]> wrote:

>
> > If I have a proj string representing a coordinate operation, is there
> > any way to determine what grids are required for using it when those
> > grids are NOT available for use on the system?
>
> I don't think so. As you find out, PROJ will refuse to instanciate a
> coordinate operation from a PROJ string that has a grids not present on the
> system.
>
> Asthere are not so many ways of expressing grids in PROJ grids, basically you
> could look for +grids=xxxxx

That's what I was originally thinking could be a fallback option. But
then, I can't see any way to get the detailed metadata relating to a
grid except for proj_coordoperation_get_grid_used. So even if I
captured the grid name from the proj string, I still can't determine
the url and package name, etc relating to that grid. (or can I?)

What I'm aiming to do is take a proj coordinate operation string
stored in a project file, and show a nice warning to users if it can't
be used locally, along with the links to download the missing grid
files.

Nyall
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Even Rouault-2
On vendredi 31 mai 2019 14:05:38 CEST Nyall Dawson wrote:
> On Tue, 28 May 2019 at 03:23, Even Rouault <[hidden email]>
wrote:

> > > If I have a proj string representing a coordinate operation, is there
> > > any way to determine what grids are required for using it when those
> > > grids are NOT available for use on the system?
> >
> > I don't think so. As you find out, PROJ will refuse to instanciate a
> > coordinate operation from a PROJ string that has a grids not present on
> > the
> > system.
> >
> > Asthere are not so many ways of expressing grids in PROJ grids, basically
> > you could look for +grids=xxxxx
>
> That's what I was originally thinking could be a fallback option. But
> then, I can't see any way to get the detailed metadata relating to a
> grid except for proj_coordoperation_get_grid_used. So even if I
> captured the grid name from the proj string, I still can't determine
> the url and package name, etc relating to that grid. (or can I?)

In the C++ API, you have DatabaseContext::lookForGridInfo() that would cover
this. Would need to be exposed in the C API.

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Obtaining possible transformation pipelines

Nyall Dawson
On Fri, 31 May 2019 at 19:06, Even Rouault <[hidden email]> wrote:

>
> On vendredi 31 mai 2019 14:05:38 CEST Nyall Dawson wrote:
> > On Tue, 28 May 2019 at 03:23, Even Rouault <[hidden email]>
> wrote:
> > > > If I have a proj string representing a coordinate operation, is there
> > > > any way to determine what grids are required for using it when those
> > > > grids are NOT available for use on the system?
> > >
> > > I don't think so. As you find out, PROJ will refuse to instanciate a
> > > coordinate operation from a PROJ string that has a grids not present on
> > > the
> > > system.
> > >
> > > Asthere are not so many ways of expressing grids in PROJ grids, basically
> > > you could look for +grids=xxxxx
> >
> > That's what I was originally thinking could be a fallback option. But
> > then, I can't see any way to get the detailed metadata relating to a
> > grid except for proj_coordoperation_get_grid_used. So even if I
> > captured the grid name from the proj string, I still can't determine
> > the url and package name, etc relating to that grid. (or can I?)
>
> In the C++ API, you have DatabaseContext::lookForGridInfo() that would cover
> this. Would need to be exposed in the C API.

Thanks, PR at https://github.com/OSGeo/proj.4/pull/1494

Nyall
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj