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
-- 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
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'
+#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
\ No newline at end of file
projects / postgis / commitdiff