Coordinate transformation function, transform() added in this file.

Adds requirement for linking the proj4 library if non-null version of
function is requested.


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

Makefile

@@ -21,17 +21,18 @@ NAME=postgis
 SO_MAJOR_VERSION=0
 SO_MINOR_VERSION=6
 
+
 #override CPPFLAGS := -I$(srcdir) $(CPPFLAGS)
 # Altered for Cynwin
-override CPPFLAGS := -g  -I$(srcdir) $(CPPFLAGS) -DFRONTEND -DSYSCONFDIR='"$(sysconfdir)"'
+override CPPFLAGS := -g  -I$(srcdir) $(CPPFLAGS) -DFRONTEND -DSYSCONFDIR='"$(sysconfdir)"' -DWANT_PROJECTION 
 override DLLLIBS := $(BE_DLLLIBS) $(DLLLIBS)
 
-OBJS=postgis_debug.o postgis_ops.o postgis_fn.o postgis_inout.o postgis_proj.o postgis_chip.o
+OBJS=postgis_debug.o postgis_ops.o postgis_fn.o postgis_inout.o postgis_proj.o postgis_chip.o postgis_transform.o
 
 # Add libraries that libpq depends (or might depend) on into the
 # shared library link.  (The order in which you list them here doesn't
 # matter.)
-SHLIB_LINK = $(filter -L%, $(LDFLAGS))
+SHLIB_LINK = $(filter -L%, $(LDFLAGS)) -lproj
 
 all: all-lib $(NAME).sql $(NAME)_undef.sql loaderdumper
 

postgis.h

@@ -347,6 +347,9 @@ void deparse_hex(unsigned char str, unsigned char *result);
 
 char *geometry_to_text(GEOMETRY  *geometry);
 
+
+
+
 //exposed to psql
 
 Datum box3d_in(PG_FUNCTION_ARGS);
@@ -456,6 +459,10 @@ Datum height_chip(PG_FUNCTION_ARGS);
 Datum datatype_chip(PG_FUNCTION_ARGS);
 Datum compression_chip(PG_FUNCTION_ARGS);
 
+
+Datum transform_geom(PG_FUNCTION_ARGS);
+
+
 //for GIST index
 typedef char* (*BINARY_UNION)(char*, char*, int*);
 typedef float (*SIZE_BOX)(char*);

postgis.sql.in

@@ -24,7 +24,8 @@ create table spatial_ref_sys (
 	 srid integer not null primary key,
 	 auth_name varchar(256), 
 	 auth_srid integer, 
-	 srtext varchar(2048) 
+	 srtext varchar(2048),
+	 proj4text varchar(2048) 
 					);
 
 -- create the metadata table.  spec, section 3.2.2.1
@@ -52,6 +53,35 @@ CREATE FUNCTION find_SRID(varchar,varchar,varchar) returns int4 as
 -- select find_srid('','geometry_test','mygeom');
 
 
+-- given an SRID, find the proj4 definition string
+CREATE FUNCTION get_proj4_from_srid(integer) returns text as
+'SELECT proj4text::text FROM spatial_ref_sys WHERE srid= $1' 
+LANGUAGE 'sql' with (iscachable,isstrict);
+
+-- select get_proj4_from_srid(1);
+
+
+
+--- drop function transform(geometry,integer);
+-- given a geometry and a SRID, convert the geometry to the new SRID
+--     transform(geometry,new_srid)
+CREATE FUNCTION transform(geometry,integer) returns geometry as
+'BEGIN
+  RETURN transform_geometry( $1 , get_proj4_from_srid(SRID( $1 ) ), get_proj4_from_srid( $2 ), $2  );
+ END;
+'
+LANGUAGE 'plpgsql' with (iscachable,isstrict);
+
+-- test:
+--- trans=# select * from spatial_ref_sys ;
+--- srid |   auth_name   | auth_srid | srtext |                                proj4text                                 
+--- ------+---------------+-----------+--------+--------------------------------------------------------------------------
+---    1 | latlong WGS84 |         1 |        | +proj=longlat +datum=WGS84
+---    2 | BC albers     |         2 |        | proj=aea ellps=GRS80 lon_0=-126 lat_0=45 lat_1=50 lat_2=58.5 x_0=1000000
+
+-- select transform( 'SRID=1;POINT(-120.8 50.3)', 2);
+---   -> 'SRID=2;POINT(1370033.37046971 600755.810968684)'
+
 
 --- DropGeometryColumn(<db name>,<table name>,<column name>)
 --- There is no ALTER TABLE DROP COLUMN command in postgresql
@@ -240,6 +270,12 @@ create function WKB_in(opaque)
    AS '@MODULE_FILENAME@','WKB_in'
    LANGUAGE 'c' with (isstrict);
 
+create function transform_geometry(geometry,text,text,int)
+	RETURNS geometry
+   AS '@MODULE_FILENAME@','transform_geom'
+   LANGUAGE 'c' with (isstrict,iscachable);
+
+
 create function WKB_out(opaque)
 	RETURNS opaque
    AS '@MODULE_FILENAME@','WKB_out'

postgis_proj.c

@@ -162,6 +162,8 @@ double	distance_ellipse(double lat1, double long1,
 	double u2,A,B;
 	double dsigma;
 
+	int iterations;  
+
 	L1 = atan((1.0 - sphere->f ) * tan( lat1) );
 	L2 = atan((1.0 - sphere->f ) * tan( lat2) );
 	sinU1 = sin(L1);
@@ -173,6 +175,7 @@ double	distance_ellipse(double lat1, double long1,
 	dl1 = dl;
 	cosdl1 = cos(dl);
 	sindl1 = sin(dl);
+	iterations = 0;
 	do {
 		cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1;
 		sigma = acos(cosSigma);
@@ -183,7 +186,8 @@ double	distance_ellipse(double lat1, double long1,
 		dl1 = dl + dl2;
 		cosdl1 = cos(dl1);
 		sindl1 = sin(dl1);
-	} while (fabs(dl3) > 1.0e-032);
+		iterations++;
+	} while ( (iterations<999) && (fabs(dl3) > 1.0e-032));
 
 	   // compute expansions A and B 
 	u2 = mu2(azimuthEQ,sphere);

postgis_transform.c

@@ -0,0 +1,265 @@
+#include "postgres.h"
+
+
+#include <math.h>
+#include <float.h>
+#include <string.h>
+#include <stdio.h>
+#include <errno.h>
+
+#include "access/gist.h"
+#include "access/itup.h"
+#include "access/rtree.h"
+
+
+#include "fmgr.h"
+
+
+#include "postgis.h"
+#include "utils/elog.h"
+
+#include "projects.h"
+
+
+#define SHOW_DIGS_DOUBLE 15
+#define MAX_DIGS_DOUBLE (SHOW_DIGS_DOUBLE + 6 + 1 + 3 +1)
+
+
+// if WANT_PROJECTION undefined, we get a do-nothing transform() function
+#ifdef WANT_PROJECTION
+
+PJ *make_project(char *str1);
+void to_rad(POINT3D *pts, int num_points);
+void to_dec(POINT3D *pts, int num_points);
+
+
+// convert decimal degress to radians
+void to_rad(POINT3D *pts, int num_points)
+{
+	int t;
+	for(t=0;t<num_points;t++)
+	{
+		pts[t].x *= PI/180.0;
+		pts[t].y *= PI/180.0;
+	}
+}
+
+// convert radians to decimal degress
+void to_dec(POINT3D *pts, int num_points)
+{
+	int t;
+	for(t=0;t<num_points;t++)
+	{
+		pts[t].x *= 180.0/PI;
+		pts[t].y *= 180.0/PI;
+	}
+
+}
+
+	//given a string, make a PJ object
+PJ *make_project(char *str1)
+{
+	int t;
+	char	*params[1024];  //one for each parameter
+	char  *loc;
+	char 	*str;
+	PJ *result;
+	
+
+	if (str1 == NULL)
+		return NULL;
+
+	if (strlen(str1) ==0)
+		return NULL;
+
+	str = palloc(1+strlen(str1) );
+	strcpy(str,str1);
+
+	//first we split the string into a bunch of smaller strings, based on the " " separator
+	
+	params[0] = str; //1st param, we'll null terminate at the " " soon
+
+	loc = str;
+	t =1;
+	while  ((loc != NULL) && (*loc != 0) )
+	{
+		loc = strchr( loc,' ');
+		if (loc != NULL)
+		{
+			*loc = 0; // null terminate
+			params[t] = loc +1;
+			loc++; // next char
+			t++; //next param 
+		}
+	}
+	
+	if (!(result= pj_init ( t , params)))
+	{
+		pfree(str);	
+		return NULL;
+	}
+	pfree(str);	
+	return result;
+}
+
+
+	//tranform_geom( GEOMETRY, TEXT (input proj4), TEXT (input proj4), INT (output srid)
+PG_FUNCTION_INFO_V1(transform_geom);
+Datum transform_geom(PG_FUNCTION_ARGS)
+{
+	GEOMETRY		   *geom = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+	GEOMETRY		   *result;
+	PJ			   *input_pj,*output_pj;
+
+	char				*o1;
+	int32				*offsets1;
+	int				j,type1,i,poly_points;
+
+	POLYGON3D			*poly;
+	LINE3D			*line;
+	POINT3D			*pt,*poly_pts;
+
+	char				*input_proj4, *output_proj4;
+
+	
+	text	   *input_proj4_text  = (PG_GETARG_TEXT_P(1));
+	text	   *output_proj4_text = (PG_GETARG_TEXT_P(2));
+	int32	   result_srid   = PG_GETARG_INT32(3);
+
+	input_proj4 = (char *) palloc(input_proj4_text->vl_len +1-4);
+      memcpy(input_proj4,input_proj4_text->vl_dat, input_proj4_text->vl_len-4);
+	input_proj4[input_proj4_text->vl_len-4] = 0; //null terminate
+
+	output_proj4 = (char *) palloc(output_proj4_text->vl_len +1-4);
+      memcpy(output_proj4,output_proj4_text->vl_dat, output_proj4_text->vl_len-4);
+	output_proj4[output_proj4_text->vl_len-4] = 0; //null terminate
+
+	if (geom->SRID == -1)
+	{
+		pfree(input_proj4); pfree(output_proj4);
+		elog(ERROR,"tranform: source SRID = -1");
+		PG_RETURN_NULL();			// no srid, cannot convert
+	}
+	
+	if (result_srid   == -1)
+	{
+		pfree(input_proj4); pfree(output_proj4);
+		elog(ERROR,"tranform: destination SRID = -1");	
+		PG_RETURN_NULL();			// no srid, cannot convert
+	}
+
+	//make input and output projection objects
+	input_pj = make_project(input_proj4);
+	if ( (input_pj == NULL) || pj_errno) 
+	{
+		pfree(input_proj4); pfree(output_proj4);
+		elog(ERROR,"tranform: couldnt parse proj4 input string");
+		PG_RETURN_NULL();
+	}
+
+	output_pj = make_project(output_proj4);
+	if ((output_pj == NULL)|| pj_errno)
+	{
+		pfree(input_proj4); pfree(output_proj4);
+		pj_free(input_pj);
+		elog(ERROR,"tranform: couldnt parse proj4 output string");
+		
+		PG_RETURN_NULL();
+	}
+	//great, now we have a geometry, and input/output PJ* structs.  Excellent.
+
+		//copy the geometry structure - we're only going to change the points, not the structures
+	result = (GEOMETRY *) palloc (geom->size);
+	memcpy(result,geom, geom->size);
+
+		//handle each sub-geometry
+		offsets1 = (int32 *) ( ((char *) &(result->objType[0] ))+ sizeof(int32) * result->nobjs ) ;
+		for (j=0; j< result->nobjs; j++)		//for each object in geom1
+		{
+			o1 = (char *) result +offsets1[j] ;  
+			type1=  result->objType[j];
+
+			if (type1 == POINTTYPE)	//point
+			{
+				pt = (POINT3D *) o1;
+				if (input_pj->is_latlong)
+					to_rad(pt,1);
+				pj_transform(input_pj,output_pj, 1,3, &pt->x,&pt->y, &pt->z);
+				if (pj_errno )
+				{
+					pfree(input_proj4); pfree(output_proj4);
+					pj_free(input_pj); pj_free(output_pj);
+					elog(ERROR,"tranform: couldnt project point");
+					PG_RETURN_NULL();	
+				}
+				if (output_pj->is_latlong)
+					to_dec(pt,1);
+			}
+			if (type1 == LINETYPE)	//line
+			{
+				line = (LINE3D *) o1;
+				if (input_pj->is_latlong)
+					to_rad(&line->points[0],line->npoints);
+				pj_transform(input_pj,output_pj, line->npoints ,3, 
+							&line->points[0].x,&line->points[0].y, &line->points[0].z);
+				if (pj_errno )
+				{
+					pfree(input_proj4); pfree(output_proj4);
+					pj_free(input_pj); pj_free(output_pj);
+					elog(ERROR,"tranform: couldnt project line");
+					PG_RETURN_NULL();	
+				}
+				if (output_pj->is_latlong)
+					to_dec(&line->points[0],line->npoints);
+			}
+			if (type1 == POLYGONTYPE)	//POLYGON
+			{
+				poly = (POLYGON3D *) o1;
+				poly_points = 0;
+
+				for (i=0; i<poly->nrings;i++)
+				{
+					poly_points  += poly->npoints[i];
+				}
+				poly_pts = (POINT3D *) ( (char *)&(poly->npoints[poly->nrings] )  );
+				if (input_pj->is_latlong)
+					to_rad(poly_pts , poly_points);
+
+				pj_transform(input_pj,output_pj, poly_points,3, 
+							&poly_pts[0].x,&poly_pts[0].y, &poly_pts[0].z);
+				if (pj_errno)
+				{
+					pfree(input_proj4); pfree(output_proj4);
+					pj_free(input_pj); pj_free(output_pj);
+					elog(ERROR,"tranform: couldnt project polygon");
+					PG_RETURN_NULL();	
+				}
+				if (output_pj->is_latlong)
+					to_dec(poly_pts , poly_points);
+			}
+		}
+
+	// clean up
+	pj_free(input_pj);
+	pj_free(output_pj);
+	pfree(input_proj4); pfree(output_proj4);
+
+	result->SRID = result_srid   ;
+	PG_RETURN_POINTER(result); // new geometry
+}
+
+
+#else   
+
+	// return the original geometry
+PG_FUNCTION_INFO_V1(transform_geom);
+Datum transform_geom(PG_FUNCTION_ARGS)
+{
+	GEOMETRY		   *geom1 = (GEOMETRY *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+	GEOMETRY		   *result;
+
+	result = (GEOMETRY *) palloc (geom1->size);
+	memcpy(result,geom1, geom1->size);
+	PG_RETURN_POINTER(result);
+}
+#endif