+
+#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
+
projects / postgis / commitdiff