ST_Transform: Support PROJ pipelines

[cdestigter]

I hoped to be able to use a ​PROJ pipeline string to apply a non-default transform between two coordinate systems.

Usage might look like this:

select st_asewkt(
    ST_TransformPipeline(
        'SRID=4283;POINT(143.0 -37.0)'::GEOMETRY,
        '+proj=pipeline +step +proj=axisswap +order=2,1 +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 +step +proj=axisswap +order=2,1'
    )
);

Ideally the EPSG code for the transform would be usable too, e.g.

select st_asewkt(
    ST_TransformPipeline(
        'SRID=4283;POINT(143.0 -37.0)'::GEOMETRY,
        'EPSG:8447'
    )
);

I'm not sure if supporting pipelines defined by other authorities is a concern (maybe?), and if so I'd probably steer clear of supporting an integer parameter for the EPSG code, otherwise we'll need a lookup table:

-- this would probably need a lookup table, maybe best to avoid?
select st_asewkt(
    ST_TransformPipeline(
        'SRID=4283;POINT(143.0 -37.0)'::GEOMETRY,
        8447
    )
);

[strk]

I guess we might want to also allow a third parameter specifying a SRID for the target (defaulting to 0)

[rcoup]

@strk — is this what you're thinking?

geometry ST_TransformPipeline(geometry geom, text pipeline);
geometry ST_TransformPipeline(geometry geom, text pipeline, integer to_srid);

Where pipeline can be anything ​`proj_create()` accepts that produces a coordinate operation / transformation:

• EPSG:1671 (referencing a transformation)
• +proj=pipeline ...
• urn:ogc:def:coordinateOperation:EPSG::1671
• concatenated operations urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618

I agree with Craig, I don't think integer codes are a great idea. For a start, spatial_ref_sys only holds CRS codes, not transformations.

Which would, underneath:

1. proj_create(pipeline)
2. check the result is a coordinate transform and not a CRS (!proj_is_crs())
3. use it via proj_trans*()
4. ST_SetSRID() on the result if we have a to_srid parameter

Things we might need to figure out / be careful with: PostGIS always calls proj_normalize_for_visualization() since it uses an X,Y axis order convention. Maybe we document that any axis swapping happening in the pipeline is subject to an additional X,Y normalisation at the end?

@rouault any comments from a Proj perspective?

[rouault]

Using a string to express is probably the best approach here (you could also pass transformations as WKT or PROJJSON, although that would be super verbose). The integer as being a code of an EPSG registered transformation would of course restrict it to EPSG, but also would prevent pipelines from/to a projected CRS. I'd note that proj_create() with a "EPSG:XXXX" string is an ambiguous operation. As there are EPSG objects with the same code but different type (CRS, transformations, etc.), proj_create() makes the arbitrary choice to return the CRS in case of ambiguity. So EPSG:XXXX might return if you're lucky a transformation, but if at the next update of the database that code is also attributed to a new CRS, game over. The URN syntax "urn:ogc:def:coordinateOperation:EPSG::1671" would have to be used. PostGIs could still decide to parse FOO:BAR and transform that to urn:ogc:def:coordinateOperation:FOO::BAR in the context of that function.

The to_srid argument is probably useful indeed. If using a registered transformation, you can query the target CRS on it with the proj API, but that don't necessarily give you the corresponding PostGIS SRID if the target is not a EPSG CRS. And if you use a PROJ pipeline, PROJ won't be able to know the target CRS at all.

The axis order might be a tricky. If the user pass a PROJ pipeline that does axis switching (like the GDA94 to GDA2020), then I'm not sure what PostGIS can do. The proj_normalize_for_visualization() accept PROJ coordinate operation objects, but that won't work with a PROJ pipeline, as it, at least currently, requires a PROJ coordinate operation for which the source and target CRS are known, and as mentionned above this is not the case for a PROJ pipeline. Three solutions I can think of:

• document that users must strip any leading and trailing +proj=axisswap steps from their pipelines
• or implement that in PostGIS
• or implement that in proj_normalize_for_visualization()

[strk]

My idea of an integer ID was NOT to query spatial_ref_sys but a different, new table. Sometihng like spatial_ref_sys_transforms or similar. Such table could have a column expressing source SRID, a column expressiong target SRID and a proj4 string expressing the actual transformation/pipeline. The key to access it could even not be an integer but a descriptive name.

I guess users could build such table theirselves and query it dynamically with their queries, but if we had an official table maybe we could do some neat things like verifying the input geometry has a valid source SRID for that transformation… (just an idea).

In any case the free-form signature would still be the workhorse so still the best first step.

[rcoup]

Thanks

I'd note that proj_create() with a "EPSG:XXXX" string is an ambiguous operation. As there are EPSG objects with the same code but different type (CRS, transformations, etc.)

Ah, of course there are. <sarcasm>Because we can't waste unique numbers, there's so few of them…</sarcasm>

The URN syntax "urn:ogc:def:coordinateOperation:EPSG::1671" would have to be used. PostGIs could still decide to parse FOO:BAR and transform that to urn:ogc:def:coordinateOperation:FOO::BAR in the context of that function.

PostGIS uses 1234 in context of CRS lookups, so even EPSG:1234 would be a "new" way of expressing CRS-related things. Might as well just document "use urn:ogc:def:coordinateOperation:FOO::BAR or use a pipeline/WKT string".

OT: @rouault is un-overloading proj_create() maybe a sensible API progression for Proj? To differentiate "EPSG:1234 and I want a CRS" from "EPSG:1234 and I want a transformation"? ie: "I expect an object of type X from this string".

My idea of an integer ID was NOT to query spatial_ref_sys but a different, new table. Sometihng like spatial_ref_sys_transforms or similar. Such table could have a column expressing source SRID, a column expressiong target SRID and a proj4 string expressing the actual transformation/pipeline. The key to access it could even not be an integer but a descriptive name.

Seems to me like a lot of (ongoing) effort when 99.9% of the time Proj will use the right transform automatically? From my point of view this is more being able to escape the automatic route if there are some special/specific needs. As you say, can always be enhanced in the future with an integer lookup if it proves popular.

The proj_normalize_for_visualization() accept PROJ coordinate operation objects, but that won't work with a PROJ pipeline, as it, at least currently, requires a PROJ coordinate operation for which the source and target CRS are known, and as mentionned above this is not the case for a PROJ pipeline.

But if the pipeline was specified as urn:ogc:def:coordinateOperation:EPSG::8447 Proj does know the source/target CRS and it'll work?

Q2: Is source/target CRS settable on a transformation after creating it? If so and the user passes a pipeline string, PostGIS does know the source CRS, and the user can provide a target CRS ala ST_TransformPipeline(geom, '+proj=pipeline ...', 1234) (after mapping SRIDs to CRS via AUTH:CODE).

I guess might be complicated somewhat with horizontal+vertical compound CRS/transforms, though normalize_for_visualization() only applies to horizontal anyway.

[rouault]

@rouault is un-overloading proj_create() maybe a sensible API progression for Proj? To differentiate "EPSG:1234 and I want a CRS" from "EPSG:1234 and I want a transformation"? ie: "I expect an object of type X from this string".

This is ​https://proj.org/development/reference/functions.html#c.proj_create_from_database where you specify (authority_name, code, object_type).

But if the pipeline was specified as urn:ogc:def:coordinateOperation:EPSG::8447 Proj does know the source/target CRS and it'll work?

Yes, that will work with anything that is not a PROJ pipeline string

Is source/target CRS settable on a transformation after creating it?

No, we could potentially have a proj_alter_source_crs / target_crs to do that, but they are not yet available

[rcoup]

ok, so revising a little further after the above discussions, how about this as a first cut to make it possible to use:

geometry ST_TransformPipeline(geometry geom, text pipeline)
geometry ST_TransformPipeline(geometry geom, text pipeline, integer to_srid)

Where pipeline is documented as:

• urn:ogc:def:coordinateOperation:AUTHORITY::CODE
• a Proj pipeline (+proj=pipeline ...). Document that automatic axis normalisation doesn't happen, the caller would need to do it manually.
• concatenated operations: urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618
• though anything proj_create() accepts would be passed through.

Which would roughly:

1. proj_create(pipeline)
2. check the result is a coordinate transform and not a CRS (!proj_is_crs())
3. if the source & target CRS of the operation are known (proj_get_id_auth_name()/proj_get_id_auth_code()) then call normalize_for_visualization()
4. use it via proj_trans*()
5. set the Postgis srid on the resulting geometry if we have a to_srid parameter

[strk]

Sounds good to me. I'm only wondering at this point if re-using ST_Transform(geometry, text) and ST_Transform(geometry, text, int) as signatures would be problematic. Since we can check with proj_is_crs() I guess we could have that drive the behaviour ?

For reference, here's the documentation of the current signatures: ​https://postgis.net/docs/ST_Transform.html

The only problem *might* be the parameter *name* (ie: to_proj would be misleading)

[rcoup]

Implementation PR at ​https://github.com/postgis/postgis/pull/705

[Regina Obe <lr@…>]

In 961552f/git:

Merge remote-tracking branch 'rcoup/rc-transform-pipeline-5006'
Support for ST_Transform: Support PROJ pipelines
Closes #5006
Closes GH705