#2991, Enable ST_Transform to use PROJ.4 text (Mike Toews)

git-svn-id: http://svn.osgeo.org/postgis/trunk@14698 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/NEWS b/NEWS
index eb5bb3c920689651d34f6f4beca03f883509ac59..661084b5a10f75929e6ae9afefe98ee5b694afdb 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -10,6 +10,7 @@ PostGIS 2.3.0
   - 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)

diff --git a/doc/reference_editor.xml b/doc/reference_editor.xml
index 5356660545dc833737e384078ea759826600a030..7449d6342af6663d20c63eb21f68d58cac3cd43f 100644 (file)
--- a/doc/reference_editor.xml
+++ b/doc/reference_editor.xml
@@ -1855,7 +1855,7 @@ LINESTRING(26 125,54 84,101 100)
                <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 @@ LINESTRING(26 125,54 84,101 100)
                        <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 @@ LINESTRING(26 125,54 84,101 100)
                <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 @@ LINESTRING(26 125,54 84,101 100)
 
          <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 @@ CREATE INDEX idx_the_geom_26986_parcels
   (ST_Transform(the_geom, 26986))
   WHERE the_geom IS NOT NULL;
                </programlisting>
-         </refsection>
 
+               <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,@alaska,@ntv2_0.gsb,@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>

diff --git a/postgis/lwgeom_transform.c b/postgis/lwgeom_transform.c
index 8cefbb89598e4cc100c4df9c6bc06c25067903a3..c93dcb439e9d47327dff9d0acc6b6ad5a6c5f109 100644 (file)
--- a/postgis/lwgeom_transform.c
+++ b/postgis/lwgeom_transform.c
@@ -122,22 +122,8 @@ Datum transform_geom(PG_FUNCTION_ARGS)
        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();

diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
index 07fe463cfd90d961f0601df55b045f30c7d9d7f1..8af3ad9ddf1726f5320eeba7663e7babadbc5bb6 100644 (file)
--- a/postgis/postgis.sql.in
+++ b/postgis/postgis.sql.in
@@ -2661,7 +2661,6 @@ END;
 $$
 LANGUAGE 'plpgsql' IMMUTABLE STRICT;
 
-
 ---------------------------------------------------------------
 -- PROJ support
 ---------------------------------------------------------------
@@ -2674,6 +2673,17 @@ END;
 $$
 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,6 +2695,25 @@ CREATE OR REPLACE FUNCTION ST_Transform(geometry,integer)
        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,17 +4369,6 @@ CREATE OR REPLACE FUNCTION ST_IsEmpty(geometry)
        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

diff --git a/regress/regress_proj.sql b/regress/regress_proj.sql
index de3b55e08278085d146ccb6d837d1256982c29c9..fb125350f3f31a17182ebb13c67fc30b54cd23b2 100644 (file)
--- a/regress/regress_proj.sql
+++ b/regress/regress_proj.sql
@@ -36,5 +36,25 @@ SELECT ST_transform(ST_GeomFromEWKT('SRID=0;POINT(0 0)'),100002);
 --- 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;
 

diff --git a/regress/regress_proj_expected b/regress/regress_proj_expected
index 9b17f6eaea215e7c7b799c5340aed5bf57613807..b59b2ed27283cb2a723eba63e0de7de159decca0 100644 (file)
--- a/regress/regress_proj_expected
+++ b/regress/regress_proj_expected
@@ -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