From: Daniel Baston Date: Thu, 25 Feb 2016 18:03:25 +0000 (+0000) Subject: #2991, Enable ST_Transform to use PROJ.4 text (Mike Toews) X-Git-Tag: 2.3.0beta1~224 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=afa23eb02f979c0092b47c6d131dc9e6b4a1f174;p=postgis #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 eb5bb3c92..661084b5a 100644 --- 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 535666054..7449d6342 100644 --- a/doc/reference_editor.xml +++ b/doc/reference_editor.xml @@ -1855,7 +1855,7 @@ LINESTRING(26 125,54 84,101 100) ST_Transform Returns a new geometry with its coordinates transformed to - the SRID referenced by the integer parameter. + a different spatial reference. @@ -1865,18 +1865,50 @@ LINESTRING(26 125,54 84,101 100) geometry g1 integer srid + + + geometry ST_Transform + geometry geom + text to_proj + + + + geometry ST_Transform + geometry geom + text from_proj + text to_proj + + + + geometry ST_Transform + geometry geom + text from_proj + integer to_srid + + Description - 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 SPATIAL_REF_SYS table. + Returns a new geometry with its coordinates transformed to + a different spatial reference system. The destination spatial + reference to_srid may be identified by a valid + SRID integer parameter (i.e. it must exist in the + spatial_ref_sys table). + Alternatively, a spatial reference defined as a PROJ.4 string + can be used for to_proj and/or + from_proj, 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 + from_proj, input geometries must have a defined SRID. + + 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 + the geometry. Requires PostGIS be compiled with Proj support. Use to confirm you have proj support compiled in. @@ -1890,6 +1922,7 @@ LINESTRING(26 125,54 84,101 100) Prior to 1.3.4, this function crashes if used with geometries that contain CURVES. This is fixed in 1.3.4+ Enhanced: 2.0.0 support for Polyhedral surfaces was introduced. + Enhanced: 2.3.0 support for direct PROJ.4 text was introduced. &sqlmm_compliant; SQL-MM 3: 5.1.6 &curve_support; &P_support; @@ -1898,7 +1931,7 @@ LINESTRING(26 125,54 84,101 100) Examples - Change Mass state plane US feet geometry to WGS 84 long lat + Change Massachusetts state plane US feet geometry to WGS 84 long lat 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; - + Examples of using PROJ.4 text to transform with custom spatial references. + +-- 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)) + + + Configuring transformation behaviour - 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. + 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 to_proj values of + PROJ.4 text, or altering the proj4text value + within the spatial_ref_sys table. + For example, the proj4text parameter +datum=NAD87 is a shorthand form for the following +nadgrids parameter: +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat 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. diff --git a/postgis/lwgeom_transform.c b/postgis/lwgeom_transform.c index 8cefbb895..c93dcb439 100644 --- 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 07fe463cf..8af3ad9dd 100644 --- 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 de3b55e08..fb125350f 100644 --- 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 9b17f6eae..b59b2ed27 100644 --- 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