Add ST_TransformPipeline #705

rcoup · 2022-08-29T15:34:30Z

As discussed in Trac #5006, add the ability to use a specific PROJ conversion pipeline instead of the automatically selected one via ST_Transform() (proj_crs_to_crs()).

Thanks to @rouault for all his help at the FOSS4G 2022 code sprint getting it working 😃 Commit messages should be relatively self-explanatory.

geometry ST_TransformPipeline(geom geometry, pipeline text[, to_srid integer])
geometry ST_InverseTransformPipeline(geom geometry, pipeline text[, to_srid integer])

Transformation pipelines are defined using any of the following string formats:

urn:ogc:def:coordinateOperation:AUTHORITY::CODE. Note that a simple AUTHORITY:CODE string does not uniquely identify a coordinate operation: the same code can be used for a CRS definition.
a PROJ pipeline string of the form: +proj=pipeline .... Automatic axis normalisation will not be applied, and if necessary the caller will need to add an additional pipeline step.
concatenated operations of the form: urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618

The SRID of the input geometry is ignored, and the SRID of the output geometry will be set to zero unless a value is provided via the optional to_srid parameter.

When using ST_TransformPipeline() the pipeline is executed in a forward direction. Using ST_InverseTransformPipeline() the pipeline is executed in the inverse direction.

GDA2020 example from the docs:

-- using ST_Transform with automatic selection of a conversion pipeline.
SELECT ST_AsText(ST_Transform('SRID=4939;POINT(143.0 -37.0)'::geometry, 7844)) AS gda2020_auto;

                 gda2020_auto
-----------------------------------------------
 POINT(143.00000635638918 -36.999986706128176)
(1 row)

-- using a defined conversion (EPSG:8447)
SELECT ST_AsText(ST_TransformPipeline('SRID=4939;POINT(143.0 -37.0)'::geometry,
  'urn:ogc:def:coordinateOperation:EPSG::8447')) AS gda2020_code;

                   gda2020_code
----------------------------------------------
 POINT(143.0000063280214 -36.999986718287545)
(1 row)

-- using a PROJ pipeline definition matching EPSG:8447, as returned from
-- 'projinfo -s EPSG:4939 -t EPSG:7844'.
-- NOTE: any 'axisswap' steps must be removed.
SELECT ST_AsText(ST_TransformPipeline('SRID=4939;POINT(143.0 -37.0)'::geometry,
  '+proj=pipeline
   +step +proj=unitconvert +xy_in=deg +xy_out=rad
   +step +proj=hgridshift +grids=au_icsm_GDA94_GDA2020_conformal_and_distortion.tif
   +step +proj=unitconvert +xy_in=rad +xy_out=deg')) AS gda2020_pipeline;

                   gda2020_pipeline
----------------------------------------------
 POINT(143.0000063280214 -36.999986718287545)
(1 row)

robe2 · 2022-08-30T05:08:19Z

@rcoup I have a couple of issues with your commit:

You are missing an availability notice on the new function.

It should read:

<para>Availability: 3.4.0 - Requires PostGIS be compiled with PROJ v6.1+.  Use <xref linkend="PostGIS_Full_Version" /> to confirm you have PROJ support compiled in.</para>

Your tests also don't seem to guard against running on PROJ < 6.1, but that might be okay.
I think come 3.4.0 we might dispense with PROJ < 6.1. As I recall @pramsey is the only one that had qualms about dropping PROJ < 6. @pramsey you still feel on the fence on that.
Your sql functions have too many prototypes and have the wrong availability, they should read

Availability: 3.4.0

not 3.3.0

I would also reduce down to 1 SQL interface of:

geometry ST_TransformPipeline(geom geometry, pipeline text, forward boolean DEFAULT true, to_srid integer DEFAULT 0)

Also our preferred datatype is boolean instead of bool for the sql api. I know we still have some dangling bool's in there. so sorry for that.

rcoup · 2022-08-30T08:44:43Z

@robe2 thanks! Can easily make those changes.

Your tests also don't seem to guard against running on PROJ < 6.1, but that might be okay.

Is this the best example of such a guard, or is there an easier mechanism?

strk · 2022-08-30T12:41:08Z

I agree with Regina about reducing to single signature. Is the "forward" parameter really needed ? I mean: can a pipeline not always be expressed differently, to be backward ? I'm thinking of a function that, taken a pipeline, returns it's "reverse", so for non-forward you'd do: SELECT ST_TransformPipeline( geom, ST_ReverseTransformPipeline(:pipeline), :target_srid ); I keep thinking we should have some Transform _TYPE_ at some point, to do these things.

…

--strk;

rcoup · 2022-08-30T12:50:51Z

Is the "forward" parameter really needed ?
I mean: can a pipeline not always be expressed differently,
to be backward ? I'm thinking of a function that, taken a pipeline,
returns it's "reverse", so for non-forward you'd do:

@rouault can explain better than me, but my understanding is:

the EPSG DB defines only forward conversions (eg: EPSG:8447)
PROJ can perform a transform in either forward or reverse
while if you take the proj pipeline text (+proj=pipeline) and reverse the steps, you'll generally get the same answers, the PROJ pipeline definition can't express everything that the ISO spec / defined CRS do
so having "execute in reverse" is a better approach

rouault · 2022-08-30T13:00:53Z

so having "execute in reverse" is a better approach

what is better or not is a master of taste, but it is indeed more straightforward as implementation. That said, PJ* proj_coordoperation_create_inverse(PJ_CONTEXT *ctx, const PJ *obj) exists and could in theory be used to implement a ST_ReverseTransformPipeline().
To be noted that WKT obtained from such an inversed transformation object is generally not interoperable with other implementations, due to falling back to transformation method names like "Inverse of {forward_method_name}" to indicate that the method should be taken in the reverse path (not sure if there are other implementations than PROJ that accept WKT for transformations anyway...), when the transformation cannot be reversed with an existing foward method. This remarks only matter if PostGIS would have a WKT output of transformation pipelines.

robe2 · 2022-08-30T15:29:30Z

@robe2 thanks! Can easily make those changes.

Your tests also don't seem to guard against running on PROJ < 6.1, but that might be okay.

Is this the best example of such a guard, or is there an easier mechanism?

Our general mechanism is to exclude the test if version is not a specific version.
We have those for GEOS / SFCGAL - such as:

postgis/regress/core/tests.mk.in

Line 151 in a17e612

ifeq ($(shell expr "$(POSTGIS_GEOS_VERSION)" ">=" 30700),1)

But you'd need to add at the top:

POSTGIS_PROJ_VERSION=@POSTGIS_PROJ_VERSION@

I think this variable is picked up at configure time -- here

postgis/configure.ac

Line 898 in a17e612

    
           			POSTGIS_PROJ_VERSION=`$PKG_CONFIG proj --modversion | sed 's/\([[0-9]]\).*\([[0-9]]\).*\([[0-9]]\)/\1\2/'`

So should just work if you add to the top. @strk am I missing anything here?

strk · 2022-08-30T15:43:02Z

I didn't look at the link (sorry) but just make sure you do get a specific user error if trying to use the function, something like "PostGIS was built with proj < $version so this functionlity is not available, please rebuild" and not something like "function not found". It's important that the SQL and C functions do exist, for easier upgrades later (or this is our policy at the moment I think).

robe2 · 2022-08-31T04:46:16Z

I didn't look at the link (sorry) but just make sure you do get a specific user error if trying to use the function, something like "PostGIS was built with proj < $version so this functionlity is not available, please rebuild" and not something like "function not found". It's important that the SQL and C functions do exist, for easier upgrades later (or this is our policy at the moment I think).

@rcoup Regarding what @strk said above. We always define the SQL function even when compiling with a lower version, but the safeguard is put in place and stubbed with an error if PostGIS is compiled with a lower version.

Such as this -- https://github.com/postgis/postgis/blob/master/postgis/lwgeom_transform.c#L155
and then throw an informative error such as this:

https://github.com/postgis/postgis/blob/master/postgis/lwgeom_geos.c#L347

But if I get my wish -- all of this will be unnecessary, cause we'll drop support for Proj < 6.1 so no ifdef will be needed - as noted here:

https://lists.osgeo.org/pipermail/postgis-devel/2022-August/029786.html

BTW opinions on the vote are welcome from anybody and will sway the PSC decision.

rcoup · 2022-08-31T10:04:01Z

But if I get my wish -- all of this will be unnecessary, cause we'll drop support for Proj < 6.1 so no ifdef will be needed - as noted here:

https://lists.osgeo.org/pipermail/postgis-devel/2022-August/029786.html

BTW opinions on the vote are welcome from anybody and will sway the PSC decision.

Ah, I'm only subscribed to postgis-users, so I didn't see this. Certainly makes sense to me: seems unlikely people will want latest-greatest PostGIS from 2022 but aren't prepared to bring along a vaguely modern PROJ (6.1 was released May 2019).

Will wait to update this PR until that decision is settled, but if it passes it's probably worth rebasing this on a PROJ-removal PR. I'm happy to do that if I can find some time.

rcoup · 2022-08-31T11:54:44Z

could in theory be used to implement a ST_ReverseTransformPipeline()

Feel to me like ST_TransformPipeline() is already extra-for-experts territory: PostGIS has got a long way with ST_Transform(), and the people dealing with more complicated transformations probably know a bit more about what they're trying to do.

robe2 · 2022-09-02T22:51:13Z

Ah, I'm only subscribed to postgis-users, so I didn't see this. Certainly makes sense to me: seems unlikely people will want latest-greatest PostGIS from 2022 but aren't prepared to bring along a vaguely modern PROJ (6.1 was released May 2019).

Will wait to update this PR until that decision is settled, but if it passes it's probably worth rebasing this on a PROJ-removal PR. I'm happy to do that if I can find some time.

Motion passed.
https://lists.osgeo.org/pipermail/postgis-devel/2022-September/029793.html

But let me first start removing bots testing with < 6.1 and then we can rerun.
We also voted to remove PG 11 too.

robe2 · 2022-09-03T00:22:39Z

@rcoup Okay done ripping out support of PROJ < 6.1 and PostgreSQL < 12. So now you can rebase.

https://trac.osgeo.org/postgis/ticket/5229

robe2 · 2022-10-01T02:22:27Z

@rcoup Just checking if you are ready to rebase and reapply. As mentioned no need for the conditional proj since we no longer support PROJ < 6.1 for 3.4

strk · 2022-10-11T07:09:43Z

Thanks for the explanation. Would encoding inverse in the function name be more elegant then ? ST_ReverseTransformPipeline() ?

rcoup · 2022-10-18T14:14:18Z

@robe2 @strk I've made those changes now.

I went with adding ST_InverseTransformPipeline() instead of the forward parameter.

Feedback welcome 😃

While the code was there for doing angular unit conversion in ptarray_transform() on both input & output (which is correct and matches proj's `cct` behaviour); converted points were never written back to the coordinate array. In practise this only applies in edge cases - normally with transforms created via proj_crs_to_crs() PROJ does it automatically, but is exercised by the upcoming ST_TransformPipeline() changes.

EPSG/etc define forward pipeline definitions (eg: EPSG:16031), but there's no inverse definitions, and reversing the PROJ pipeline steps in a proj-string pipeline definition will go via a very different codepath in PROJ than changing the pipeline direction (though in general will produce the same result coordinates). Add the concept of a transform direction to the LWPROJ structure. For the existing CRS to CRS conversion used with ST_Transform() we always do a forward transform, but with the new transform pipeline method the user will also have the option to execute inverse pipelines.

There can be reasons to apply a specific CRS transformation pipeline, rather than leaving it to PROJ to select a candidate pipeline based on the source & target CRS. Introduce `lwproj_from_str_pipeline()` to instantiate a LWPROJ operation from a PROJ pipeline definition, using either of these formats: - a defined coordinate operation: `urn:ogc:def:coordinateOperation:AUTHORITY::CODE` - a PROJ pipeline: `+proj=pipeline ...` - concatenated operations: `urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618` - any other pipeline accepted via `proj_create()` Pipelines can be instantiated in a forward (default) or inverse direction. Introduce `lwgeom_transform_pipeline()` to both instantiate and execute such pipeline transformations.

Add the capability to perform CRS transformations using a specific/defined pipeline rather than leaving PROJ to automatically select the transform to use via ST_Transform(). There are two forms of this function: geometry ST_TransformPipeline(geom geometry, pipeline text[, to_srid integer]) geometry ST_InverseTransformPipeline(geom geometry, pipeline text[, to_srid integer]) Returns a new geometry with its coordinates transformed to a different coordinate reference system using a defined coordinate transformation pipeline. Transformation pipelines are defined using any of the following string formats: - `urn:ogc:def:coordinateOperation:AUTHORITY::CODE`. Note that a simple `AUTHORITY:CODE` string does not uniquely identify a coordinate operation: the same code can be used for a CRS definition. - a PROJ pipeline string of the form: `+proj=pipeline ...`. Automatic axis normalisation will not be applied, and if necessary the caller will need to add an additional pipeline step. - concatenated operations of the form: `urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618` The SRID of the input geometry is ignored, and the SRID of the output geometry will be set to zero unless a value is provided via the optional to_srid parameter. When using `ST_TransformPipeline()` the pipeline is executed in a forward direction. Using `ST_InverseTransformPipeline()` the pipeline is executed in the inverse direction.

Use spaces consistently

rcoup · 2022-10-20T12:43:12Z

ok, so make check seems to work locally (and exits with 0). But looking closer:

Loading PostGIS into 'postgis_reg-3.4'
 failed (Error encountered loading /Users/rcoup/code/postgis/regress/00-regress-install/share/contrib/postgis/postgis_comments.sql: /tmp/pgis_reg/regress_log)

$ cat /tmp/pgis_reg/regress_log
NOTICE:  schema "public" already exists, skipping
NOTICE:  schema "public" already exists, skipping
psql:/Users/rcoup/code/postgis/regress/00-regress-install/share/contrib/postgis/postgis_comments.sql:282: ERROR:  function st_transformpipeline(geometry, text) does not exist
NOTICE:  schema "public" already exists, skipping
...

Which leads to a syntax issue in the docs. Easy fix 👍

I can't spot the area in the twisty maze of Makefiles where it's actually ignoring the error.

strk · 2022-10-20T13:08:04Z

Could you please file a trac ticket for the false negative in make check ? I suspect it's to do with run_test.pl itself not catching psql error (or new versions of psql returning errors in a different way)

See #705 (comment)

strk · 2022-10-20T13:26:35Z

The false success should be fixed with 0090048

rcoup · 2022-10-20T14:21:35Z

@strk thanks. I re-added a similar error locally & make check fails (returns 2) now. 👍

robe2 requested review from pramsey and strk August 30, 2022 04:45

rcoup force-pushed the rc-transform-pipeline-5006 branch 2 times, most recently from d7011f3 to 78f898d Compare October 18, 2022 14:11

rcoup added 7 commits October 20, 2022 12:20

regress: add tests for transform pipelines

1bc9ace

doc: reference_srs: convert whitespace

01bbb2c

Use spaces consistently

doc: add documentation for transform pipelines

f9db6ac

rcoup force-pushed the rc-transform-pipeline-5006 branch from 78f898d to f9db6ac Compare October 20, 2022 12:42

strk added a commit that referenced this pull request Oct 20, 2022

Fix false success when db preparation fails

0090048

See #705 (comment)

strk merged commit 961552f into postgis:master Oct 21, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add ST_TransformPipeline #705

Add ST_TransformPipeline #705

rcoup commented Aug 29, 2022 •

edited

robe2 commented Aug 30, 2022 •

edited

rcoup commented Aug 30, 2022

strk commented Aug 30, 2022 via email

rcoup commented Aug 30, 2022

rouault commented Aug 30, 2022

robe2 commented Aug 30, 2022

strk commented Aug 30, 2022 via email

robe2 commented Aug 31, 2022

rcoup commented Aug 31, 2022

rcoup commented Aug 31, 2022

robe2 commented Sep 2, 2022

robe2 commented Sep 3, 2022

robe2 commented Oct 1, 2022

strk commented Oct 11, 2022 via email

rcoup commented Oct 18, 2022

rcoup commented Oct 20, 2022

strk commented Oct 20, 2022

strk commented Oct 20, 2022

rcoup commented Oct 20, 2022

Add ST_TransformPipeline #705

Add ST_TransformPipeline #705

Conversation

rcoup commented Aug 29, 2022 • edited

robe2 commented Aug 30, 2022 • edited

rcoup commented Aug 30, 2022

strk commented Aug 30, 2022 via email

rcoup commented Aug 30, 2022

rouault commented Aug 30, 2022

robe2 commented Aug 30, 2022

strk commented Aug 30, 2022 via email

robe2 commented Aug 31, 2022

rcoup commented Aug 31, 2022

rcoup commented Aug 31, 2022

robe2 commented Sep 2, 2022

robe2 commented Sep 3, 2022

robe2 commented Oct 1, 2022

strk commented Oct 11, 2022 via email

rcoup commented Oct 18, 2022

rcoup commented Oct 20, 2022

strk commented Oct 20, 2022

strk commented Oct 20, 2022

rcoup commented Oct 20, 2022

rcoup commented Aug 29, 2022 •

edited

robe2 commented Aug 30, 2022 •

edited