From 57664ef7799bbc03e6208837ec6d5ba555be5e9b Mon Sep 17 00:00:00 2001
From: Sandro Santilli <strk@keybit.net>
Date: Fri, 20 Aug 2004 16:35:38 +0000
Subject: [PATCH] initial skel for transform()

git-svn-id: http://svn.osgeo.org/postgis/trunk@714 b70326c6-7e19-0410-871a-916f4a2858ee
---
 lwgeom/lwgeom_transform.c | 275 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 275 insertions(+)

diff --git a/lwgeom/lwgeom_transform.c b/lwgeom/lwgeom_transform.c
index e69de29bb..f92f2669d 100644
--- a/lwgeom/lwgeom_transform.c
+++ b/lwgeom/lwgeom_transform.c
@@ -0,0 +1,275 @@
+
+#include <math.h>
+#include <float.h>
+#include <string.h>
+#include <stdio.h>
+#include <errno.h>
+
+#include "postgres.h"
+#include "fmgr.h"
+
+#include "lwgeom.h"
+
+// if USE_PROJECTION undefined, we get a do-nothing transform() function
+#ifndef USE_PROJ
+
+PG_FUNCTION_INFO_V1(transform_geom);
+Datum transform_geom(PG_FUNCTION_ARGS)
+{
+
+	elog(ERROR,"PostGIS transform() called, but support not compiled in.  Modify your makefile to add proj support, remake and re-install");
+	PG_RETURN_NULL();
+}
+
+PG_FUNCTION_INFO_V1(postgis_proj_version);
+Datum postgis_proj_version(PG_FUNCTION_ARGS)
+{
+	PG_RETURN_NULL();
+}
+
+#else // defined USE_PROJ 
+
+#include "projects.h"
+
+PJ *make_project(char *str1);
+void to_rad(POINT3D *pts, int num_points);
+void to_dec(POINT3D *pts, int num_points);
+int pj_transform_nodatum(PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset, double *x, double *y, double *z );
+
+//this is *exactly* the same as PROJ.4's pj_transform(), but it doesnt do the
+// datum shift.
+int
+pj_transform_nodatum(PJ *srcdefn, PJ *dstdefn, long point_count,
+int point_offset, double *x, double *y, double *z )
+{
+    long      i;
+    //int       need_datum_shift;
+
+    pj_errno = 0;
+
+    if( point_offset == 0 )
+        point_offset = 1;
+
+    if( !srcdefn->is_latlong )
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            XY         projected_loc;
+            LP	       geodetic_loc;
+
+            projected_loc.u = x[point_offset*i];
+            projected_loc.v = y[point_offset*i];
+
+            geodetic_loc = pj_inv( projected_loc, srcdefn );
+            if( pj_errno != 0 )
+                return pj_errno;
+
+            x[point_offset*i] = geodetic_loc.u;
+            y[point_offset*i] = geodetic_loc.v;
+        }
+    }
+
+    if( !dstdefn->is_latlong )
+    {
+        for( i = 0; i < point_count; i++ )
+        {
+            XY         projected_loc;
+            LP	       geodetic_loc;
+
+            geodetic_loc.u = x[point_offset*i];
+            geodetic_loc.v = y[point_offset*i];
+
+            projected_loc = pj_fwd( geodetic_loc, dstdefn );
+            if( pj_errno != 0 )
+                return pj_errno;
+
+            x[point_offset*i] = projected_loc.u;
+            y[point_offset*i] = projected_loc.v;
+        }
+    }
+
+    return 0;
+}
+
+
+
+// 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 (output proj4), INT (output srid)
+// tmpPts - if there is a nadgrid error (-38), we re-try the transform on a copy of points.  The transformed points
+//          are in an indeterminate state after the -38 error is thrown.
+PG_FUNCTION_INFO_V1(transform_geom);
+Datum transform_geom(PG_FUNCTION_ARGS)
+{
+	LWGEOM *geom;
+
+	LWGEOM *result=NULL;
+	PJ *input_pj,*output_pj;
+
+/*
+	LWGEOM_INSPECTED *geom_inspected;
+	char *o1;
+	int32 *offsets1;
+	int j,type1, i,poly_points;
+	unsigned char gtype;
+
+	LWPOLY *poly;
+	LWLINE *line;
+	LWPOINT *pt;
+	POINT3D *poly_pts;
+
+
+	BOX3D *bbox;
+	POINT3D *tmpPts;
+*/
+	char *input_proj4, *output_proj4;
+
+	text	   *input_proj4_text;
+	text	   *output_proj4_text;
+	int32	   result_srid ;
+
+	elog(ERROR, "Not implemented yet");
+
+	result_srid   = PG_GETARG_INT32(3);
+	if (result_srid   == -1)
+	{
+		elog(ERROR,"tranform: destination SRID = -1");
+		PG_RETURN_NULL();
+	}
+
+	geom = (LWGEOM *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+	if (lwgeom_getSRID(geom) == -1)
+	{
+		PG_FREE_IF_COPY(geom, 0);
+		elog(ERROR,"tranform: source SRID = -1");
+		PG_RETURN_NULL();
+	}
+
+	input_proj4_text  = (PG_GETARG_TEXT_P(1));
+	output_proj4_text = (PG_GETARG_TEXT_P(2));
+
+	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
+
+	//make input and output projection objects
+	input_pj = make_project(input_proj4);
+	if ( (input_pj == NULL) || pj_errno)
+	{
+		pfree(input_proj4); pfree(output_proj4);
+		PG_FREE_IF_COPY(geom, 0);
+		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);
+		PG_FREE_IF_COPY(geom, 0);
+		elog(ERROR,"tranform: couldnt parse proj4 output string");
+		PG_RETURN_NULL();
+	}
+
+	//great, now we have a geometry, and input/output PJ* structs. 
+	//Excellent.
+
+
+	// clean up
+	pj_free(input_pj);
+	pj_free(output_pj);
+	pfree(input_proj4); pfree(output_proj4);
+
+	elog(ERROR, "Not implemented yet");
+
+	PG_RETURN_POINTER(result); // new geometry
+}
+
+PG_FUNCTION_INFO_V1(postgis_proj_version);
+Datum postgis_proj_version(PG_FUNCTION_ARGS)
+{
+	const char *ver = pj_release;
+	text *result;
+	result = (text *) palloc(VARHDRSZ  + strlen(ver));
+	VARATT_SIZEP(result) = VARHDRSZ + strlen(ver) ;
+	memcpy(VARDATA(result), ver, strlen(ver));
+	PG_RETURN_POINTER(result);
+}
+
+#endif
+
-- 
2.40.0