[postgis-tickets] [SCM] PostGIS branch master updated. 3.3.0rc2-241-g961552f6f
git at osgeo.org
git at osgeo.org
Fri Oct 21 09:41:51 PDT 2022
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, master has been updated
via 961552f6f8191db87448ab6c93c1cd750fdc34da (commit)
via 347f2ea746b41b030119ed81676f3539f4ce961a (commit)
via f9db6acfdd5aa02b86418282c13d997a67eedfa5 (commit)
via 01bbb2c39f0dbe599d46bb04f77920d7811cb495 (commit)
via 1bc9ace3593bd3711cc388682b9891f68652418e (commit)
via d33ebde11f679cbc100aecf6d98c608d211820c6 (commit)
via 5606caad8fde789045d6b7cdaca5340cb14ccc10 (commit)
via 3d6c8256034aec3dcb91d7c36a4bbe86ff2697ef (commit)
via 3cefc95100088f98c69b7f5984d119728fb88f9e (commit)
from 7f54d8ce0805917aad05ea4fb3d57bc6afaa55e4 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 961552f6f8191db87448ab6c93c1cd750fdc34da
Merge: 347f2ea74 f9db6acfd
Author: Regina Obe <lr at pcorp.us>
Date: Fri Oct 21 12:35:48 2022 -0400
Merge remote-tracking branch 'rcoup/rc-transform-pipeline-5006'
Support for ST_Transform: Support PROJ pipelines
Closes #5006
Closes GH705
diff --cc NEWS
index b2d931c69,c1a1d6114..065b0c536
--- a/NEWS
+++ b/NEWS
@@@ -6,7 -6,7 +6,8 @@@ xxxx/xx/x
* New Features *
- New install-extension-upgrades command in postgis script (Sandro Santilli)
- - #5257, #5261 Support changes for PG16 (Regina Obe)
- - #5257, Remove MemoryContextContains use for PG16 (Regina Obe)
++ - #5257, #5261, Support changes for PG16 (Regina Obe)
++ - #5006, GH705, ST_Transform: Support PROJ pipelines (Robert Coup, Koordinates)
* Enhancements *
- #5194, do not update system catalogs from postgis_extensions_upgrade (Sandro Santilli)
diff --cc doc/introduction.xml
index 1ed7fa4f2,1ed7fa4f2..61db2a0c4
--- a/doc/introduction.xml
+++ b/doc/introduction.xml
@@@ -299,6 -299,6 +299,7 @@@
<member>Ralph Mason</member>
<member>Rémi Cura</member>
<member>Richard Greenwood</member>
++ <member>Robert Coup</member>
<member>Roger Crew</member>
<member>Ron Mayer</member>
<member>Sebastiaan Couwenberg</member>
commit 347f2ea746b41b030119ed81676f3539f4ce961a
Author: Regina Obe <lr at pcorp.us>
Date: Fri Oct 21 12:22:07 2022 -0400
Change TWKB arg names to agree with the actual names in the function
diff --git a/doc/reference_output.xml b/doc/reference_output.xml
index b11955cae..52b192603 100644
--- a/doc/reference_output.xml
+++ b/doc/reference_output.xml
@@ -1549,22 +1549,22 @@ FROM mvtgeom;
<funcsynopsis>
<funcprototype>
<funcdef>bytea <function>ST_AsTWKB</function></funcdef>
- <paramdef><type>geometry </type> <parameter>g1</parameter></paramdef>
- <paramdef><type>integer </type> <parameter>decimaldigits_xy=0</parameter></paramdef>
- <paramdef><type>integer </type> <parameter>decimaldigits_z=0</parameter></paramdef>
- <paramdef><type>integer </type> <parameter>decimaldigits_m=0</parameter></paramdef>
- <paramdef><type>boolean </type> <parameter>include_sizes=false</parameter></paramdef>
- <paramdef><type>boolean </type> <parameter>include_bounding boxes=false</parameter></paramdef>
+ <paramdef><type>geometry </type> <parameter>geom</parameter></paramdef>
+ <paramdef><type>integer </type> <parameter>prec=0</parameter></paramdef>
+ <paramdef><type>integer </type> <parameter>prec_z=0</parameter></paramdef>
+ <paramdef><type>integer </type> <parameter>prec_m=0</parameter></paramdef>
+ <paramdef><type>boolean </type> <parameter>with_sizes=false</parameter></paramdef>
+ <paramdef><type>boolean </type> <parameter>with_boxes=false</parameter></paramdef>
</funcprototype>
<funcprototype>
<funcdef>bytea <function>ST_AsTWKB</function></funcdef>
- <paramdef><type>geometry[] </type> <parameter>geometries</parameter></paramdef>
- <paramdef><type>bigint[] </type> <parameter>unique_ids</parameter></paramdef>
- <paramdef><type>integer </type> <parameter>decimaldigits_xy=0</parameter></paramdef>
- <paramdef><type>integer </type> <parameter>decimaldigits_z=0</parameter></paramdef>
- <paramdef><type>integer </type> <parameter>decimaldigits_m=0</parameter></paramdef>
- <paramdef><type>boolean </type> <parameter>include_sizes=false</parameter></paramdef>
- <paramdef><type>boolean </type> <parameter>include_bounding_boxes=false</parameter></paramdef>
+ <paramdef><type>geometry[] </type> <parameter>geom</parameter></paramdef>
+ <paramdef><type>bigint[] </type> <parameter>ids</parameter></paramdef>
+ <paramdef><type>integer </type> <parameter>prec=0</parameter></paramdef>
+ <paramdef><type>integer </type> <parameter>prec_z=0</parameter></paramdef>
+ <paramdef><type>integer </type> <parameter>prec_m=0</parameter></paramdef>
+ <paramdef><type>boolean </type> <parameter>with_sizes=false</parameter></paramdef>
+ <paramdef><type>boolean </type> <parameter>with_boxes=false</parameter></paramdef>
</funcprototype>
</funcsynopsis>
</refsynopsisdiv>
commit f9db6acfdd5aa02b86418282c13d997a67eedfa5
Author: Robert Coup <robert at coup.net.nz>
Date: Mon Aug 29 16:25:24 2022 +0100
doc: add documentation for transform pipelines
diff --git a/doc/reference_srs.xml b/doc/reference_srs.xml
index 2436d3d55..f03776c7d 100644
--- a/doc/reference_srs.xml
+++ b/doc/reference_srs.xml
@@ -171,6 +171,9 @@ SRID=3785;POINT(-13732990.8753491 6178458.96425423)
of a geometry from one spatial reference system to another, while ST_SetSRID() simply changes the SRID identifier of
the geometry.</para>
+ <para>ST_Transform automatically selects a suitable conversion pipeline given the source and target spatial
+ reference systems. To use a specific conversion method, use <xref linkend="ST_TransformPipeline" />.</para>
+
<note>
<para>Requires PostGIS be compiled with PROJ support. Use <xref linkend="PostGIS_Full_Version" /> to confirm you have PROJ support compiled in.</para>
</note>
@@ -270,7 +273,147 @@ CREATE INDEX idx_geom_26986_parcels
<refsection>
<title>See Also</title>
- <para><xref linkend="spatial_ref_sys" />, <xref linkend="ST_SetSRID" />, <xref linkend="ST_SRID" />, <xref linkend="UpdateGeometrySRID"/></para>
+ <para><xref linkend="spatial_ref_sys" />, <xref linkend="ST_SetSRID" />, <xref linkend="ST_SRID" />, <xref linkend="UpdateGeometrySRID"/>, <xref linkend="ST_TransformPipeline" /></para>
+ </refsection>
+ </refentry>
+
+ <refentry id="ST_TransformPipeline">
+ <refnamediv>
+ <refname>ST_TransformPipeline</refname>
+
+ <refpurpose>Return a new geometry with coordinates transformed to
+ a different spatial reference system using a defined coordinate
+ transformation pipeline.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_TransformPipeline</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>g1</parameter></paramdef>
+ <paramdef><type>text </type> <parameter>pipeline</parameter></paramdef>
+ <paramdef choice="opt"><type>integer </type> <parameter>to_srid</parameter></paramdef>
+ </funcprototype>
+
+ <funcprototype>
+ <funcdef>geometry <function>ST_InverseTransformPipeline</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>geom</parameter></paramdef>
+ <paramdef><type>text </type> <parameter>pipeline</parameter></paramdef>
+ <paramdef choice="opt"><type>integer </type> <parameter>to_srid</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>
+ Return a new geometry with coordinates transformed to a different spatial reference system
+ using a defined coordinate transformation pipeline.
+ </para>
+
+ <para>
+ Transformation pipelines are defined using any of the following string formats:
+ <itemizedlist>
+ <listitem>
+ <para>
+ <varname>urn:ogc:def:coordinateOperation:AUTHORITY::CODE</varname>. Note that a simple
+ <varname>EPSG:CODE</varname> string does not uniquely identify a coordinate operation:
+ the same EPSG code can be used for a CRS definition.
+ </para>
+ </listitem>
+
+ <listitem>
+ <para>
+ A PROJ pipeline string of the form: <varname>+proj=pipeline ...</varname>. Automatic
+ axis normalisation will not be applied, and if necessary the caller will need to add an
+ additional pipeline step, or remove <varname>axisswap</varname> steps.
+ </para>
+ </listitem>
+
+ <listitem>
+ <para>
+ Concatenated operations of the form: <varname>urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618</varname>.
+ </para>
+ </listitem>
+ </itemizedlist>
+ </para>
+
+ <para>
+ 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 <varname>to_srid</varname> parameter. When
+ using `ST_TransformPipeline()` the pipeline is executed in a forward direction. Using
+ `ST_InverseTransformPipeline()` the pipeline is executed in the inverse direction.</para>
+
+ <para>Transforms using pipelines are a specialised version of <xref linkend="ST_Transform" />.
+ In most cases `ST_Transform` will choose the correct operations to convert between coordinate
+ systems, and should be preferred.</para>
+
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+ <para>Change WGS 84 long lat to UTM 31N using the EPSG:16031 conversion</para>
+ <programlisting>
+-- Forward direction
+SELECT ST_AsText(ST_TransformPipeline('SRID=4326;POINT(2 49)'::geometry,
+ 'urn:ogc:def:coordinateOperation:EPSG::16031') AS utm_geom);
+
+ utm_geom
+--------------------------------------------
+ POINT(426857.9877165967 5427937.523342293)
+(1 row)
+
+-- Inverse direction
+SELECT ST_AsText(ST_InverseTransformPipeline('POINT(426857.9877165967 5427937.523342293)'::geometry,
+ 'urn:ogc:def:coordinateOperation:EPSG::16031')) AS wgs_geom;
+
+ wgs_geom
+----------------------------
+ POINT(2 48.99999999999999)
+(1 row)
+ </programlisting>
+
+ <para>GDA2020 example.</para>
+ <programlisting>
+-- 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)
+ </programlisting>
+ </refsection>
+
+ <!-- Optionally add a "See Also" section -->
+ <refsection>
+ <title>See Also</title>
+
+ <para><xref linkend="ST_Transform" /></para>
</refsection>
</refentry>
diff --git a/doc/using_postgis_dataman.xml b/doc/using_postgis_dataman.xml
index 8d7f6ea7e..61350ba3f 100644
--- a/doc/using_postgis_dataman.xml
+++ b/doc/using_postgis_dataman.xml
@@ -1381,7 +1381,7 @@ SELECT valid, reason, ST_AsText(location) AS location
<link linkend="Spatial_Relationships">relationship</link> functions)
require that the input geometries are in the same spatial reference system (have the same SRID).
Geometry data can be transformed into a different spatial reference system using
- <xref linkend="ST_Transform" />.
+ <xref linkend="ST_Transform" /> and <xref linkend="ST_TransformPipeline" />.
Geometry returned from functions has the same SRS as the input geometries.
</para>
commit 01bbb2c39f0dbe599d46bb04f77920d7811cb495
Author: Robert Coup <robert at coup.net.nz>
Date: Mon Aug 29 16:23:48 2022 +0100
doc: reference_srs: convert whitespace
Use spaces consistently
diff --git a/doc/reference_srs.xml b/doc/reference_srs.xml
index b6aa6d118..2436d3d55 100644
--- a/doc/reference_srs.xml
+++ b/doc/reference_srs.xml
@@ -2,130 +2,130 @@
<sect1 id="SRS_Functions">
<sect1info>
<abstract>
- <para>These functions work with the Spatial Reference System of geometries
+ <para>These functions work with the Spatial Reference System of geometries
as defined in the <varname>spatial_ref_sys</varname> table.</para>
</abstract>
</sect1info>
- <title>Spatial Reference System Functions</title>
-
- <refentry id="ST_SetSRID">
- <refnamediv>
- <refname>ST_SetSRID</refname>
-
- <refpurpose>Set the SRID on a geometry.</refpurpose>
- </refnamediv>
-
- <refsynopsisdiv>
- <funcsynopsis>
- <funcprototype>
- <funcdef>geometry <function>ST_SetSRID</function></funcdef>
-
- <paramdef><type>geometry </type>
- <parameter>geom</parameter></paramdef>
-
- <paramdef><type>integer </type>
- <parameter>srid</parameter></paramdef>
- </funcprototype>
- </funcsynopsis>
- </refsynopsisdiv>
-
- <refsection>
- <title>Description</title>
-
- <para>Sets the SRID on a geometry to a particular integer value.
- Useful in constructing bounding boxes for queries.</para>
-
- <note>
- <para>This function does not transform the geometry coordinates in any way -
- it simply sets the meta data defining the spatial reference system the geometry is assumed to be in.
- Use <xref linkend="ST_Transform"/> if you want to transform the
- geometry into a new projection.</para>
- </note>
- <para>&sfs_compliant;</para>
- <para>&curve_support;</para>
- </refsection>
-
- <refsection>
- <title>Examples</title>
- <para>-- Mark a point as WGS 84 long lat --</para>
- <programlisting>SELECT ST_SetSRID(ST_Point(-123.365556, 48.428611),4326) As wgs84long_lat;
+ <title>Spatial Reference System Functions</title>
+
+ <refentry id="ST_SetSRID">
+ <refnamediv>
+ <refname>ST_SetSRID</refname>
+
+ <refpurpose>Set the SRID on a geometry.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_SetSRID</function></funcdef>
+
+ <paramdef><type>geometry </type>
+ <parameter>geom</parameter></paramdef>
+
+ <paramdef><type>integer </type>
+ <parameter>srid</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>Sets the SRID on a geometry to a particular integer value.
+ Useful in constructing bounding boxes for queries.</para>
+
+ <note>
+ <para>This function does not transform the geometry coordinates in any way -
+ it simply sets the meta data defining the spatial reference system the geometry is assumed to be in.
+ Use <xref linkend="ST_Transform"/> if you want to transform the
+ geometry into a new projection.</para>
+ </note>
+ <para>&sfs_compliant;</para>
+ <para>&curve_support;</para>
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+ <para>-- Mark a point as WGS 84 long lat --</para>
+ <programlisting>SELECT ST_SetSRID(ST_Point(-123.365556, 48.428611),4326) As wgs84long_lat;
-- the ewkt representation (wrap with ST_AsEWKT) -
SRID=4326;POINT(-123.365556 48.428611)
- </programlisting>
- <para>-- Mark a point as WGS 84 long lat and then transform to web mercator (Spherical Mercator) --</para>
- <programlisting>SELECT ST_Transform(ST_SetSRID(ST_Point(-123.365556, 48.428611),4326),3785) As spere_merc;
+ </programlisting>
+ <para>-- Mark a point as WGS 84 long lat and then transform to web mercator (Spherical Mercator) --</para>
+ <programlisting>SELECT ST_Transform(ST_SetSRID(ST_Point(-123.365556, 48.428611),4326),3785) As spere_merc;
-- the ewkt representation (wrap with ST_AsEWKT) -
SRID=3785;POINT(-13732990.8753491 6178458.96425423)
- </programlisting>
- </refsection>
-
- <refsection>
- <title>See Also</title>
-
- <para><xref linkend="spatial_ref_sys" />, <xref linkend="ST_SRID"/>, <xref linkend="ST_Transform"/>, <xref linkend="UpdateGeometrySRID"/></para>
- </refsection>
-
- </refentry>
-
- <refentry id="ST_SRID">
- <refnamediv>
- <refname>ST_SRID</refname>
- <refpurpose>Returns the spatial reference identifier for a geometry.</refpurpose>
- </refnamediv>
-
- <refsynopsisdiv>
- <funcsynopsis>
- <funcprototype>
- <funcdef>integer <function>ST_SRID</function></funcdef>
- <paramdef><type>geometry </type> <parameter>g1</parameter></paramdef>
- </funcprototype>
- </funcsynopsis>
- </refsynopsisdiv>
-
- <refsection>
- <title>Description</title>
-
- <para>Returns the spatial reference identifier for the ST_Geometry as defined in spatial_ref_sys table. <xref linkend="spatial_ref_sys" /></para>
- <para><note><para>spatial_ref_sys
- table is a table that catalogs all spatial reference systems known to PostGIS and is used for transformations from one spatial
- reference system to another. So verifying you have the right spatial reference system identifier is important if you plan to ever transform your geometries.</para></note></para>
- <para>&sfs_compliant; s2.1.1.1</para>
- <para>&sqlmm_compliant; SQL-MM 3: 5.1.5</para>
- <para>&curve_support;</para>
-
- </refsection>
-
- <refsection>
- <title>Examples</title>
-
- <programlisting>SELECT ST_SRID(ST_GeomFromText('POINT(-71.1043 42.315)',4326));
- --result
- 4326
- </programlisting>
- </refsection>
- <refsection>
- <title>See Also</title>
-
- <para><xref linkend="spatial_ref_sys" />, <xref linkend="ST_SetSRID" />, <xref linkend="ST_Transform" />, <xref linkend="RT_ST_SRID" />, <xref linkend="TG_ST_SRID" /></para>
- </refsection>
- </refentry>
-
- <refentry id="ST_Transform">
- <refnamediv>
- <refname>ST_Transform</refname>
-
- <refpurpose>Return a new geometry with coordinates transformed to
- a different spatial reference system.</refpurpose>
- </refnamediv>
-
- <refsynopsisdiv>
- <funcsynopsis>
- <funcprototype>
- <funcdef>geometry <function>ST_Transform</function></funcdef>
- <paramdef><type>geometry </type> <parameter>g1</parameter></paramdef>
- <paramdef><type>integer </type> <parameter>srid</parameter></paramdef>
- </funcprototype>
+ </programlisting>
+ </refsection>
+
+ <refsection>
+ <title>See Also</title>
+
+ <para><xref linkend="spatial_ref_sys" />, <xref linkend="ST_SRID"/>, <xref linkend="ST_Transform"/>, <xref linkend="UpdateGeometrySRID"/></para>
+ </refsection>
+
+ </refentry>
+
+ <refentry id="ST_SRID">
+ <refnamediv>
+ <refname>ST_SRID</refname>
+ <refpurpose>Returns the spatial reference identifier for a geometry.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>integer <function>ST_SRID</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>g1</parameter></paramdef>
+ </funcprototype>
+ </funcsynopsis>
+ </refsynopsisdiv>
+
+ <refsection>
+ <title>Description</title>
+
+ <para>Returns the spatial reference identifier for the ST_Geometry as defined in spatial_ref_sys table. <xref linkend="spatial_ref_sys" /></para>
+ <para><note><para>spatial_ref_sys
+ table is a table that catalogs all spatial reference systems known to PostGIS and is used for transformations from one spatial
+ reference system to another. So verifying you have the right spatial reference system identifier is important if you plan to ever transform your geometries.</para></note></para>
+ <para>&sfs_compliant; s2.1.1.1</para>
+ <para>&sqlmm_compliant; SQL-MM 3: 5.1.5</para>
+ <para>&curve_support;</para>
+
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+
+ <programlisting>SELECT ST_SRID(ST_GeomFromText('POINT(-71.1043 42.315)',4326));
+ --result
+ 4326
+ </programlisting>
+ </refsection>
+ <refsection>
+ <title>See Also</title>
+
+ <para><xref linkend="spatial_ref_sys" />, <xref linkend="ST_SetSRID" />, <xref linkend="ST_Transform" />, <xref linkend="RT_ST_SRID" />, <xref linkend="TG_ST_SRID" /></para>
+ </refsection>
+ </refentry>
+
+ <refentry id="ST_Transform">
+ <refnamediv>
+ <refname>ST_Transform</refname>
+
+ <refpurpose>Return a new geometry with coordinates transformed to
+ a different spatial reference system.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <funcdef>geometry <function>ST_Transform</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>g1</parameter></paramdef>
+ <paramdef><type>integer </type> <parameter>srid</parameter></paramdef>
+ </funcprototype>
<funcprototype>
<funcdef>geometry <function>ST_Transform</function></funcdef>
@@ -147,55 +147,55 @@ SRID=3785;POINT(-13732990.8753491 6178458.96425423)
<paramdef><type>integer </type> <parameter>to_srid</parameter></paramdef>
</funcprototype>
- </funcsynopsis>
- </refsynopsisdiv>
+ </funcsynopsis>
+ </refsynopsisdiv>
- <refsection>
- <title>Description</title>
+ <refsection>
+ <title>Description</title>
<para>Returns a new geometry with its coordinates transformed to
a different spatial reference system. The destination spatial
- reference <varname>to_srid</varname> may be identified by a valid
- SRID integer parameter (i.e. it must exist in the
- <varname>spatial_ref_sys</varname> table).
- Alternatively, a spatial reference defined as a PROJ.4 string
- can be used for <varname>to_proj</varname> and/or
- <varname>from_proj</varname>, however these methods are not
- optimized. If the destination spatial reference system is
- expressed with a PROJ.4 string instead of an SRID, the SRID of the
- output geometry will be set to zero. With the exception of functions with
- <varname>from_proj</varname>, input geometries must have a defined SRID.
- </para>
-
- <para>ST_Transform is often confused with <xref linkend="ST_SetSRID" />. ST_Transform actually changes the coordinates
- of a geometry from one spatial reference system to another, while ST_SetSRID() simply changes the SRID identifier of
- the geometry.</para>
-
- <note>
- <para>Requires PostGIS be compiled with PROJ support. Use <xref linkend="PostGIS_Full_Version" /> to confirm you have PROJ support compiled in.</para>
- </note>
-
- <note>
- <para>If using more than one transformation, it is useful to have a functional index on the commonly used
- transformations to take advantage of index usage.</para>
- </note>
-
- <note><para>Prior to 1.3.4, this function crashes if used with geometries that contain CURVES. This is fixed in 1.3.4+</para></note>
-
- <para>Enhanced: 2.0.0 support for Polyhedral surfaces was introduced.</para>
- <para>Enhanced: 2.3.0 support for direct PROJ.4 text was introduced.</para>
- <para>&sqlmm_compliant; SQL-MM 3: 5.1.6</para>
- <para>&curve_support;</para>
- <para>&P_support;</para>
-
- </refsection>
-
- <refsection>
- <title>Examples</title>
- <para>Change Massachusetts state plane US feet geometry to WGS 84 long lat</para>
- <programlisting>
+ reference <varname>to_srid</varname> may be identified by a valid
+ SRID integer parameter (i.e. it must exist in the
+ <varname>spatial_ref_sys</varname> table).
+ Alternatively, a spatial reference defined as a PROJ.4 string
+ can be used for <varname>to_proj</varname> and/or
+ <varname>from_proj</varname>, however these methods are not
+ optimized. If the destination spatial reference system is
+ expressed with a PROJ.4 string instead of an SRID, the SRID of the
+ output geometry will be set to zero. With the exception of functions with
+ <varname>from_proj</varname>, input geometries must have a defined SRID.
+ </para>
+
+ <para>ST_Transform is often confused with <xref linkend="ST_SetSRID" />. ST_Transform actually changes the coordinates
+ of a geometry from one spatial reference system to another, while ST_SetSRID() simply changes the SRID identifier of
+ the geometry.</para>
+
+ <note>
+ <para>Requires PostGIS be compiled with PROJ support. Use <xref linkend="PostGIS_Full_Version" /> to confirm you have PROJ support compiled in.</para>
+ </note>
+
+ <note>
+ <para>If using more than one transformation, it is useful to have a functional index on the commonly used
+ transformations to take advantage of index usage.</para>
+ </note>
+
+ <note><para>Prior to 1.3.4, this function crashes if used with geometries that contain CURVES. This is fixed in 1.3.4+</para></note>
+
+ <para>Enhanced: 2.0.0 support for Polyhedral surfaces was introduced.</para>
+ <para>Enhanced: 2.3.0 support for direct PROJ.4 text was introduced.</para>
+ <para>&sqlmm_compliant; SQL-MM 3: 5.1.6</para>
+ <para>&curve_support;</para>
+ <para>&P_support;</para>
+
+ </refsection>
+
+ <refsection>
+ <title>Examples</title>
+ <para>Change Massachusetts state plane US feet geometry to WGS 84 long lat</para>
+ <programlisting>
SELECT ST_AsText(ST_Transform(ST_GeomFromText('POLYGON((743238 2967416,743238 2967450,
- 743265 2967450,743265.625 2967416,743238 2967416))',2249),4326)) As wgs_geom;
+ 743265 2967450,743265.625 2967416,743238 2967416))',2249),4326)) As wgs_geom;
wgs_geom
---------------------------
@@ -207,25 +207,25 @@ SELECT ST_AsText(ST_Transform(ST_GeomFromText('POLYGON((743238 2967416,743238 29
--3D Circular String example
SELECT ST_AsEWKT(ST_Transform(ST_GeomFromEWKT('SRID=2249;CIRCULARSTRING(743238 2967416 1,743238 2967450 2,743265 2967450 3,743265.625 2967416 3,743238 2967416 4)'),4326));
- st_asewkt
+ st_asewkt
--------------------------------------------------------------------------------------
SRID=4326;CIRCULARSTRING(-71.1776848522251 42.3902896512902 1,-71.1776843766326 42.3903829478009 2,
-71.1775844305465 42.3903826677917 3,
-71.1775825927231 42.3902893647987 3,-71.1776848522251 42.3902896512902 4)
- </programlisting>
- <para>Example of creating a partial functional index. For tables where you are not sure all the geometries
- will be filled in, its best to use a partial index that leaves out null geometries which will both conserve space and make your index smaller and more efficient.</para>
- <programlisting>
+ </programlisting>
+ <para>Example of creating a partial functional index. For tables where you are not sure all the geometries
+ will be filled in, its best to use a partial index that leaves out null geometries which will both conserve space and make your index smaller and more efficient.</para>
+ <programlisting>
CREATE INDEX idx_geom_26986_parcels
ON parcels
USING gist
(ST_Transform(geom, 26986))
WHERE geom IS NOT NULL;
- </programlisting>
+ </programlisting>
- <para>Examples of using PROJ.4 text to transform with custom spatial references.</para>
- <programlisting>
+ <para>Examples of using PROJ.4 text to transform with custom spatial references.</para>
+ <programlisting>
-- Find intersection of two polygons near the North pole, using a custom Gnomic projection
-- See http://boundlessgeo.com/2012/02/flattening-the-peel/
WITH data AS (
@@ -242,36 +242,36 @@ CREATE INDEX idx_geom_26986_parcels
st_astext
--------------------------------------------------------------------------------
POLYGON((-170 74.053793645338,-141 73.4268621378904,-141 68,-170 68,-170 74.053793645338))
- </programlisting>
-
- </refsection>
- <refsection>
- <title>Configuring transformation behavior</title>
- <para>Sometimes coordinate transformation involving a grid-shift
- can fail, for example if PROJ.4 has not been built with
- grid-shift files or the coordinate does not lie within the
- range for which the grid shift is defined. By default, PostGIS
- will throw an error if a grid shift file is not present, but
- this behavior can be configured on a per-SRID basis either
- by testing different <varname>to_proj</varname> values of
- PROJ.4 text, or altering the <varname>proj4text</varname> value
- within the <varname>spatial_ref_sys</varname> table.
- </para>
- <para>For example, the proj4text parameter +datum=NAD87 is a shorthand form for the following +nadgrids parameter:</para>
- <programlisting>+nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat</programlisting>
- <para>The @ prefix means no error is reported if the files are not present, but if the end of the list is reached with no file having been appropriate (ie. found and overlapping) then an error is issued.</para>
- <para>If, conversely, you wanted to ensure that at least the standard files were present, but that if all files were scanned without a hit a null transformation is applied you could use:</para>
- <programlisting>+nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat,null</programlisting>
- <para>The null grid shift file is a valid grid shift file covering the whole world and applying no shift. So for a complete example, if you wanted to alter PostGIS so that transformations to SRID 4267 that didn't lie within the correct range did not throw an ERROR, you would use the following:</para>
- <programlisting>UPDATE spatial_ref_sys SET proj4text = '+proj=longlat +ellps=clrk66 +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat,null +no_defs' WHERE srid = 4267;</programlisting>
- </refsection>
-
- <!-- Optionally add a "See Also" section -->
- <refsection>
- <title>See Also</title>
-
- <para><xref linkend="spatial_ref_sys" />, <xref linkend="ST_SetSRID" />, <xref linkend="ST_SRID" />, <xref linkend="UpdateGeometrySRID"/></para>
- </refsection>
- </refentry>
+ </programlisting>
+
+ </refsection>
+ <refsection>
+ <title>Configuring transformation behavior</title>
+ <para>Sometimes coordinate transformation involving a grid-shift
+ can fail, for example if PROJ.4 has not been built with
+ grid-shift files or the coordinate does not lie within the
+ range for which the grid shift is defined. By default, PostGIS
+ will throw an error if a grid shift file is not present, but
+ this behavior can be configured on a per-SRID basis either
+ by testing different <varname>to_proj</varname> values of
+ PROJ.4 text, or altering the <varname>proj4text</varname> value
+ within the <varname>spatial_ref_sys</varname> table.
+ </para>
+ <para>For example, the proj4text parameter +datum=NAD87 is a shorthand form for the following +nadgrids parameter:</para>
+ <programlisting>+nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat</programlisting>
+ <para>The @ prefix means no error is reported if the files are not present, but if the end of the list is reached with no file having been appropriate (ie. found and overlapping) then an error is issued.</para>
+ <para>If, conversely, you wanted to ensure that at least the standard files were present, but that if all files were scanned without a hit a null transformation is applied you could use:</para>
+ <programlisting>+nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat,null</programlisting>
+ <para>The null grid shift file is a valid grid shift file covering the whole world and applying no shift. So for a complete example, if you wanted to alter PostGIS so that transformations to SRID 4267 that didn't lie within the correct range did not throw an ERROR, you would use the following:</para>
+ <programlisting>UPDATE spatial_ref_sys SET proj4text = '+proj=longlat +ellps=clrk66 +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat,null +no_defs' WHERE srid = 4267;</programlisting>
+ </refsection>
+
+ <!-- Optionally add a "See Also" section -->
+ <refsection>
+ <title>See Also</title>
+
+ <para><xref linkend="spatial_ref_sys" />, <xref linkend="ST_SetSRID" />, <xref linkend="ST_SRID" />, <xref linkend="UpdateGeometrySRID"/></para>
+ </refsection>
+ </refentry>
</sect1>
commit 1bc9ace3593bd3711cc388682b9891f68652418e
Author: Robert Coup <robert at coup.net.nz>
Date: Sun Aug 28 13:53:24 2022 +0200
regress: add tests for transform pipelines
diff --git a/regress/core/regress_proj_pipeline.sql b/regress/core/regress_proj_pipeline.sql
new file mode 100644
index 000000000..2b5b98998
--- /dev/null
+++ b/regress/core/regress_proj_pipeline.sql
@@ -0,0 +1,23 @@
+-- Tests for ST_TransformPipeline()
+SELECT '1', ST_AsText(ST_SnapToGrid(ST_TransformPipeline('SRID=4326;POINT(174 -37)'::geometry, '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80'), 10));
+
+-- also checks deg->rad on input
+SELECT '2', ST_AsText(ST_SnapToGrid(ST_TransformPipeline('SRID=4326;POINT(2 49)'::geometry, 'urn:ogc:def:coordinateOperation:EPSG::16031'), 10));
+
+-- error: invalid pipeline definition
+SELECT '3', ST_AsText(ST_TransformPipeline('SRID=4326;POINT(0 0)'::geometry, 'bad coordinate transform pipeline'));
+
+-- error: code defines a CRS, not a coordinateOperation
+SELECT '4', ST_AsText(ST_TransformPipeline('SRID=4326;POINT(0 0)'::geometry, 'EPSG:2193'));
+
+-- assigning/checking SRID on results
+SELECT '5', ST_SRID(ST_TransformPipeline('SRID=4326;POINT(174 -37)'::geometry, '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80'));
+SELECT '6', ST_SRID(ST_TransformPipeline('SRID=4326;POINT(174 -37)'::geometry, '+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80', 12345));
+
+-- GDA2020 is complicated, check it works
+SELECT '7', ST_AsText(ST_SnapToGrid(ST_TransformPipeline('SRID=4283;POINT(151.2 -33.8)'::geometry, 'urn:ogc:def:coordinateOperation:EPSG::8048'), 0.1));
+
+-- testing inverse pipelines; also checks rad->deg on results
+SELECT '8', ST_AsText(ST_SnapToGrid(ST_InverseTransformPipeline('SRID=32631;POINT(426857 5427937)'::geometry, 'urn:ogc:def:coordinateOperation:EPSG::16031'), 0.1));
+SELECT '9', ST_SRID(ST_InverseTransformPipeline('SRID=32631;POINT(426857 5427937)'::geometry, 'urn:ogc:def:coordinateOperation:EPSG::16031', 12345));
+
diff --git a/regress/core/regress_proj_pipeline_expected b/regress/core/regress_proj_pipeline_expected
new file mode 100644
index 000000000..43834cbd0
--- /dev/null
+++ b/regress/core/regress_proj_pipeline_expected
@@ -0,0 +1,9 @@
+1|POINT(1688980 5904660)
+2|POINT(426860 5427940)
+ERROR: could not parse coordinate operation 'bad coordinate transform pipeline'
+ERROR: lwgeom_transform_pipeline: Failed to transform
+5|0
+6|12345
+7|POINT(151.20000000000002 -33.800000000000004)
+8|POINT(2 49)
+9|12345
diff --git a/regress/core/tests.mk.in b/regress/core/tests.mk.in
index 1ceae8573..60bd6d23c 100644
--- a/regress/core/tests.mk.in
+++ b/regress/core/tests.mk.in
@@ -87,6 +87,7 @@ TESTS += \
$(top_srcdir)/regress/core/regress_proj_adhoc \
$(top_srcdir)/regress/core/regress_proj_cache_overflow \
$(top_srcdir)/regress/core/regress_proj_4890 \
+ $(top_srcdir)/regress/core/regress_proj_pipeline \
$(top_srcdir)/regress/core/relate \
$(top_srcdir)/regress/core/remove_repeated_points \
$(top_srcdir)/regress/core/removepoint \
commit d33ebde11f679cbc100aecf6d98c608d211820c6
Author: Robert Coup <robert at coup.net.nz>
Date: Sun Aug 28 13:41:28 2022 +0200
postgis: Add transform pipeline SQL functions.
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.
diff --git a/postgis/lwgeom_transform.c b/postgis/lwgeom_transform.c
index 1f59b306e..646f0e6dd 100644
--- a/postgis/lwgeom_transform.c
+++ b/postgis/lwgeom_transform.c
@@ -34,6 +34,7 @@
Datum transform(PG_FUNCTION_ARGS);
Datum transform_geom(PG_FUNCTION_ARGS);
+Datum transform_pipeline_geom(PG_FUNCTION_ARGS);
Datum postgis_proj_version(PG_FUNCTION_ARGS);
Datum LWGEOM_asKML(PG_FUNCTION_ARGS);
@@ -149,6 +150,51 @@ Datum transform_geom(PG_FUNCTION_ARGS)
}
+/**
+ * transform_pipeline_geom( GEOMETRY, TEXT (input proj pipeline),
+ * BOOL (is forward), INT (output srid)
+ */
+PG_FUNCTION_INFO_V1(transform_pipeline_geom);
+Datum transform_pipeline_geom(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *gser, *gser_result=NULL;
+ LWGEOM *geom;
+ char *input_pipeline;
+ bool is_forward;
+ int32 result_srid;
+ int rv;
+
+ /* Take a copy, since we will be altering the coordinates */
+ gser = PG_GETARG_GSERIALIZED_P_COPY(0);
+
+ /* Convert from text to cstring for libproj */
+ input_pipeline = text_to_cstring(PG_GETARG_TEXT_P(1));
+ is_forward = PG_GETARG_BOOL(2);
+ result_srid = PG_GETARG_INT32(3);
+
+ geom = lwgeom_from_gserialized(gser);
+ rv = lwgeom_transform_pipeline(geom, input_pipeline, is_forward);
+ pfree(input_pipeline);
+
+ if (rv == LW_FAILURE)
+ {
+ elog(ERROR, "coordinate transformation failed");
+ PG_RETURN_NULL();
+ }
+
+ /* Re-compute bbox if input had one (COMPUTE_BBOX TAINTING) */
+ geom->srid = result_srid;
+ if (geom->bbox)
+ lwgeom_refresh_bbox(geom);
+
+ gser_result = geometry_serialize(geom);
+ lwgeom_free(geom);
+ PG_FREE_IF_COPY(gser, 0);
+
+ PG_RETURN_POINTER(gser_result); /* new geometry */
+}
+
+
PG_FUNCTION_INFO_V1(postgis_proj_version);
Datum postgis_proj_version(PG_FUNCTION_ARGS)
{
diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
index 2e50a77f5..dc5e93fd9 100644
--- a/postgis/postgis.sql.in
+++ b/postgis/postgis.sql.in
@@ -2876,6 +2876,30 @@ CREATE OR REPLACE FUNCTION ST_Transform(geom geometry, from_proj text, to_srid i
LANGUAGE 'sql' IMMUTABLE STRICT PARALLEL SAFE
_COST_HIGH;
+CREATE OR REPLACE FUNCTION postgis_transform_pipeline_geometry(geom geometry, pipeline text, forward boolean, to_srid integer)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME','transform_pipeline_geom'
+ LANGUAGE 'c' IMMUTABLE STRICT PARALLEL SAFE
+ _COST_MEDIUM;
+
+-- Availability: 3.4.0
+
+CREATE OR REPLACE FUNCTION ST_TransformPipeline(geom geometry, pipeline text, to_srid integer DEFAULT 0)
+ RETURNS geometry AS
+ 'SELECT @extschema at .postgis_transform_pipeline_geometry($1, $2, TRUE, $3)'
+ LANGUAGE 'sql' IMMUTABLE STRICT PARALLEL SAFE
+ _COST_HIGH;
+
+-- Availability: 3.4.0
+
+CREATE OR REPLACE FUNCTION ST_InverseTransformPipeline(geom geometry, pipeline text, to_srid integer DEFAULT 0)
+ RETURNS geometry AS
+ 'SELECT @extschema at .postgis_transform_pipeline_geometry($1, $2, FALSE, $3)'
+ LANGUAGE 'sql' IMMUTABLE STRICT PARALLEL SAFE
+ _COST_HIGH;
+
+-- Availability: 3.4.0
+
-----------------------------------------------------------------------
-- POSTGIS_VERSION()
-----------------------------------------------------------------------
commit 5606caad8fde789045d6b7cdaca5340cb14ccc10
Author: Robert Coup <robert at coup.net.nz>
Date: Sun Aug 28 13:32:33 2022 +0200
liblwgeom: implement transforms via pipeline
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.
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index 8a8b25349..aad9398d7 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -2453,6 +2453,32 @@ int ptarray_transform(POINTARRAY *pa, LWPROJ* pj);
*/
LWPROJ *lwproj_from_str(const char* str_in, const char* str_out);
+/**
+ * Transform (reproject) a geometry in-place using a PROJ pipeline.
+ * @param geom the geometry to transform
+ * @param pipeline the coordinate operation pipeline string. Either:
+ * - `urn:ogc:def:coordinateOperation:AUTHORITY::CODE`
+ * - a PROJ pipeline string: `+proj=pipeline ...`
+ * - concatenated operations: `urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618`
+ * @param is_forward whether to execute the pipeline in a forward or inverse
+ * direction.
+ */
+int lwgeom_transform_pipeline(LWGEOM *geom, const char* pipeline, bool is_forward);
+
+/**
+ * Allocate a new LWPROJ containing the reference to the PROJ's PJ using a
+ * PROJ pipeline definition:
+ * @param str_pipeline the coordinate operation pipeline string. Either:
+ * - `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 coordinate operation accepted via `proj_create()`
+ * @param is_foward whether to execute the the pipeline in a forward or inverse
+ * direction.
+ */
+LWPROJ *lwproj_from_str_pipeline(const char* str_pipeline, bool is_forward);
+
+
/*******************************************************************************
* GEOS-dependent extra functions on LWGEOM
******************************************************************************/
diff --git a/liblwgeom/lwgeom_transform.c b/liblwgeom/lwgeom_transform.c
index 19461aad8..64e2a3ecc 100644
--- a/liblwgeom/lwgeom_transform.c
+++ b/liblwgeom/lwgeom_transform.c
@@ -118,6 +118,41 @@ lwproj_from_str(const char* str_in, const char* str_out)
return lp;
}
+LWPROJ *
+lwproj_from_str_pipeline(const char* str_pipeline, bool is_forward)
+{
+ /* Usable inputs? */
+ if (!str_pipeline)
+ return NULL;
+
+ PJ* pj = proj_create(PJ_DEFAULT_CTX, str_pipeline);
+ if (!pj)
+ return NULL;
+
+ /* check we have a transform, not a crs */
+ if (proj_is_crs(pj))
+ return NULL;
+
+ /* Add in an axis swap if necessary */
+ PJ* pj_norm = proj_normalize_for_visualization(PJ_DEFAULT_CTX, pj);
+ if (!pj_norm)
+ pj_norm = pj;
+ /* Swap is not a copy of input? Clean up input */
+ else if (pj != pj_norm)
+ proj_destroy(pj);
+
+ /* Allocate and populate return value */
+ LWPROJ *lp = lwalloc(sizeof(LWPROJ));
+ lp->pj = pj_norm; /* Caller is going to have to explicitly proj_destroy this */
+ lp->pipeline_is_forward = is_forward;
+
+ /* this is stuff for geography calculations; doesn't matter here */
+ lp->source_is_latlong = LW_FALSE;
+ lp->source_semi_major_metre = DBL_MAX;
+ lp->source_semi_minor_metre = DBL_MAX;
+ return lp;
+}
+
int
lwgeom_transform_from_str(LWGEOM *geom, const char* instr, const char* outstr)
{
@@ -148,6 +183,28 @@ lwgeom_transform_from_str(LWGEOM *geom, const char* instr, const char* outstr)
return ret;
}
+int
+lwgeom_transform_pipeline(LWGEOM *geom, const char* pipelinestr, bool is_forward)
+{
+ LWPROJ *lp = lwproj_from_str_pipeline(pipelinestr, is_forward);
+ if (!lp)
+ {
+ PJ *pj_in = proj_create(PJ_DEFAULT_CTX, pipelinestr);
+ if (!pj_in)
+ {
+ proj_errno_reset(NULL);
+ lwerror("could not parse coordinate operation '%s'", pipelinestr);
+ }
+ proj_destroy(pj_in);
+ lwerror("%s: Failed to transform", __func__);
+ return LW_FAILURE;
+ }
+ int ret = lwgeom_transform(geom, lp);
+ proj_destroy(lp->pj);
+ lwfree(lp);
+ return ret;
+}
+
int
ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
{
commit 3d6c8256034aec3dcb91d7c36a4bbe86ff2697ef
Author: Robert Coup <robert at coup.net.nz>
Date: Sun Aug 28 10:55:59 2022 +0200
LWPROJ: add transform direction
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.
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index e696d56de..8a8b25349 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -32,6 +32,7 @@
#include <stdarg.h>
#include <stdint.h>
#include <stdio.h>
+#include <stdbool.h>
#include "../postgis_config.h"
@@ -43,6 +44,10 @@
typedef struct LWPROJ
{
PJ* pj;
+
+ /* for pipeline transforms, whether to do a forward or inverse */
+ bool pipeline_is_forward;
+
/* Source crs is geographic: Used in geography calls (source srid == dst srid) */
uint8_t source_is_latlong;
/* Source ellipsoid parameters */
diff --git a/liblwgeom/lwgeom_transform.c b/liblwgeom/lwgeom_transform.c
index 7e7d6735f..19461aad8 100644
--- a/liblwgeom/lwgeom_transform.c
+++ b/liblwgeom/lwgeom_transform.c
@@ -111,6 +111,7 @@ lwproj_from_str(const char* str_in, const char* str_out)
/* Allocate and populate return value */
LWPROJ *lp = lwalloc(sizeof(LWPROJ));
lp->pj = pj_norm; /* Caller is going to have to explicitly proj_destroy this */
+ lp->pipeline_is_forward = true;
lp->source_is_latlong = source_is_latlong;
lp->source_semi_major_metre = semi_major_metre;
lp->source_semi_minor_metre = semi_minor_metre;
@@ -158,8 +159,10 @@ ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
int has_z = ptarray_has_z(pa);
double *pa_double = (double*)(pa->serialized_pointlist);
+ PJ_DIRECTION direction = pj->pipeline_is_forward ? PJ_FWD : PJ_INV;
+
/* Convert to radians if necessary */
- if (proj_angular_input(pj->pj, PJ_FWD))
+ if (proj_angular_input(pj->pj, direction))
{
for (i = 0; i < pa->npoints; i++)
{
@@ -175,7 +178,7 @@ ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
PJ_XYZT v = {pa_double[0], pa_double[1], has_z ? pa_double[2] : 0.0, 0.0};
PJ_COORD c;
c.xyzt = v;
- PJ_COORD t = proj_trans(pj->pj, PJ_FWD, c);
+ PJ_COORD t = proj_trans(pj->pj, direction, c);
int pj_errno_val = proj_errno_reset(pj->pj);
if (pj_errno_val)
@@ -199,7 +202,7 @@ ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
*/
n_converted = proj_trans_generic(pj->pj,
- PJ_FWD,
+ direction,
pa_double,
point_size,
n_points, /* X */
@@ -229,7 +232,7 @@ ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
}
/* Convert radians to degrees if necessary */
- if (proj_angular_output(pj->pj, PJ_FWD))
+ if (proj_angular_output(pj->pj, direction))
{
for (i = 0; i < pa->npoints; i++)
{
commit 3cefc95100088f98c69b7f5984d119728fb88f9e
Author: Robert Coup <robert at coup.net.nz>
Date: Sun Aug 28 10:49:53 2022 +0200
lwgeom: fix angular unit conversion in transform
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.
diff --git a/liblwgeom/lwgeom_transform.c b/liblwgeom/lwgeom_transform.c
index d761293bf..7e7d6735f 100644
--- a/liblwgeom/lwgeom_transform.c
+++ b/liblwgeom/lwgeom_transform.c
@@ -165,6 +165,7 @@ ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
{
getPoint4d_p(pa, i, &p);
to_rad(&p);
+ ptarray_set_point4d(pa, i, &p);
}
}
@@ -234,6 +235,7 @@ ptarray_transform(POINTARRAY *pa, LWPROJ *pj)
{
getPoint4d_p(pa, i, &p);
to_dec(&p);
+ ptarray_set_point4d(pa, i, &p);
}
}
-----------------------------------------------------------------------
Summary of changes:
NEWS | 5 +-
doc/introduction.xml | 1 +
doc/reference_output.xml | 26 +-
doc/reference_srs.xml | 541 ++++++++++++++++++----------
doc/using_postgis_dataman.xml | 2 +-
liblwgeom/liblwgeom.h.in | 31 ++
liblwgeom/lwgeom_transform.c | 70 +++-
postgis/lwgeom_transform.c | 46 +++
postgis/postgis.sql.in | 24 ++
regress/core/regress_proj_pipeline.sql | 23 ++
regress/core/regress_proj_pipeline_expected | 9 +
regress/core/tests.mk.in | 1 +
12 files changed, 560 insertions(+), 219 deletions(-)
create mode 100644 regress/core/regress_proj_pipeline.sql
create mode 100644 regress/core/regress_proj_pipeline_expected
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list