[postgis-tickets] r14698 - #2991, Enable ST_Transform to use PROJ.4 text (Mike Toews)
Daniel Baston
dbaston at gmail.com
Thu Feb 25 10:03:26 PST 2016
Author: dbaston
Date: 2016-02-25 10:03:25 -0800 (Thu, 25 Feb 2016)
New Revision: 14698
Modified:
trunk/NEWS
trunk/doc/reference_editor.xml
trunk/postgis/lwgeom_transform.c
trunk/postgis/postgis.sql.in
trunk/regress/regress_proj.sql
trunk/regress/regress_proj_expected
Log:
#2991, Enable ST_Transform to use PROJ.4 text (Mike Toews)
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2016-02-25 16:02:29 UTC (rev 14697)
+++ trunk/NEWS 2016-02-25 18:03:25 UTC (rev 14698)
@@ -10,6 +10,7 @@
- TopoGeom_addElement, TopoGeom_remElement (Sandro Santilli)
- populate_topology_layer (Sandro Santilli)
- #2259 ST_Voronoi (Dan Baston)
+ - #2991 Enable ST_Transform to use PROJ.4 text (Mike Toews)
- #3339 ST_GeneratePoints (Paul Ramsey)
- #3362 ST_ClusterDBSCAN (Dan Baston)
- #3428 ST_Points (Dan Baston)
Modified: trunk/doc/reference_editor.xml
===================================================================
--- trunk/doc/reference_editor.xml 2016-02-25 16:02:29 UTC (rev 14697)
+++ trunk/doc/reference_editor.xml 2016-02-25 18:03:25 UTC (rev 14698)
@@ -1855,7 +1855,7 @@
<refname>ST_Transform</refname>
<refpurpose>Returns a new geometry with its coordinates transformed to
- the SRID referenced by the integer parameter.</refpurpose>
+ a different spatial reference.</refpurpose>
</refnamediv>
<refsynopsisdiv>
@@ -1865,18 +1865,50 @@
<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>
+ <paramdef><type>geometry </type> <parameter>geom</parameter></paramdef>
+ <paramdef><type>text </type> <parameter>to_proj</parameter></paramdef>
+ </funcprototype>
+
+ <funcprototype>
+ <funcdef>geometry <function>ST_Transform</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>geom</parameter></paramdef>
+ <paramdef><type>text </type> <parameter>from_proj</parameter></paramdef>
+ <paramdef><type>text </type> <parameter>to_proj</parameter></paramdef>
+ </funcprototype>
+
+ <funcprototype>
+ <funcdef>geometry <function>ST_Transform</function></funcdef>
+ <paramdef><type>geometry </type> <parameter>geom</parameter></paramdef>
+ <paramdef><type>text </type> <parameter>from_proj</parameter></paramdef>
+ <paramdef><type>integer </type> <parameter>to_srid</parameter></paramdef>
+ </funcprototype>
+
</funcsynopsis>
</refsynopsisdiv>
<refsection>
<title>Description</title>
- <para>Returns a new geometry with its coordinates transformed to
- spatial reference system referenced by the SRID integer parameter. The destination SRID
- must exist in the <varname>SPATIAL_REF_SYS</varname> table.</para>
+ <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 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>
+ 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>
@@ -1890,6 +1922,7 @@
<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>
@@ -1898,7 +1931,7 @@
<refsection>
<title>Examples</title>
- <para>Change Mass state plane US feet geometry to WGS 84 long lat</para>
+ <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;
@@ -1929,11 +1962,40 @@
(ST_Transform(the_geom, 26986))
WHERE the_geom IS NOT NULL;
</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 (
+ SELECT
+ ST_GeomFromText('POLYGON((170 50,170 72,-130 72,-130 50,170 50))', 4326) AS p1,
+ ST_GeomFromText('POLYGON((-170 68,-170 90,-141 90,-141 68,-170 68))', 4326) AS p2,
+ '+proj=gnom +ellps=WGS84 +lat_0=70 +lon_0=-160 +no_defs'::text AS gnom
+ )
+ SELECT ST_AsText(
+ ST_Transform(
+ ST_Intersection(ST_Transform(p1, gnom), ST_Transform(p2, gnom)),
+ gnom, 4326))
+ FROM data;
+ st_astext
+ --------------------------------------------------------------------------------
+ POLYGON((-170 74.053793645338,-141 73.4268621378904,-141 68,-170 68,-170 74.053793645338))
+ </programlisting>
+
</refsection>
-
<refsection>
<title>Configuring transformation behaviour</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 behaviour can be configured on a per-SRID basis by altering the proj4text value within the spatial_ref_sys table.</para>
+ <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 behaviour 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>
Modified: trunk/postgis/lwgeom_transform.c
===================================================================
--- trunk/postgis/lwgeom_transform.c 2016-02-25 16:02:29 UTC (rev 14697)
+++ trunk/postgis/lwgeom_transform.c 2016-02-25 18:03:25 UTC (rev 14698)
@@ -122,22 +122,8 @@
int32 result_srid ;
char *pj_errstr;
-
-
result_srid = PG_GETARG_INT32(3);
- if (result_srid == SRID_UNKNOWN)
- {
- elog(ERROR,"transform: destination SRID = %d",SRID_UNKNOWN);
- PG_RETURN_NULL();
- }
-
geom = PG_GETARG_GSERIALIZED_P_COPY(0);
- if (gserialized_get_srid(geom) == SRID_UNKNOWN)
- {
- pfree(geom);
- elog(ERROR,"transform_geom: source SRID = %d",SRID_UNKNOWN);
- PG_RETURN_NULL();
- }
/* Set the search path if we haven't already */
SetPROJ4LibPath();
Modified: trunk/postgis/postgis.sql.in
===================================================================
--- trunk/postgis/postgis.sql.in 2016-02-25 16:02:29 UTC (rev 14697)
+++ trunk/postgis/postgis.sql.in 2016-02-25 18:03:25 UTC (rev 14698)
@@ -2661,7 +2661,6 @@
$$
LANGUAGE 'plpgsql' IMMUTABLE STRICT;
-
---------------------------------------------------------------
-- PROJ support
---------------------------------------------------------------
@@ -2674,6 +2673,17 @@
$$
LANGUAGE 'plpgsql' IMMUTABLE STRICT;
+-- Availability: 1.2.2
+CREATE OR REPLACE FUNCTION ST_SetSRID(geometry,int4)
+ RETURNS geometry
+ AS 'MODULE_PATHNAME','LWGEOM_set_srid'
+ LANGUAGE 'c' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION ST_SRID(geometry)
+ RETURNS int4
+ AS 'MODULE_PATHNAME','LWGEOM_get_srid'
+ LANGUAGE 'c' IMMUTABLE STRICT;
+
CREATE OR REPLACE FUNCTION postgis_transform_geometry(geometry,text,text,int)
RETURNS geometry
AS 'MODULE_PATHNAME','transform_geom'
@@ -2685,7 +2695,26 @@
AS 'MODULE_PATHNAME','transform'
LANGUAGE 'c' IMMUTABLE STRICT;
+-- Availability: 2.3.0
+CREATE OR REPLACE FUNCTION ST_Transform(geom geometry, to_proj text)
+ RETURNS geometry AS
+'SELECT postgis_transform_geometry($1, proj4text, $2, 0)
+FROM spatial_ref_sys WHERE srid=ST_SRID($1);'
+ LANGUAGE sql IMMUTABLE STRICT;
+-- Availability: 2.3.0
+CREATE OR REPLACE FUNCTION ST_Transform(geom geometry, from_proj text, to_proj text)
+ RETURNS geometry AS
+'SELECT postgis_transform_geometry($1, $2, $3, 0)'
+ LANGUAGE sql IMMUTABLE STRICT;
+
+-- Availability: 2.3.0
+CREATE OR REPLACE FUNCTION ST_Transform(geom geometry, from_proj text, to_srid integer)
+ RETURNS geometry AS
+'SELECT postgis_transform_geometry($1, $2, proj4text, $3)
+FROM spatial_ref_sys WHERE srid=$3;'
+ LANGUAGE sql IMMUTABLE STRICT;
+
-----------------------------------------------------------------------
-- POSTGIS_VERSION()
-----------------------------------------------------------------------
@@ -4340,18 +4369,7 @@
AS 'MODULE_PATHNAME', 'LWGEOM_isempty'
LANGUAGE 'c' IMMUTABLE STRICT;
-CREATE OR REPLACE FUNCTION ST_SRID(geometry)
- RETURNS int4
- AS 'MODULE_PATHNAME','LWGEOM_get_srid'
- LANGUAGE 'c' IMMUTABLE STRICT;
-
-- Availability: 1.2.2
-CREATE OR REPLACE FUNCTION ST_SetSRID(geometry,int4)
- RETURNS geometry
- AS 'MODULE_PATHNAME','LWGEOM_set_srid'
- LANGUAGE 'c' IMMUTABLE STRICT;
-
--- Availability: 1.2.2
CREATE OR REPLACE FUNCTION ST_AsBinary(geometry,text)
RETURNS bytea
AS 'MODULE_PATHNAME','LWGEOM_asBinary'
Modified: trunk/regress/regress_proj.sql
===================================================================
--- trunk/regress/regress_proj.sql 2016-02-25 16:02:29 UTC (rev 14697)
+++ trunk/regress/regress_proj.sql 2016-02-25 18:03:25 UTC (rev 14698)
@@ -36,5 +36,25 @@
--- test #8: Transforming to same SRID
SELECT 8,ST_AsEWKT(ST_transform(ST_GeomFromEWKT('SRID=100002;POINT(0 0)'),100002));
+SELECT 9, ST_AsEWKT(ST_SnapToGrid(ST_Transform(
+ ST_GeomFromEWKT('SRID=100002;POINT(16 48)'),
+ '+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs '), 10));
+
+--- test #10: Transform from_proj to_proj
+SELECT 10, ST_AsEWKT(ST_SnapToGrid(ST_Transform(
+ ST_GeomFromEWKT('POINT(16 48)'),
+ '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ',
+ '+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs '), 10));
+
+--- test #11: Transform from_proj to_srid
+SELECT 11, ST_AsEWKT(ST_SnapToGrid(ST_Transform(
+ ST_GeomFromEWKT('POINT(16 48)'),
+ '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ', 100001), 10));
+
+--- test #12: Transform with bad to_proj
+SELECT 12, ST_AsEWKT(ST_Transform(
+ ST_GeomFromEWKT('SRID=100002;POINT(16 48)'),
+ 'invalid projection'));
+
DELETE FROM spatial_ref_sys WHERE srid >= 100000;
Modified: trunk/regress/regress_proj_expected
===================================================================
--- trunk/regress/regress_proj_expected 2016-02-25 16:02:29 UTC (rev 14697)
+++ trunk/regress/regress_proj_expected 2016-02-25 18:03:25 UTC (rev 14698)
@@ -7,3 +7,7 @@
6|16.00000000|48.00000000
ERROR: Input geometry has unknown (0) SRID
8|SRID=100002;POINT(0 0)
+9|POINT(574600 5316780)
+10|POINT(574600 5316780)
+11|SRID=100001;POINT(574600 5316780)
+ERROR: transform_geom: couldn't parse proj4 output string: 'invalid projection': projection not named
More information about the postgis-tickets
mailing list