--- /dev/null
+Introduction to the Geography Type
+==================================
+
+The geography type provides native support for spatial features represented on "geographic" coordinates (sometimes called "geodetic" coordinates, or "lat/lon", or "lon/lat"). Geographic coordinates are spherical coordinates expressed in angular units (degrees).
+
+The basis for the PostGIS geometry type is a plane. The shortest path between two points on the plane is a straight line. That means calculations on geometries (areas, distances, lengths, intersections, etc) can be calculated using cartesian mathematics and straight line vectors.
+
+The basis for the PostGIS geographic type is a sphere. The shortest path between two points on the sphere is a great circle arc. That means that calculations on geographies (areas, distances, lengths, intersections, etc) must be calculated on the sphere, using more complicated mathematics. For more accurate measurements, the calculations must take the actual spheroidal shape of the world into account, and the mathematics becomes very complicated indeed.
+
+Because the underlying mathematics is much more complicated, there are fewer functions defined for the geography type than for the geometry type. Over time, as new algorithms are added, the capabilities of the geography type will expand.
+
+
+Geography Basics
+================
+
+The geography type only supports the simplest of simple features:
+
+* POINT
+* LINESTRING
+* POLYGON
+* MULTIPOINT
+* MULTILINESTRING
+* MULTIPOLYGON
+* GEOMETRYCOLLECTION
+
+Properly handing curves on spherical surface would be very difficult indeed.
+
+The text representations for geography objects are the same as for geometries:
+
+* POINT(0 0)
+* LINESTRING(0 0, 1 1)
+* POLYGON((0 0, 1 0, 1 1, 1 0, 0 0))
+* ..etc..
+
+The coordinates for geography objects are expected to be in decimal degrees, with longitude as the first ordinate and latitude as the second:
+
+* POINT(-126 45) is a legal geography point
+* POINT(45 -126) is an illegal geography point
+
+The geography implementation currently assumes that all coordinates are relative to the WGS84 spheroid (using an SRID of 4326), and does not allow any other SRID values to be used. Future enhancements may allow multiple spheroids to be supported and transformations between spheroids. An SRID value of 0 (undefined) will be treated as implicitly using WGS84 (4326).
+
+
+Installing the Geography Extension
+==================================
+
+* Pull the source from the geography sandbox.
+
+ svn checkout http://svn.osgeo.org/postgis/spike/pramsey/geodetic postgis-geodetic-svn
+
+* Compile and install PostGIS as usual.
+
+ cd postgis-geodetic-svn
+ ./configure
+ make
+ make install
+
+* Then, create a database and install PostGIS into that database. (Assuming your PostgreSQL is installed in /usr/local/pgsql.)
+
+ createdb geography
+ psql geography < /usr/local/pgsql/share/contrib/postgis.sql
+ psql geography < /usr/local/pgsql/share/contrib/spatial_ref_sys
+
+* Then, install the geography extension.
+
+ psql geography < /usr/local/pgsql/share/contrib/geography.sql
+
+* That's it!
+
+
+Creating a Geography Table
+==========================
+
+You can create a geography-enabled table using the CREATE TABLE statement as follows:
+
+ CREATE TABLE gtable ( id integer primary key, geog geography );
+
+Check out the contents of "gtable":
+
+ \d gtable
+
+ Table "public.gtable"
+ Column | Type | Modifiers
+ --------+-----------+-----------
+ id | integer | not null
+ geog | geography |
+ Indexes:
+ "gtable_pkey" PRIMARY KEY, btree (id)
+
+
+Now, check the "geography_columns" view and see that your table is listed:
+
+ SELECT * FROM geography_columns;
+
+Note that "geography_columns" metadata is a view against the system tables and kept up to date automatically. This is a big improvement on the manually maintained "geometry_columns" table.
+
+Insert some data into the table and then pull it out:
+
+ INSERT INTO gtable1 VALUES (1, 'POINT(0 0)');
+ SELECT id,ST_AsText(geog) FROM gtable2;
+
+Try inserting some invalid geography data and see the database complain:
+
+ INSERT INTO gtable1 values (1, 'POINT(1000 0)');
+
+
+Using Type Restrictions
+=======================
+
+You can restrict the geography types allowed in a column by adding a type restriction when you create the table:
+
+ CREATE TABLE gtable2 ( id integer, geog geography(linestring) );
+
+ \d gtable2
+
+ Table "public.gtable2"
+ Column | Type | Modifiers
+ --------+-----------------------+-----------
+ id | integer |
+ geog | geography(LineString) |
+
+
+Now, if you try to insert a point, the database will complain:
+
+ INSERT INTO gtable2 VALUES (1, 'POINT(0 0)');
+
+You can also add SRID restrictions to a column, though at this point (with only one SRID supported) there is not a lot of utility to the feature:
+
+ CREATE TABLE gtable3 ( id integer, geog geography(polygon, 4326) );
+
+ \d gtable3
+
+ Table "public.gtable3"
+ Column | Type | Modifiers
+ --------+-------------------------+-----------
+ id | integer |
+ geog | geography(Polygon,4326) |
+
+
+
+Using Input/Output Functions
+============================
+
+There are only four intput/output functions at this time supporting the OGC well-known text (WKT) and well-known binary (WKB) formats. Adding further output formats (GML, GeoJSON) should be straight-forward, borrowing code from the geometry implementations.
+
+* ST_AsText(geography) returns text
+* ST_AsBinary(geography) returns bytea
+* ST_GeographyFromText(text) returns geography
+* ST_GeographyFromBinary(bytea) returns geography
+
+You can test that they are bi-directional by stringing them together:
+
+ SELECT ST_AsText(ST_GeographyFromBinary(ST_AsBinary(ST_GeographyFromText('LINESTRING(0 0, 1 1)'))));
+
+
+Casting from Geometry
+=====================
+
+There is currently a simple cast from geometry to geography, which can be useful for importing larger selections of data into geography until the data loader is upgraded to allow direct shape file imports. In the future, it is possible the cast could do coordinate conversions, and other magic, but for now it is a direct conversion -- if your coordinates are out of range it will error out.
+
+ CREATE TABLE geomtable ( id integer, geom geometry );
+ INSERT INTO geomtable VALUES ( 2, 'POINT(0 0)' );
+ CREATE TABLE geogtable AS SELECT id, geom::geography AS geog FROM geomtable;
+ SELECT ST_AsText(geog), id FROM geogtable;
+
+
lwsegmentize.o \
wktparse.tab.o \
lex.yy.o \
- vsprintf.o
+ vsprintf.o \
+ g_box.o \
+ g_coord.o \
+ g_geometry.o \
+ g_ptarray.o \
+ g_serialized.o \
+ g_util.o \
+ lwgeodetic.o
SA_HEADERS = \
liblwgeom.h \
OBJS= cu_tester.o \
cu_algorithm.o \
- cu_measures.o
+ cu_measures.o \
+ cu_libgeom.o
# If we couldn't find the cunit library then display a helpful message
ifeq ($(CUNIT_LDFLAGS),)
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "cu_libgeom.h"
+
+/*
+** Called from test harness to register the tests in this file.
+*/
+CU_pSuite register_libgeom_suite(void)
+{
+ CU_pSuite pSuite;
+ pSuite = CU_add_suite("LibGeom Suite", init_libgeom_suite, clean_libgeom_suite);
+ if (NULL == pSuite)
+ {
+ CU_cleanup_registry();
+ return NULL;
+ }
+
+ if (
+ (NULL == CU_add_test(pSuite, "test_typmod_macros()", test_typmod_macros)) ||
+ (NULL == CU_add_test(pSuite, "test_flags_macros()", test_flags_macros)) ||
+ (NULL == CU_add_test(pSuite, "test_serialized_srid()", test_serialized_srid)) ||
+ (NULL == CU_add_test(pSuite, "test_gserialized_from_lwgeom_size()", test_gserialized_from_lwgeom_size)) ||
+ (NULL == CU_add_test(pSuite, "test_gserialized_from_lwgeom()", test_gserialized_from_lwgeom)) ||
+ (NULL == CU_add_test(pSuite, "test_lwgeom_from_gserialized()", test_lwgeom_from_gserialized)) ||
+ (NULL == CU_add_test(pSuite, "test_geometry_type_from_string()", test_geometry_type_from_string)) ||
+ (NULL == CU_add_test(pSuite, "test_lwgeom_check_geodetic()", test_lwgeom_check_geodetic)) ||
+ (NULL == CU_add_test(pSuite, "test_lwgeom_count_vertices()", test_lwgeom_count_vertices)) ||
+ (NULL == CU_add_test(pSuite, "test_on_gser_lwgeom_count_vertices()", test_on_gser_lwgeom_count_vertices)) ||
+ (NULL == CU_add_test(pSuite, "test_gbox_from_spherical_coordinates()", test_gbox_from_spherical_coordinates)) ||
+ (NULL == CU_add_test(pSuite, "test_gserialized_get_gbox_geocentric()", test_gserialized_get_gbox_geocentric)) ||
+ (NULL == CU_add_test(pSuite, "test_gbox_calculation()", test_gbox_calculation))
+ )
+ {
+ CU_cleanup_registry();
+ return NULL;
+ }
+ return pSuite;
+}
+
+/*
+** The suite initialization function.
+** Create any re-used objects.
+*/
+int init_libgeom_suite(void)
+{
+ return 0;
+}
+
+/*
+** The suite cleanup function.
+** Frees any global objects.
+*/
+int clean_libgeom_suite(void)
+{
+ return 0;
+}
+
+void test_typmod_macros(void)
+{
+ uint32 typmod = 0;
+ int srid = 4326;
+ int type = 6;
+ int z = 1;
+ int rv;
+
+ typmod = TYPMOD_SET_SRID(typmod,srid);
+ rv = TYPMOD_GET_SRID(typmod);
+ CU_ASSERT_EQUAL(rv, srid);
+
+ typmod = TYPMOD_SET_TYPE(typmod,type);
+ rv = TYPMOD_GET_TYPE(typmod);
+ CU_ASSERT_EQUAL(rv,type);
+
+ typmod = TYPMOD_SET_Z(typmod);
+ rv = TYPMOD_GET_Z(typmod);
+ CU_ASSERT_EQUAL(rv,z);
+
+ rv = TYPMOD_GET_M(typmod);
+ CU_ASSERT_EQUAL(rv,0);
+
+}
+
+void test_flags_macros(void)
+{
+ uchar flags = 0;
+
+ CU_ASSERT_EQUAL(0, FLAGS_GET_Z(flags));
+ flags = FLAGS_SET_Z(flags, 1);
+ CU_ASSERT_EQUAL(1, FLAGS_GET_Z(flags));
+ flags = FLAGS_SET_Z(flags, 0);
+ CU_ASSERT_EQUAL(0, FLAGS_GET_Z(flags));
+ CU_ASSERT_EQUAL(0, FLAGS_GET_BBOX(flags));
+
+ CU_ASSERT_EQUAL(0, FLAGS_GET_M(flags));
+ flags = FLAGS_SET_M(flags, 1);
+ CU_ASSERT_EQUAL(1, FLAGS_GET_M(flags));
+
+ CU_ASSERT_EQUAL(0, FLAGS_GET_BBOX(flags));
+ flags = FLAGS_SET_BBOX(flags, 1);
+ CU_ASSERT_EQUAL(1, FLAGS_GET_BBOX(flags));
+
+ CU_ASSERT_EQUAL(0, FLAGS_GET_GEODETIC(flags));
+ flags = FLAGS_SET_GEODETIC(flags, 1);
+ CU_ASSERT_EQUAL(1, FLAGS_GET_GEODETIC(flags));
+
+ flags = 0;
+ flags = gflags(1, 0, 1); /* z=1, m=0, geodetic=1 */
+
+ CU_ASSERT_EQUAL(1, FLAGS_GET_GEODETIC(flags));
+ CU_ASSERT_EQUAL(1, FLAGS_GET_Z(flags));
+ CU_ASSERT_EQUAL(0, FLAGS_GET_M(flags));
+}
+
+void test_serialized_srid(void)
+{
+ GSERIALIZED s;
+ uint32 srid = 4326;
+ uint32 rv;
+
+ gserialized_set_srid(&s, srid);
+ rv = gserialized_get_srid(&s);
+ CU_ASSERT_EQUAL(rv, srid);
+
+ srid = 0;
+ gserialized_set_srid(&s, srid);
+ rv = gserialized_get_srid(&s);
+ CU_ASSERT_EQUAL(rv, srid);
+
+ srid = 1000000;
+ gserialized_set_srid(&s, srid);
+ rv = gserialized_get_srid(&s);
+ CU_ASSERT_EQUAL(rv, srid);
+}
+
+void test_gserialized_from_lwgeom_size(void)
+{
+ LWGEOM *g;
+ size_t size = 0;
+
+ g = lwgeom_from_ewkt("POINT(0 0)", PARSER_CHECK_NONE);
+ size = gserialized_from_lwgeom_size(g, 0);
+ CU_ASSERT_EQUAL( size, 32 );
+ lwgeom_free(g);
+
+ g = lwgeom_from_ewkt("POINT(0 0 0)", PARSER_CHECK_NONE);
+ size = gserialized_from_lwgeom_size(g, 0);
+ CU_ASSERT_EQUAL( size, 40 );
+ lwgeom_free(g);
+
+ g = lwgeom_from_ewkt("MULTIPOINT(0 0 0, 1 1 1)", PARSER_CHECK_NONE);
+ size = gserialized_from_lwgeom_size(g, 0);
+ CU_ASSERT_EQUAL( size, 80 );
+ lwgeom_free(g);
+
+ g = lwgeom_from_ewkt("LINESTRING(0 0, 1 1)", PARSER_CHECK_NONE);
+ size = gserialized_from_lwgeom_size(g, 0);
+ CU_ASSERT_EQUAL( size, 48 );
+ lwgeom_free(g);
+
+ g = lwgeom_from_ewkt("MULTILINESTRING((0 0, 1 1),(0 0, 1 1))", PARSER_CHECK_NONE);
+ size = gserialized_from_lwgeom_size(g, 0);
+ CU_ASSERT_EQUAL( size, 96 );
+ lwgeom_free(g);
+
+ g = lwgeom_from_ewkt("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", PARSER_CHECK_NONE);
+ size = gserialized_from_lwgeom_size(g, 0);
+ CU_ASSERT_EQUAL( size, 104 );
+ lwgeom_free(g);
+
+ g = lwgeom_from_ewkt("POLYGON((-1 -1, -1 2, 2 2, 2 -1, -1 -1), (0 0, 0 1, 1 1, 1 0, 0 0))", PARSER_CHECK_NONE);
+ size = gserialized_from_lwgeom_size(g, 0);
+ CU_ASSERT_EQUAL( size, 184 );
+ lwgeom_free(g);
+
+}
+
+void test_gserialized_from_lwgeom(void)
+{
+ LWGEOM *geom;
+ GSERIALIZED *g;
+ uint32 type;
+ double *inspect; /* To poke right into the blob. */
+
+ geom = lwgeom_from_ewkt("POINT(0 0.2)", PARSER_CHECK_NONE);
+ g = gserialized_from_lwgeom(geom, 1, 0);
+ type = gserialized_get_type(g);
+ CU_ASSERT_EQUAL( type, POINTTYPE );
+ inspect = (double*)g;
+ CU_ASSERT_EQUAL(inspect[3], 0.2);
+ lwgeom_free(geom);
+ lwfree(g);
+
+ geom = lwgeom_from_ewkt("POLYGON((-1 -1, -1 2.5, 2 2, 2 -1, -1 -1), (0 0, 0 1, 1 1, 1 0, 0 0))", PARSER_CHECK_NONE);
+ g = gserialized_from_lwgeom(geom, 1, 0);
+ type = gserialized_get_type(g);
+ CU_ASSERT_EQUAL( type, POLYGONTYPE );
+ inspect = (double*)g;
+ CU_ASSERT_EQUAL(inspect[9], 2.5);
+ lwgeom_free(geom);
+ lwfree(g);
+
+ geom = lwgeom_from_ewkt("MULTILINESTRING((0 0, 1 1),(0 0.1, 1 1))", PARSER_CHECK_NONE);
+ g = gserialized_from_lwgeom(geom, 1, 0);
+ type = gserialized_get_type(g);
+ CU_ASSERT_EQUAL( type, MULTILINETYPE );
+ inspect = (double*)g;
+ CU_ASSERT_EQUAL(inspect[12], 0.1);
+ lwgeom_free(geom);
+ lwfree(g);
+
+}
+
+
+void test_lwgeom_from_gserialized(void)
+{
+ LWGEOM *geom;
+ GSERIALIZED *g;
+ char *in_ewkt;
+ char *out_ewkt;
+ int i = 0;
+
+ char ewkt[][512] = {
+ "POINT(0 0.2)",
+ "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
+ "MULTIPOINT(0.9 0.9,0.9 0.9,0.9 0.9,0.9 0.9,0.9 0.9,0.9 0.9)",
+ "SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
+ "SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
+ "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
+ "SRID=4326;POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
+ "SRID=4326;POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))",
+ "SRID=100000;POLYGON((-1 -1 3,-1 2.5 3,2 2 3,2 -1 3,-1 -1 3),(0 0 3,0 1 3,1 1 3,1 0 3,0 0 3),(-0.5 -0.5 3,-0.5 -0.4 3,-0.4 -0.4 3,-0.4 -0.5 3,-0.5 -0.5 3))",
+ "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
+ "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
+ "MULTICURVE((5 5 1 3,3 5 2 2,3 3 3 1,0 3 1 1),CIRCULARSTRING(0 0 0 0,0.26794 1 3 -2,0.5857864 1.414213 1 2))",
+ "MULTISURFACE(CURVEPOLYGON(CIRCULARSTRING(-2 0,-1 -1,0 0,1 -1,2 0,0 2,-2 0),(-1 0,0 0.5,1 0,0 1,-1 0)),((7 8,10 10,6 14,4 11,7 8)))",
+ };
+
+ for( i = 0; i < 13; i++ )
+ {
+ in_ewkt = ewkt[i];
+ geom = lwgeom_from_ewkt(in_ewkt, PARSER_CHECK_NONE);
+ g = gserialized_from_lwgeom(geom, 0, 0);
+ lwgeom_free(geom);
+ geom = lwgeom_from_gserialized(g);
+ out_ewkt = lwgeom_to_ewkt(geom, PARSER_CHECK_NONE);
+ //printf("\n in = %s\nout = %s\n", in_ewkt, out_ewkt);
+ CU_ASSERT_STRING_EQUAL(in_ewkt, out_ewkt);
+ lwgeom_release(geom);
+ lwfree(g);
+ lwfree(out_ewkt);
+ }
+
+}
+
+void test_geometry_type_from_string(void)
+{
+ int rv;
+ int type = 0, z = 0, m = 0;
+ char *str;
+
+ str = " POINTZ";
+ rv = geometry_type_from_string(str, &type, &z, &m);
+ //printf("\n in type: %s\nout type: %d\n out z: %d\n out m: %d", str, type, z, m);
+ CU_ASSERT_EQUAL(rv, G_SUCCESS);
+ CU_ASSERT_EQUAL(type, POINTTYPE);
+ CU_ASSERT_EQUAL(z, 1);
+ CU_ASSERT_EQUAL(m, 0);
+
+ str = "LINESTRINGM ";
+ rv = geometry_type_from_string(str, &type, &z, &m);
+ //printf("\n in type: %s\nout type: %d\n out z: %d\n out m: %d", str, type, z, m);
+ CU_ASSERT_EQUAL(rv, G_SUCCESS);
+ CU_ASSERT_EQUAL(type, LINETYPE);
+ CU_ASSERT_EQUAL(z, 0);
+ CU_ASSERT_EQUAL(m, 1);
+
+ str = "MULTIPOLYGONZM";
+ rv = geometry_type_from_string(str, &type, &z, &m);
+ //printf("\n in type: %s\nout type: %d\n out z: %d\n out m: %d", str, type, z, m);
+ CU_ASSERT_EQUAL(rv, G_SUCCESS);
+ CU_ASSERT_EQUAL(type, MULTIPOLYGONTYPE);
+ CU_ASSERT_EQUAL(z, 1);
+ CU_ASSERT_EQUAL(m, 1);
+
+ str = " GEOMETRYCOLLECTIONZM ";
+ rv = geometry_type_from_string(str, &type, &z, &m);
+ //printf("\n in type: %s\nout type: %d\n out z: %d\n out m: %d", str, type, z, m);
+ CU_ASSERT_EQUAL(rv, G_SUCCESS);
+ CU_ASSERT_EQUAL(type, COLLECTIONTYPE);
+ CU_ASSERT_EQUAL(z, 1);
+ CU_ASSERT_EQUAL(m, 1);
+
+ str = " GEOMERYCOLLECTIONZM ";
+ rv = geometry_type_from_string(str, &type, &z, &m);
+ //printf("\n in type: %s\nout type: %d\n out z: %d\n out m: %d", str, type, z, m);
+ CU_ASSERT_EQUAL(rv, G_FAILURE);
+
+}
+
+/*
+* Build LWGEOM on top of *aligned* structure so we can use the read-only
+* point access methods on them.
+*/
+static LWGEOM* lwgeom_over_gserialized(char *wkt)
+{
+ LWGEOM *lwg;
+ GSERIALIZED *g;
+
+ lwg = lwgeom_from_ewkt(wkt, PARSER_CHECK_NONE);
+ g = gserialized_from_lwgeom(lwg, 1, 0);
+ lwgeom_free(lwg);
+ return lwgeom_from_gserialized(g);
+}
+
+void test_lwgeom_check_geodetic(void)
+{
+ LWGEOM *geom;
+ int i = 0;
+
+ char ewkt[][512] = {
+ "POINT(0 0.2)",
+ "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
+ "SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
+ "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
+ "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
+ "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
+ "POINT(0 220.2)",
+ "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
+ "SRID=1;MULTILINESTRING((-1 -131,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
+ "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
+ "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,211 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
+ "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1111 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
+ };
+
+ for( i = 0; i < 6; i++ )
+ {
+ geom = lwgeom_over_gserialized(ewkt[i]);
+ CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), G_TRUE);
+ lwgeom_free(geom);
+ }
+
+ for( i = 6; i < 12; i++ )
+ {
+ //char *out_ewkt;
+ geom = lwgeom_over_gserialized(ewkt[i]);
+ CU_ASSERT_EQUAL(lwgeom_check_geodetic(geom), G_FALSE);
+ //out_ewkt = lwgeom_to_ewkt(geom, PARSER_CHECK_NONE);
+ //printf("%s\n", out_ewkt);
+ lwgeom_free(geom);
+ }
+
+}
+
+void test_lwgeom_count_vertices(void)
+{
+ LWGEOM *geom;
+
+ geom = lwgeom_from_ewkt("MULTIPOINT(-1 -1,-1 2.5,2 2,2 -1)", PARSER_CHECK_NONE);
+ CU_ASSERT_EQUAL(lwgeom_count_vertices(geom),4);
+ lwgeom_free(geom);
+
+ geom = lwgeom_from_ewkt("SRID=1;MULTILINESTRING((-1 -131,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))", PARSER_CHECK_NONE);
+ CU_ASSERT_EQUAL(lwgeom_count_vertices(geom),16);
+ lwgeom_free(geom);
+
+ geom = lwgeom_from_ewkt("SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,211 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))", PARSER_CHECK_NONE);
+ CU_ASSERT_EQUAL(lwgeom_count_vertices(geom),30);
+ lwgeom_free(geom);
+
+}
+
+
+void test_on_gser_lwgeom_count_vertices(void)
+{
+ LWGEOM *lwgeom;
+ GSERIALIZED *g_ser1;
+ size_t ret_size;
+
+ lwgeom = lwgeom_from_ewkt("MULTIPOINT(-1 -1,-1 2.5,2 2,2 -1,1 1,2 2,4 5)", PARSER_CHECK_NONE);
+ CU_ASSERT_EQUAL(lwgeom_count_vertices(lwgeom),7);
+ g_ser1 = gserialized_from_lwgeom(lwgeom, 1, &ret_size);
+ g_ser1->flags = FLAGS_SET_GEODETIC(g_ser1->flags, 1);
+ lwgeom_free(lwgeom);
+
+ lwgeom = lwgeom_from_gserialized(g_ser1);
+ CU_ASSERT_EQUAL(lwgeom_count_vertices(lwgeom),7);
+ lwgeom_release(lwgeom);
+
+ lwgeom = lwgeom_from_gserialized(g_ser1);
+
+ CU_ASSERT_EQUAL(lwgeom_count_vertices(lwgeom),7);
+ lwgeom_release(lwgeom);
+
+ lwfree(g_ser1);
+
+}
+
+
+void test_gbox_from_spherical_coordinates(void)
+{
+ double ll[64];
+ GBOX *box = gbox_new(gflags(0, 0, 1));
+ int rv;
+ POINTARRAY *pa;
+
+ ll[0] = 0.0;
+ ll[1] = 0.0;
+ ll[2] = 45.0;
+ ll[3] = 45.0;
+
+ pa = pointArray_construct((uchar*)ll, 0, 0, 2);
+ rv = ptarray_calculate_gbox_geodetic(pa, box);
+ //printf("%s\n", gbox_to_string(box, gflags(0, 0, 1)));
+ CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.707106781187, 0.001);
+ lwfree(pa);
+
+ ll[0] = 45.0;
+ ll[1] = 45.0;
+ ll[2] = 45.5;
+ ll[3] = 45.5;
+ ll[4] = 46.0;
+ ll[5] = 46.0;
+
+ pa = pointArray_construct((uchar*)ll, 0, 0, 3);
+ rv = ptarray_calculate_gbox_geodetic(pa, box);
+ //printf("%s\n", gbox_to_string(box, gflags(0, 0, 1)));
+ CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.694658370459, 0.001);
+ lwfree(pa);
+
+
+ ll[0] = 45.0;
+ ll[1] = 45.0;
+
+ pa = pointArray_construct((uchar*)ll, 0, 0, 1);
+ rv = ptarray_calculate_gbox_geodetic(pa, box);
+ //printf("%s\n", gbox_to_string(box, gflags(0, 0, 1)));
+ CU_ASSERT_DOUBLE_EQUAL(box->xmin, 0.707106781187, 0.001);
+ lwfree(pa);
+
+ lwfree(box);
+
+}
+
+void test_gserialized_get_gbox_geocentric(void)
+{
+ LWGEOM *lwg;
+ GSERIALIZED *g;
+ GBOX *gbox;
+ int i;
+
+ char wkt[][512] = {
+ "POINT(45 45)",
+ "MULTILINESTRING((45 45, 45 45),(45 45, 45 45))",
+ "MULTILINESTRING((45 45, 45 45),(45 45, 45.5 45.5),(45 45, 46 46))",
+ "GEOMETRYCOLLECTION(MULTILINESTRING((45 45, 45 45),(45 45, 45.5 45.5),(45 45, 46 46)))",
+ "MULTIPOLYGON(((45 45, 45 45, 45 45, 45 45),(45 45, 45 45, 45 45, 45 45)))",
+ "LINESTRING(45 45, 45.1 45.1)",
+ "MULTIPOINT(45 45, 45.05 45.05, 45.1 45.1)"
+ };
+
+ double xmin[] = {
+ 0.707106781187,
+ 0.707106781187,
+ 0.694658370459,
+ 0.694658370459,
+ 0.707106781187,
+ 0.705871570679,
+ 0.705871570679
+ };
+
+ for ( i = 0; i < 7; i++ )
+ {
+ //printf("%s\n", wkt[i]);
+ lwg = lwgeom_from_ewkt(wkt[i], PARSER_CHECK_NONE);
+ g = gserialized_from_lwgeom(lwg, 1, 0);
+ g->flags = FLAGS_SET_GEODETIC(g->flags, 1);
+ lwgeom_free(lwg);
+ gbox = gserialized_calculate_gbox_geocentric(g);
+ //printf("%s\n", gbox_to_string(box, g->flags));
+ CU_ASSERT_DOUBLE_EQUAL(gbox->xmin, xmin[i], 0.00000001);
+ lwfree(g);
+ lwfree(gbox);
+ }
+
+}
+
+void test_gbox_calculation(void)
+{
+
+ LWGEOM *geom;
+ int i = 0;
+ GBOX *gbox = gbox_new(gflags(0,0,0));
+ BOX3D *box3d;
+
+ char ewkt[][512] = {
+ "POINT(0 0.2)",
+ "LINESTRING(-1 -1,-1 2.5,2 2,2 -1)",
+ "SRID=1;MULTILINESTRING((-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
+ "POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
+ "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
+ "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
+ "POINT(0 220.2)",
+ "LINESTRING(-1 -1,-1231 2.5,2 2,2 -1)",
+ "SRID=1;MULTILINESTRING((-1 -131,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1),(-1 -1,-1 2.5,2 2,2 -1))",
+ "POLYGON((-1 -1,-1 2.5,2 2,2 -133,-1 -1),(0 0,0 1,1 1,1 0,0 0))",
+ "SRID=4326;MULTIPOLYGON(((-1 -1,-1 2.5,211 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)),((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5)))",
+ "SRID=4326;GEOMETRYCOLLECTION(POINT(0 1),POLYGON((-1 -1,-1111 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0)),MULTIPOLYGON(((-1 -1,-1 2.5,2 2,2 -1,-1 -1),(0 0,0 1,1 1,1 0,0 0),(-0.5 -0.5,-0.5 -0.4,-0.4 -0.4,-0.4 -0.5,-0.5 -0.5))))",
+ };
+
+ for( i = 0; i < 6; i++ )
+ {
+ geom = lwgeom_over_gserialized(ewkt[i]);
+ lwgeom_calculate_gbox(geom, gbox);
+ box3d = lwgeom_compute_box3d(geom);
+ //printf("%g %g\n", gbox->xmin, box3d->xmin);
+ CU_ASSERT_EQUAL(gbox->xmin, box3d->xmin);
+ CU_ASSERT_EQUAL(gbox->xmax, box3d->xmax);
+ CU_ASSERT_EQUAL(gbox->ymin, box3d->ymin);
+ CU_ASSERT_EQUAL(gbox->ymax, box3d->ymax);
+ lwfree(box3d);
+ }
+
+ lwfree(gbox);
+
+
+}
+
+
+
+
+
+
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "CUnit/Basic.h"
+
+#include "libgeom.h"
+
+/***********************************************************************
+** For new geometry library libgeom.h
+*/
+
+/* Set-up / clean-up functions */
+CU_pSuite register_libgeom_suite(void);
+int init_libgeom_suite(void);
+int clean_libgeom_suite(void);
+
+/* Test functions */
+void test_typmod_macros(void);
+void test_flags_macros(void);
+void test_serialized_srid(void);
+void test_gserialized_from_lwgeom_size(void);
+void test_gserialized_from_lwgeom(void);
+void test_lwgeom_from_gserialized(void);
+void test_geometry_type_from_string(void);
+void test_lwgeom_check_geodetic(void);
+void test_lwgeom_count_vertices(void);
+void test_on_gser_lwgeom_count_vertices(void);
+void test_gbox_from_spherical_coordinates(void);
+void test_gserialized_get_gbox_geocentric(void);
+void test_gbox_calculation(void);
#include "cu_algorithm.h"
#include "cu_measures.h"
+#include "cu_libgeom.h"
/*
return CU_get_error();
}
+ /* Add the libgeom suite to the registry */
+ if (NULL == register_libgeom_suite())
+ {
+ CU_cleanup_registry();
+ return CU_get_error();
+ }
+
/* Run all tests using the CUnit Basic interface */
CU_basic_set_mode(CU_BRM_VERBOSE);
CU_basic_run_tests();
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "libgeom.h"
+#include <math.h>
+
+GBOX* gbox_new(uchar flags)
+{
+ GBOX *g = (GBOX*)lwalloc(sizeof(GBOX));
+ g->flags = flags;
+ return g;
+}
+
+int gbox_merge(GBOX *master_box, GBOX *new_box)
+{
+ assert(master_box);
+ assert(new_box);
+
+ if( master_box->flags != new_box->flags )
+ return G_FAILURE;
+
+ if( new_box->xmin < master_box->xmin) master_box->xmin = new_box->xmin;
+ if( new_box->ymin < master_box->ymin) master_box->ymin = new_box->ymin;
+ if( new_box->xmax > master_box->xmax) master_box->xmax = new_box->xmax;
+ if( new_box->ymax > master_box->ymax) master_box->ymax = new_box->ymax;
+
+ if( FLAGS_GET_Z(master_box->flags) || FLAGS_GET_GEODETIC(master_box->flags) )
+ {
+ if( new_box->zmin < master_box->zmin) master_box->zmin = new_box->zmin;
+ if( new_box->zmax > master_box->zmax) master_box->zmax = new_box->zmax;
+ }
+ if( FLAGS_GET_M(master_box->flags) )
+ {
+ if( new_box->mmin < master_box->mmin) master_box->mmin = new_box->mmin;
+ if( new_box->mmax > master_box->mmax) master_box->mmax = new_box->mmax;
+ }
+
+ return G_SUCCESS;
+}
+
+int gbox_overlaps(GBOX *g1, GBOX *g2)
+{
+ if( g1->flags != g2->flags )
+ lwerror("gbox_overlaps: geometries have mismatched dimensionality");
+
+ if( g1->xmax < g2->xmin || g1->ymax < g2->ymin ||
+ g1->xmin > g2->xmax || g1->ymin > g2->ymax )
+ return G_FALSE;
+ if( FLAGS_GET_Z(g1->flags) || FLAGS_GET_GEODETIC(g1->flags) )
+ {
+ if( g1->zmax < g2->zmin || g1->zmin > g2->zmax )
+ return G_FALSE;
+ }
+ if( FLAGS_GET_M(g1->flags) )
+ {
+ if( g1->mmax < g2->mmin || g1->mmin > g2->mmax )
+ return G_FALSE;
+ }
+ return G_TRUE;
+}
+
+char* gbox_to_string(GBOX *gbox)
+{
+ static int sz = 128;
+ char *str = NULL;
+
+ if( ! gbox )
+ return strdup("NULL POINTER");
+
+ str = (char*)lwalloc(sz);
+
+ if( FLAGS_GET_GEODETIC(gbox->flags) )
+ {
+ snprintf(str, sz, "GBOX((%.12g,%.12g,%.12g),(%.12g,%.12g,%.12g))", gbox->xmin, gbox->ymin, gbox->zmin, gbox->xmax, gbox->ymax, gbox->zmax);
+ return str;
+ }
+ if( FLAGS_GET_Z(gbox->flags) && FLAGS_GET_M(gbox->flags) )
+ {
+ snprintf(str, sz, "GBOX((%.12g,%.12g,%.12g,%.12g),(%.12g,%.12g,%.12g,%.12g))", gbox->xmin, gbox->ymin, gbox->zmin, gbox->mmin, gbox->xmax, gbox->ymax, gbox->zmax, gbox->mmax);
+ return str;
+ }
+ if( FLAGS_GET_Z(gbox->flags) )
+ {
+ snprintf(str, sz, "GBOX((%.12g,%.12g,%.12g),(%.12g,%.12g,%.12g))", gbox->xmin, gbox->ymin, gbox->zmin, gbox->xmax, gbox->ymax, gbox->zmax);
+ return str;
+ }
+ if( FLAGS_GET_M(gbox->flags) )
+ {
+ snprintf(str, sz, "GBOX((%.12g,%.12g,%.12g),(%.12g,%.12g,%.12g))", gbox->xmin, gbox->ymin, gbox->mmin, gbox->xmax, gbox->ymax, gbox->mmax);
+ return str;
+ }
+ snprintf(str, sz, "GBOX((%.12g,%.12g),(%.12g,%.12g))", gbox->xmin, gbox->ymin, gbox->xmax, gbox->ymax);
+ return str;
+}
+
+GBOX* gbox_copy(GBOX *box)
+{
+ GBOX *copy = (GBOX*)lwalloc(sizeof(GBOX));
+ memcpy(copy, box, sizeof(GBOX));
+ return copy;
+}
+
+void gbox_duplicate(GBOX *original, GBOX *duplicate)
+{
+ assert(original);
+ assert(duplicate);
+
+ if( original->flags != duplicate->flags )
+ lwerror("gbox_duplicate: geometries have inconsistent dimensionality");
+
+ duplicate->xmin = original->xmin;
+ duplicate->ymin = original->ymin;
+ duplicate->xmax = original->xmax;
+ duplicate->ymax = original->ymax;
+
+ if( FLAGS_GET_GEODETIC(original->flags) || FLAGS_GET_Z(original->flags) )
+ {
+ duplicate->zmin = original->zmin;
+ duplicate->zmax = original->zmax;
+ }
+ if( FLAGS_GET_M(original->flags) )
+ {
+ duplicate->mmin = original->mmin;
+ duplicate->mmax = original->mmax;
+ }
+ return;
+}
+
+size_t gbox_serialized_size(uchar flags)
+{
+ if( ! FLAGS_GET_BBOX(flags) ) return 0;
+ if( FLAGS_GET_GEODETIC(flags) )
+ return 6 * sizeof(float);
+ else
+ return 2 * FLAGS_NDIMS(flags) * sizeof(float);
+}
+
+
+/* ********************************************************************************
+** Compute cartesian bounding GBOX boxes from LWGEOM.
+*/
+
+static int lwcircle_calculate_gbox(POINT4D p1, POINT4D p2, POINT4D p3, GBOX *gbox)
+{
+ double x1, x2, y1, y2, z1, z2, m1, m2;
+ double angle, radius, sweep;
+ /* angles from center */
+ double a1, a2, a3;
+ /* angles from center once a1 is rotated to zero */
+ double r2, r3;
+ double xe = 0.0, ye = 0.0;
+ POINT4D *center;
+ int i;
+
+ LWDEBUG(2, "lwcircle_calculate_gbox called.");
+
+ radius = lwcircle_center(&p1, &p2, &p3, ¢er);
+ if (radius < 0.0) return G_FAILURE;
+
+ x1 = MAXFLOAT;
+ x2 = -1 * MAXFLOAT;
+ y1 = MAXFLOAT;
+ y2 = -1 * MAXFLOAT;
+
+ a1 = atan2(p1.y - center->y, p1.x - center->x);
+ a2 = atan2(p2.y - center->y, p2.x - center->x);
+ a3 = atan2(p3.y - center->y, p3.x - center->x);
+
+ /* Rotate a2 and a3 such that a1 = 0 */
+ r2 = a2 - a1;
+ r3 = a3 - a1;
+
+ /*
+ * There are six cases here I'm interested in
+ * Clockwise:
+ * 1. a1-a2 < 180 == r2 < 0 && (r3 > 0 || r3 < r2)
+ * 2. a1-a2 > 180 == r2 > 0 && (r3 > 0 && r3 < r2)
+ * 3. a1-a2 > 180 == r2 > 0 && (r3 > r2 || r3 < 0)
+ * Counter-clockwise:
+ * 4. a1-a2 < 180 == r2 > 0 && (r3 < 0 || r3 > r2)
+ * 5. a1-a2 > 180 == r2 < 0 && (r3 < 0 && r3 > r2)
+ * 6. a1-a2 > 180 == r2 < 0 && (r3 < r2 || r3 > 0)
+ * 3 and 6 are invalid cases where a3 is the midpoint.
+ * BBOX is fundamental, so these cannot error out and will instead
+ * calculate the sweep using a3 as the middle point.
+ */
+
+ /* clockwise 1 */
+ if (FP_LT(r2, 0) && (FP_GT(r3, 0) || FP_LT(r3, r2)))
+ {
+ sweep = (FP_GT(r3, 0)) ? (r3 - 2 * M_PI) : r3;
+ }
+ /* clockwise 2 */
+ else if (FP_GT(r2, 0) && FP_GT(r3, 0) && FP_LT(r3, r2))
+ {
+ sweep = (FP_GT(r3, 0)) ? (r3 - 2 * M_PI) : r3;
+ }
+ /* counter-clockwise 4 */
+ else if (FP_GT(r2, 0) && (FP_LT(r3, 0) || FP_GT(r3, r2)))
+ {
+ sweep = (FP_LT(r3, 0)) ? (r3 + 2 * M_PI) : r3;
+ }
+ /* counter-clockwise 5 */
+ else if (FP_LT(r2, 0) && FP_LT(r3, 0) && FP_GT(r3, r2))
+ {
+ sweep = (FP_LT(r3, 0)) ? (r3 + 2 * M_PI) : r3;
+ }
+ /* clockwise invalid 3 */
+ else if (FP_GT(r2, 0) && (FP_GT(r3, r2) || FP_LT(r3, 0)))
+ {
+ sweep = (FP_GT(r2, 0)) ? (r2 - 2 * M_PI) : r2;
+ }
+ /* clockwise invalid 6 */
+ else
+ {
+ sweep = (FP_LT(r2, 0)) ? (r2 + 2 * M_PI) : r2;
+ }
+
+ LWDEBUGF(3, "a1 %.16f, a2 %.16f, a3 %.16f, sweep %.16f", a1, a2, a3, sweep);
+
+ angle = 0.0;
+ for (i=0; i < 6; i++)
+ {
+ switch (i)
+ {
+ /* right extent */
+ case 0:
+ angle = 0.0;
+ xe = center->x + radius;
+ ye = center->y;
+ break;
+ /* top extent */
+ case 1:
+ angle = M_PI_2;
+ xe = center->x;
+ ye = center->y + radius;
+ break;
+ /* left extent */
+ case 2:
+ angle = M_PI;
+ xe = center->x - radius;
+ ye = center->y;
+ break;
+ /* bottom extent */
+ case 3:
+ angle = -1 * M_PI_2;
+ xe = center->x;
+ ye = center->y - radius;
+ break;
+ /* first point */
+ case 4:
+ angle = a1;
+ xe = p1.x;
+ ye = p1.y;
+ break;
+ /* last point */
+ case 5:
+ angle = a3;
+ xe = p3.x;
+ ye = p3.y;
+ break;
+ }
+ /* determine if the extents are outside the arc */
+ if (i < 4)
+ {
+ if (FP_GT(sweep, 0.0))
+ {
+ if (FP_LT(a3, a1))
+ {
+ if (FP_GT(angle, (a3 + 2 * M_PI)) || FP_LT(angle, a1)) continue;
+ }
+ else
+ {
+ if (FP_GT(angle, a3) || FP_LT(angle, a1)) continue;
+ }
+ }
+ else
+ {
+ if (FP_GT(a3, a1))
+ {
+ if (FP_LT(angle, (a3 - 2 * M_PI)) || FP_GT(angle, a1)) continue;
+ }
+ else
+ {
+ if (FP_LT(angle, a3) || FP_GT(angle, a1)) continue;
+ }
+ }
+ }
+
+ LWDEBUGF(3, "lwcircle_calculate_gbox: potential extreame %d (%.16f, %.16f)", i, xe, ye);
+
+ x1 = (FP_LT(x1, xe)) ? x1 : xe;
+ y1 = (FP_LT(y1, ye)) ? y1 : ye;
+ x2 = (FP_GT(x2, xe)) ? x2 : xe;
+ y2 = (FP_GT(y2, ye)) ? y2 : ye;
+ }
+
+ LWDEBUGF(3, "lwcircle_calculate_gbox: extreames found (%.16f %.16f, %.16f %.16f)", x1, y1, x2, y2);
+
+ z1 = FP_MIN(p1.z, p2.z);
+ z1 = FP_MIN(z1, p3.z);
+ z2 = FP_MAX(p1.z, p2.z);
+ z2 = FP_MAX(z2, p3.z);
+
+ m1 = FP_MIN(p1.m, p2.m);
+ m1 = FP_MIN(m1, p3.m);
+ m2 = FP_MAX(p1.m, p2.m);
+ m2 = FP_MAX(m2, p3.m);
+
+ gbox->xmin = x1;
+ gbox->xmax = x2;
+ gbox->ymin = y1;
+ gbox->ymax = y2;
+
+ if( FLAGS_GET_Z(gbox->flags) )
+ {
+ gbox->zmin = z1;
+ gbox->zmax = z2;
+ }
+ if( FLAGS_GET_M(gbox->flags) )
+ {
+ gbox->mmin = m1;
+ gbox->mmax = m2;
+ }
+
+ return G_SUCCESS;
+}
+
+static int ptarray_calculate_gbox( POINTARRAY *pa, GBOX *gbox )
+{
+ int i;
+ POINT4D p;
+ int has_z = FLAGS_GET_Z(gbox->flags);
+ int has_m = FLAGS_GET_M(gbox->flags);
+
+ if ( ! pa ) return G_FAILURE;
+ if ( pa->npoints < 1 ) return G_FAILURE;
+
+ getPoint4d_p(pa, 0, &p);
+ gbox->xmin = gbox->xmax = p.x;
+ gbox->ymin = gbox->ymax = p.y;
+ if( has_z )
+ gbox->zmin = gbox->zmax = p.z;
+ if( has_m )
+ gbox->mmin = gbox->mmax = p.m;
+
+ for ( i = 1 ; i < pa->npoints; i++ )
+ {
+ getPoint4d_p(pa, i, &p);
+ gbox->xmin = FP_MIN(gbox->xmin, p.x);
+ gbox->xmax = FP_MAX(gbox->xmax, p.x);
+ gbox->ymin = FP_MIN(gbox->ymin, p.y);
+ gbox->ymax = FP_MAX(gbox->ymax, p.y);
+ if( has_z )
+ {
+ gbox->zmin = FP_MIN(gbox->zmin, p.z);
+ gbox->zmax = FP_MAX(gbox->zmax, p.z);
+ }
+ if( has_m )
+ {
+ gbox->mmin = FP_MIN(gbox->mmin, p.m);
+ gbox->mmax = FP_MAX(gbox->mmax, p.m);
+ }
+ }
+ return G_SUCCESS;
+}
+
+static int lwcircstring_calculate_gbox(LWCIRCSTRING *curve, GBOX *gbox)
+{
+ uchar flags = gflags(TYPE_HASZ(curve->type), TYPE_HASM(curve->type), 0);
+ GBOX tmp;
+ POINT4D p1, p2, p3;
+ int i;
+
+ if ( ! curve ) return G_FAILURE;
+ if ( curve->points->npoints < 3 ) return G_FAILURE;
+
+ tmp.flags = flags;
+
+ /* Initialize */
+ gbox->xmin = gbox->ymin = gbox->zmin = gbox->mmin = MAXFLOAT;
+ gbox->xmax = gbox->ymax = gbox->zmax = gbox->mmax = -1 * MAXFLOAT;
+
+ for ( i = 2; i < curve->points->npoints; i += 2 )
+ {
+ getPoint4d_p(curve->points, i-2, &p1);
+ getPoint4d_p(curve->points, i-1, &p2);
+ getPoint4d_p(curve->points, i, &p3);
+
+ if (lwcircle_calculate_gbox(p1, p2, p3, &tmp) == G_FAILURE)
+ continue;
+
+ gbox_merge(gbox, &tmp);
+ }
+
+ return G_SUCCESS;
+}
+
+static int lwpoint_calculate_gbox(LWPOINT *point, GBOX *gbox)
+{
+ if ( ! point ) return G_FAILURE;
+ return ptarray_calculate_gbox( point->point, gbox );
+}
+
+static int lwline_calculate_gbox(LWLINE *line, GBOX *gbox)
+{
+ if ( ! line ) return G_FAILURE;
+ return ptarray_calculate_gbox( line->points, gbox );
+}
+
+static int lwpoly_calculate_gbox(LWPOLY *poly, GBOX *gbox)
+{
+ if ( ! poly ) return G_FAILURE;
+ if ( poly->nrings == 0 ) return G_FAILURE;
+ /* Just need to check outer ring */
+ return ptarray_calculate_gbox( poly->rings[0], gbox );
+}
+
+static int lwcollection_calculate_gbox(LWCOLLECTION *coll, GBOX *gbox)
+{
+ GBOX subbox;
+ int i;
+ int result = G_FAILURE;
+ int first = G_TRUE;
+ assert(coll);
+ if( coll->ngeoms == 0 )
+ return G_FAILURE;
+
+ subbox.flags = gbox->flags;
+
+ for( i = 0; i < coll->ngeoms; i++ )
+ {
+ if( lwgeom_calculate_gbox((LWGEOM*)(coll->geoms[i]), &subbox) == G_FAILURE )
+ {
+ continue;
+ }
+ else
+ {
+ if( first )
+ {
+ gbox_duplicate(gbox, &subbox);
+ first = G_FALSE;
+ }
+ else
+ {
+ gbox_merge(gbox, &subbox);
+ }
+ result = G_SUCCESS;
+ }
+ }
+ return result;
+}
+
+int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
+{
+ if ( ! lwgeom ) return G_FAILURE;
+
+ switch (TYPE_GETTYPE(lwgeom->type))
+ {
+ case POINTTYPE:
+ return lwpoint_calculate_gbox((LWPOINT *)lwgeom, gbox);
+ case LINETYPE:
+ return lwline_calculate_gbox((LWLINE *)lwgeom, gbox);
+ case CIRCSTRINGTYPE:
+ return lwcircstring_calculate_gbox((LWCIRCSTRING *)lwgeom, gbox);
+ case POLYGONTYPE:
+ return lwpoly_calculate_gbox((LWPOLY *)lwgeom, gbox);
+ case COMPOUNDTYPE:
+ case CURVEPOLYTYPE:
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTICURVETYPE:
+ case MULTIPOLYGONTYPE:
+ case MULTISURFACETYPE:
+ case COLLECTIONTYPE:
+ return lwcollection_calculate_gbox((LWCOLLECTION *)lwgeom, gbox);
+ }
+ /* Never get here, please. */
+ lwerror("unsupported type (%d)", TYPE_GETTYPE(lwgeom->type));
+ return G_FAILURE;
+}
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "libgeom.h"
+
+GCOORDINATE* gcoord_new(int ndims)
+{
+ GCOORDINATE *coord = (GCOORDINATE*)lwalloc(sizeof(GCOORDINATE));
+ assert(ndims >= 2);
+ assert(ndims <= 4);
+ if( ! coord )
+ {
+ lwerror("Out of memory!");
+ return NULL;
+ }
+ coord->ordinates = (double*)lwalloc(sizeof(double) * ndims);
+ if( ! coord->ordinates )
+ {
+ lwerror("Out of memory!");
+ return NULL;
+ }
+ /* We'll determine extra dimension is Z or M later. */
+ if( ndims == 3 )
+ {
+ coord->flags = gflags(1, 0, 0);
+ }
+ if( ndims == 4 )
+ {
+ coord->flags = gflags(1, 1, 0);
+ }
+ return coord;
+}
+
+GCOORDINATE* gcoord_new_with_flags(uchar flags)
+{
+ GCOORDINATE *coord = (GCOORDINATE*)lwalloc(sizeof(GCOORDINATE));
+ if( ! coord )
+ {
+ lwerror("Out of memory!");
+ return NULL;
+ }
+ coord->ordinates = (double*)lwalloc(sizeof(double) * FLAGS_NDIMS(flags));
+ if( ! coord->ordinates )
+ {
+ lwerror("Out of memory!");
+ return NULL;
+ }
+ coord->flags = flags;
+ return coord;
+}
+
+GCOORDINATE* gcoord_new_with_flags_and_ordinates(uchar flags, double *ordinates)
+{
+ GCOORDINATE *coord = (GCOORDINATE*)lwalloc(sizeof(GCOORDINATE));
+ assert(ordinates);
+ if( ! coord )
+ {
+ lwerror("Out of memory!");
+ return NULL;
+ }
+ coord->ordinates = ordinates;
+ if( ! coord->ordinates )
+ {
+ lwerror("Out of memory!");
+ return NULL;
+ }
+ coord->flags = FLAGS_SET_READONLY(flags, 1);
+ return coord;
+}
+
+GCOORDINATE* gcoord_copy(GCOORDINATE *coord)
+{
+ GCOORDINATE *copy = NULL;
+
+ assert(coord);
+ assert(coord->ordinates);
+
+ copy = (GCOORDINATE*)lwalloc(sizeof(GCOORDINATE));
+ if( ! copy ) return NULL;
+ copy->flags = FLAGS_SET_READONLY(coord->flags, 1);
+ copy->ordinates = (double*)lwalloc(sizeof(double) * FLAGS_NDIMS(copy->flags));
+ if( ! copy->ordinates ) return NULL;
+ memcpy(copy->ordinates, coord->ordinates, FLAGS_NDIMS(copy->flags) * sizeof(double));
+ return copy;
+}
+
+void gcoord_free(GCOORDINATE *coord)
+{
+ if ( ! coord ) return;
+ if ( ! FLAGS_GET_READONLY(coord->flags) && coord->ordinates )
+ lwfree(coord->ordinates);
+ lwfree(coord);
+}
+
+void gcoord_set_x(GCOORDINATE *coord, double x)
+{
+ assert(coord);
+ *(coord->ordinates) = x;
+}
+
+void gcoord_set_y(GCOORDINATE *coord, double y)
+{
+ assert(coord);
+ *(coord->ordinates + 1) = y;
+}
+
+void gcoord_set_z(GCOORDINATE *coord, double z)
+{
+ assert(coord);
+ assert(FLAGS_GET_Z(coord->flags));
+ *(coord->ordinates + 2) = z;
+}
+
+void gcoord_set_m(GCOORDINATE *coord, double m)
+{
+ assert(coord);
+ assert(FLAGS_GET_M(coord->flags));
+ if(FLAGS_GET_Z(coord->flags))
+ {
+ *(coord->ordinates + 3) = m;
+ }
+ else {
+ *(coord->ordinates + 2) = m;
+ }
+}
+
+void gcoord_set_ordinates(GCOORDINATE *coord, double *ordinates)
+{
+ assert(coord);
+ assert(ordinates);
+ coord->ordinates = ordinates;
+}
+
+void gcoord_set_ordinate(GCOORDINATE *coord, double ordinate, int index)
+{
+ assert(coord);
+ assert(FLAGS_NDIMS(coord->flags)>index);
+ assert(index>=0);
+ *(coord->ordinates + index) = ordinate;
+}
+
+double gcoord_get_x(GCOORDINATE *coord)
+{
+ assert(coord);
+ return *(coord->ordinates);
+}
+
+double gcoord_get_y(GCOORDINATE *coord)
+{
+ assert(coord);
+ return *(coord->ordinates + 1);
+}
+
+double gcoord_get_z(GCOORDINATE *coord)
+{
+ assert(coord);
+ assert(FLAGS_GET_Z(coord->flags));
+ return *(coord->ordinates + 2);
+}
+
+double gcoord_get_m(GCOORDINATE *coord)
+{
+ assert(coord);
+ assert(FLAGS_GET_M(coord->flags));
+ if(FLAGS_GET_Z(coord->flags))
+ {
+ return *(coord->ordinates + 3);
+ }
+ else
+ {
+ return *(coord->ordinates + 2);
+ }
+}
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "libgeom.h"
+
+G_LINESTRING* glinestring_new_from_gptarray(GPTARRAY *ptarray)
+{
+ G_LINESTRING *gline = NULL;
+ assert(ptarray);
+ gline = (G_LINESTRING*)lwalloc(sizeof(G_LINESTRING));
+ gline->flags = ptarray->flags;
+ gline->type = LINETYPE;
+ gline->bbox = NULL;
+ gline->srid = 0;
+ gline->points = ptarray;
+ return gline;
+}
+
+G_LINESTRING* glinestring_new(uchar flags)
+{
+ G_LINESTRING *gline = NULL;
+ gline = (G_LINESTRING*)lwalloc(sizeof(G_LINESTRING));
+ gline->flags = flags;
+ gline->type = LINETYPE;
+ gline->bbox = NULL;
+ gline->srid = 0;
+ gline->points = gptarray_new(gline->flags);
+ return gline;
+}
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "libgeom.h"
+
+GPTARRAY* gptarray_new(uchar flags)
+{
+ GPTARRAY *ptarr = (GPTARRAY*)lwalloc(sizeof(GPTARRAY));
+ ptarr->flags = flags;
+ ptarr->capacity = G_PT_ARRAY_DEFAULT_POINTS * FLAGS_NDIMS(flags) * sizeof(double);
+ ptarr->npoints = 0;
+ ptarr->ordinates = (double*)lwalloc(ptarr->capacity);
+ return ptarr;
+};
+
+GPTARRAY* gptarray_copy(GPTARRAY *ptarray)
+{
+ GPTARRAY *copy = NULL;
+ assert(ptarray);
+ copy = (GPTARRAY*)lwalloc(sizeof(GPTARRAY));
+ copy->flags = ptarray->flags;
+ copy->capacity = ptarray->npoints * FLAGS_NDIMS(ptarray->flags) * sizeof(double);
+ copy->npoints = ptarray->npoints;
+ copy->ordinates = (double*)lwalloc(copy->capacity);
+ if( ! copy->ordinates ) return NULL;
+ memcpy(copy->ordinates, ptarray->ordinates, copy->capacity);
+ return copy;
+};
+
+GPTARRAY* gptarray_new_with_size(uchar flags, int npoints)
+{
+ GPTARRAY *ptarray = (GPTARRAY*)lwalloc(sizeof(GPTARRAY));
+ ptarray->flags = flags;
+ ptarray->capacity = FLAGS_NDIMS(flags) * npoints * sizeof(double);
+ ptarray->npoints = npoints;
+ ptarray->ordinates = (double*)lwalloc(ptarray->capacity);
+ return ptarray;
+};
+
+GPTARRAY* gptarray_new_with_ordinates(uchar flags, int npoints, double *ordinates)
+{
+ assert(ordinates);
+ GPTARRAY *ptarray = (GPTARRAY*)lwalloc(sizeof(GPTARRAY));
+ ptarray->flags = flags;
+ ptarray->capacity = 0; /* We aren't managing this memory, watch out. */
+ ptarray->npoints = npoints;
+ ptarray->ordinates = ordinates;
+ return ptarray;
+};
+
+void gptarray_free(GPTARRAY *ptarray)
+{
+ assert(ptarray);
+ assert(ptarray->ordinates);
+ if( ptarray->capacity > 0 ) /* Only free the ordinates if we are managing them. */
+ lwfree(ptarray->ordinates);
+ lwfree(ptarray);
+};
+
+void gptarray_add_coord(GPTARRAY *ptarray, GCOORDINATE *coord)
+{
+ assert(ptarray);
+ assert(ptarray->flags == coord->flags);
+ if( FLAGS_NDIMS(ptarray->flags) * (ptarray->npoints + 1) * sizeof(double) > ptarray->capacity )
+ {
+ ptarray->capacity *= 2;
+ ptarray->ordinates = lwrealloc(ptarray->ordinates, ptarray->capacity);
+ if ( ! ptarray->ordinates )
+ {
+ lwerror("Out of memory!");
+ return;
+ }
+ }
+ memcpy(ptarray->ordinates + FLAGS_NDIMS(ptarray->flags) * ptarray->npoints * sizeof(double),
+ coord->ordinates,
+ FLAGS_NDIMS(coord->flags) * sizeof(double));
+
+ ptarray->npoints++;
+}
+
+GCOORDINATE* gptarray_get_coord_ro(GPTARRAY *ptarray, int i)
+{
+ assert(ptarray);
+
+ GCOORDINATE *coord = gcoord_new(ptarray->flags);
+
+ coord->ordinates = ptarray->ordinates + FLAGS_NDIMS(ptarray->flags) * i;
+
+ return coord;
+}
+
+GCOORDINATE* gptarray_get_coord_new(GPTARRAY *ptarray, int i)
+{
+ assert(ptarray);
+
+ GCOORDINATE *coord = gcoord_new(ptarray->flags);
+
+ memcpy(coord->ordinates,
+ ptarray->ordinates + FLAGS_NDIMS(ptarray->flags) * i,
+ FLAGS_NDIMS(ptarray->flags) * sizeof(double));
+
+ return coord;
+}
+
+
+void gptarray_set_coord(GPTARRAY *ptarray, int i, GCOORDINATE *coord)
+{
+ int dim = 0;
+ int ndims = 0;
+
+ assert(ptarray);
+ assert(coord->flags == ptarray->flags);
+
+ ndims = FLAGS_NDIMS(ptarray->flags);
+
+ for( dim = 0; dim < ndims; dim++ )
+ {
+ *(ptarray->ordinates + i * ndims + dim) = *(coord->ordinates + dim);
+ }
+
+}
+
+void gptarray_set_x(GPTARRAY *ptarray, int i, double x)
+{
+ assert(ptarray);
+
+ *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 0) = x;
+}
+
+void gptarray_set_y(GPTARRAY *ptarray, int i, double y)
+{
+ assert(ptarray);
+
+ *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 1) = y;
+}
+
+void gptarray_set_z(GPTARRAY *ptarray, int i, double z)
+{
+ assert(ptarray);
+ assert(FLAGS_GET_Z(ptarray->flags));
+
+ *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 2) = z;
+}
+
+void gptarray_set_m(GPTARRAY *ptarray, int i, double m)
+{
+ assert(ptarray);
+ assert(FLAGS_GET_M(ptarray->flags));
+
+ if( FLAGS_GET_Z(ptarray->flags))
+ { /* Four coordinates */
+ *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 3) = m;
+ }
+ else
+ { /* Three coordinates */
+ *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 2) = m;
+ }
+}
+
+double gptarray_get_x(GPTARRAY *ptarray, int i)
+{
+ assert(ptarray);
+ return *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 0);
+}
+
+double gptarray_get_y(GPTARRAY *ptarray, int i)
+{
+ assert(ptarray);
+ return *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 1);
+}
+
+double gptarray_get_z(GPTARRAY *ptarray, int i)
+{
+ assert(ptarray);
+ assert(FLAGS_GET_Z(ptarray->flags));
+ return *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 2);
+}
+
+double gptarray_get_m(GPTARRAY *ptarray, int i)
+{
+ assert(ptarray);
+ assert(FLAGS_GET_M(ptarray->flags));
+
+ if( FLAGS_GET_Z(ptarray->flags))
+ { /* Four coordinates */
+ return *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 3);
+ }
+ else
+ { /* Three coordinates */
+ return *(ptarray->ordinates + i * FLAGS_NDIMS(ptarray->flags) + 2);
+ }
+}
+
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "libgeom.h"
+
+/***********************************************************************
+* GSERIALIZED metadata utility functions.
+*/
+
+uint32 gserialized_get_type(GSERIALIZED *s)
+{
+ uint32 *ptr;
+ assert(s);
+ ptr = (uint32*)(s->data);
+ ptr += (gbox_serialized_size(s->flags) / sizeof(uint32));
+ return *ptr;
+}
+
+uint32 gserialized_get_srid(GSERIALIZED *s)
+{
+ uint32 srid = 0;
+ srid = srid | (s->srid[0] << 16);
+ srid = srid | (s->srid[1] << 8);
+ srid = srid | s->srid[2];
+ return srid;
+}
+
+void gserialized_set_srid(GSERIALIZED *s, uint32 srid)
+{
+ LWDEBUGF(3, "Called with srid = %d", srid);
+ s->srid[0] = (srid & 0x000F0000) >> 16;
+ s->srid[1] = (srid & 0x0000FF00) >> 8;
+ s->srid[2] = (srid & 0x000000FF);
+}
+
+
+/***********************************************************************
+* Calculate the GSERIALIZED size for an LWGEOM.
+*/
+
+/* Private functions */
+
+static size_t gserialized_from_any_size(LWGEOM *geom); /* Local prototype */
+
+static size_t gserialized_from_lwpoint_size(LWPOINT *point)
+{
+ size_t size = 4; /* Type number. */
+
+ assert(point);
+
+ if ( point->bbox )
+ {
+ size += TYPE_NDIMS(point->type) * sizeof(double) * 2;
+ }
+
+ size += 4; /* Number of points (one or zero (empty)). */
+ size += point->point->npoints * TYPE_NDIMS(point->type) * sizeof(double);
+
+ LWDEBUGF(3, "point size = %d", size);
+
+ return size;
+}
+
+static size_t gserialized_from_lwline_size(LWLINE *line)
+{
+ size_t size = 4; /* Type number. */
+
+ assert(line);
+
+ if ( line->bbox )
+ {
+ size += TYPE_NDIMS(line->type) * sizeof(double) * 2;
+ }
+
+ size += 4; /* Number of points (zero => empty). */
+ size += line->points->npoints * TYPE_NDIMS(line->type) * sizeof(double);
+
+ LWDEBUGF(3, "linestring size = %d", size);
+
+ return size;
+}
+
+static size_t gserialized_from_lwpoly_size(LWPOLY *poly)
+{
+ size_t size = 4; /* Type number. */
+ int i = 0;
+
+ assert(poly);
+
+ if ( poly->bbox )
+ {
+ size += TYPE_NDIMS(poly->type) * sizeof(double) * 2;
+ }
+
+ size += 4; /* Number of rings (zero => empty). */
+ if( poly->nrings % 2 )
+ size += 4; /* Padding to double alignment. */
+
+ for( i = 0; i < poly->nrings; i++ )
+ {
+ size += 4; /* Number of points in ring. */
+ size += poly->rings[i]->npoints * TYPE_NDIMS(poly->type) * sizeof(double);
+ }
+
+ LWDEBUGF(3, "polygon size = %d", size);
+
+ return size;
+}
+
+static size_t gserialized_from_lwcircstring_size(LWCIRCSTRING *curve)
+{
+ size_t size = 4; /* Type number. */
+
+ assert(curve);
+
+ if ( curve->bbox )
+ {
+ size += TYPE_NDIMS(curve->type) * sizeof(double) * 2;
+ }
+
+ size += 4; /* Number of points (zero => empty). */
+ size += curve->points->npoints * TYPE_NDIMS(curve->type) * sizeof(double);
+
+ LWDEBUGF(3, "circstring size = %d", size);
+
+ return size;
+}
+
+static size_t gserialized_from_lwcollection_size(LWCOLLECTION *col)
+{
+ size_t size = 4; /* Type number. */
+ int i = 0;
+
+ assert(col);
+
+ if ( col->bbox )
+ {
+ size += TYPE_NDIMS(col->type) * sizeof(double) * 2;
+ }
+
+ size += 4; /* Number of sub-geometries (zero => empty). */
+
+ for( i = 0; i < col->ngeoms; i++ )
+ {
+ size_t subsize = gserialized_from_any_size(col->geoms[i]);
+ size += subsize;
+ LWDEBUGF(3, "lwcollection subgeom(%d) size = %d", i, subsize);
+ }
+
+ LWDEBUGF(3, "lwcollection size = %d", size);
+
+ return size;
+}
+
+static size_t gserialized_from_any_size(LWGEOM *geom)
+{
+ int type = TYPE_GETTYPE(geom->type);
+
+ LWDEBUGF(2, "Input type: %s", lwgeom_typename(type));
+
+ switch (type)
+ {
+ case POINTTYPE:
+ return gserialized_from_lwpoint_size((LWPOINT *)geom);
+ case LINETYPE:
+ return gserialized_from_lwline_size((LWLINE *)geom);
+ case POLYGONTYPE:
+ return gserialized_from_lwpoly_size((LWPOLY *)geom);
+ case CIRCSTRINGTYPE:
+ return gserialized_from_lwcircstring_size((LWCIRCSTRING *)geom);
+ case CURVEPOLYTYPE:
+ case COMPOUNDTYPE:
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTICURVETYPE:
+ case MULTIPOLYGONTYPE:
+ case MULTISURFACETYPE:
+ case COLLECTIONTYPE:
+ return gserialized_from_lwcollection_size((LWCOLLECTION *)geom);
+ default:
+ lwerror("Unknown geometry type: %d", type);
+ return 0;
+ }
+}
+
+/* Public function */
+
+size_t gserialized_from_lwgeom_size(LWGEOM *geom, GBOX *gbox)
+{
+ size_t size = 8; /* Header overhead. */
+ assert(geom);
+ size += gserialized_from_any_size(geom);
+ if( gbox )
+ size += gbox_serialized_size(gbox->flags);
+ LWDEBUGF(3, "g_serialize size = %d", size);
+ return size;
+}
+
+/***********************************************************************
+* Serialize an LWGEOM into GSERIALIZED.
+*/
+
+/* Private functions */
+
+static size_t gserialized_from_lwgeom_any(LWGEOM *geom, uchar *buf);
+
+static size_t gserialized_from_lwpoint(LWPOINT *point, uchar *buf)
+{
+ uchar *loc;
+ int ptsize = pointArray_ptsize(point->point);
+ int type = POINTTYPE;
+
+ assert(point);
+ assert(buf);
+
+ if ( TYPE_GETZM(point->type) != TYPE_GETZM(point->point->dims) )
+ lwerror("Dimensions mismatch in lwpoint");
+
+ LWDEBUGF(2, "lwpoint_to_gserialized(%p, %p) called", point, buf);
+
+ loc = buf;
+
+ /* Write in the type. */
+ memcpy(loc, &type, sizeof(uint32));
+ loc += sizeof(uint32);
+ /* Write in the number of points (0 => empty). */
+ memcpy(loc, &(point->point->npoints), sizeof(uint32));
+ loc += sizeof(uint32);
+
+ /* Copy in the ordinates. */
+ if( point->point->npoints > 0 )
+ {
+ memcpy(loc, getPoint_internal(point->point, 0), ptsize);
+ loc += ptsize;
+ }
+
+ return (size_t)(loc - buf);
+}
+
+static size_t gserialized_from_lwline(LWLINE *line, uchar *buf)
+{
+ uchar *loc;
+ int ptsize;
+ size_t size;
+ int type = LINETYPE;
+
+ assert(line);
+ assert(buf);
+
+ LWDEBUGF(2, "lwline_to_gserialized(%p, %p) called", line, buf);
+
+ if ( TYPE_GETZM(line->type) != TYPE_GETZM(line->points->dims) )
+ lwerror("Dimensions mismatch in lwline");
+
+ ptsize = pointArray_ptsize(line->points);
+
+ loc = buf;
+
+ /* Write in the type. */
+ memcpy(loc, &type, sizeof(uint32));
+ loc += sizeof(uint32);
+
+ /* Write in the npoints. */
+ memcpy(loc, &(line->points->npoints), sizeof(uint32));
+ loc += sizeof(uint32);
+
+ LWDEBUGF(3, "lwline_to_gserialized added npoints (%d)", line->points->npoints);
+
+ /* Copy in the ordinates. */
+ if( line->points->npoints > 0 )
+ {
+ size = line->points->npoints * ptsize;
+ memcpy(loc, getPoint_internal(line->points, 0), size);
+ loc += size;
+ }
+ LWDEBUGF(3, "lwline_to_gserialized copied serialized_pointlist (%d bytes)", ptsize * line->points->npoints);
+
+ return (size_t)(loc - buf);
+}
+
+static size_t gserialized_from_lwpoly(LWPOLY *poly, uchar *buf)
+{
+ int i;
+ uchar *loc;
+ int ptsize;
+ int type = POLYGONTYPE;
+
+ assert(poly);
+ assert(buf);
+
+ LWDEBUG(2, "lwpoly_to_gserialized called");
+
+ ptsize = sizeof(double) * TYPE_NDIMS(poly->type);
+ loc = buf;
+
+ /* Write in the type. */
+ memcpy(loc, &type, sizeof(uint32));
+ loc += sizeof(uint32);
+
+ /* Write in the nrings. */
+ memcpy(loc, &(poly->nrings), sizeof(uint32));
+ loc += sizeof(uint32);
+
+ /* Write in the npoints per ring. */
+ for( i = 0; i < poly->nrings; i++ )
+ {
+ memcpy(loc, &(poly->rings[i]->npoints), sizeof(uint32));
+ loc += sizeof(uint32);
+ }
+
+ /* Add in padding if necessary to remain double aligned. */
+ if( poly->nrings % 2 )
+ {
+ loc += sizeof(uint32);
+ }
+
+ /* Copy in the ordinates. */
+ for( i = 0; i < poly->nrings; i++ )
+ {
+ POINTARRAY *pa = poly->rings[i];
+ size_t pasize;
+
+ if ( TYPE_GETZM(poly->type) != TYPE_GETZM(pa->dims) )
+ lwerror("Dimensions mismatch in lwpoly");
+
+ pasize = pa->npoints * ptsize;
+ memcpy(loc, getPoint_internal(pa, 0), pasize);
+ loc += pasize;
+ }
+ return (size_t)(loc - buf);
+}
+
+static size_t gserialized_from_lwcircstring(LWCIRCSTRING *curve, uchar *buf)
+{
+ uchar *loc;
+ int ptsize;
+ size_t size;
+ int type = CIRCSTRINGTYPE;
+
+ assert(curve);
+ assert(buf);
+
+ if (TYPE_GETZM(curve->type) != TYPE_GETZM(curve->points->dims))
+ lwerror("Dimensions mismatch in lwcircstring");
+
+
+ ptsize = pointArray_ptsize(curve->points);
+ loc = buf;
+
+ /* Write in the type. */
+ memcpy(loc, &type, sizeof(uint32));
+ loc += sizeof(uint32);
+
+ /* Write in the npoints. */
+ memcpy(loc, &curve->points->npoints, sizeof(uint32));
+ loc += sizeof(uint32);
+
+ /* Copy in the ordinates. */
+ if( curve->points->npoints > 0 )
+ {
+ size = curve->points->npoints * ptsize;
+ memcpy(loc, getPoint_internal(curve->points, 0), size);
+ loc += size;
+ }
+
+ return (size_t)(loc - buf);
+}
+
+static size_t gserialized_from_lwcollection(LWCOLLECTION *coll, uchar *buf)
+{
+ size_t subsize = 0;
+ uchar *loc;
+ int i;
+ int type;
+
+ assert(coll);
+ assert(buf);
+
+ type = TYPE_GETTYPE(coll->type);
+ loc = buf;
+
+ /* Write in the type. */
+ memcpy(loc, &type, sizeof(uint32));
+ loc += sizeof(uint32);
+
+ /* Write in the number of subgeoms. */
+ memcpy(loc, &coll->ngeoms, sizeof(uint32));
+ loc += sizeof(uint32);
+
+ /* Serialize subgeoms. */
+ for ( i=0; i<coll->ngeoms; i++ )
+ {
+ if (TYPE_GETZM(coll->type) != TYPE_GETZM(coll->geoms[i]->type))
+ lwerror("Dimensions mismatch in lwcollection");
+ subsize = gserialized_from_lwgeom_any(coll->geoms[i], loc);
+ loc += subsize;
+ }
+
+ return (size_t)(loc - buf);
+}
+
+static size_t gserialized_from_lwgeom_any(LWGEOM *geom, uchar *buf)
+{
+ int type = 0;
+
+ assert(geom);
+ assert(buf);
+
+ type = TYPE_GETTYPE(geom->type);
+
+ LWDEBUGF(2, "Input type (%d) %s", type, lwgeom_typename(type));
+ LWDEBUGF(2, "LWGEOM(%p) uchar(%p)", geom, buf);
+
+ switch (type)
+ {
+ case POINTTYPE:
+ return gserialized_from_lwpoint((LWPOINT *)geom, buf);
+ case LINETYPE:
+ return gserialized_from_lwline((LWLINE *)geom, buf);
+ case POLYGONTYPE:
+ return gserialized_from_lwpoly((LWPOLY *)geom, buf);
+ case CIRCSTRINGTYPE:
+ return gserialized_from_lwcircstring((LWCIRCSTRING *)geom, buf);
+ case CURVEPOLYTYPE:
+ case COMPOUNDTYPE:
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTICURVETYPE:
+ case MULTIPOLYGONTYPE:
+ case MULTISURFACETYPE:
+ case COLLECTIONTYPE:
+ return gserialized_from_lwcollection((LWCOLLECTION *)geom, buf);
+ default:
+ lwerror("Unknown geometry type: %d", geom->type);
+ return 0;
+ }
+ return 0;
+}
+
+static size_t gserialized_from_gbox(GBOX *gbox, uchar *buf)
+{
+ uchar *loc;
+ float xmin, xmax, ymin, ymax;
+
+ assert(buf);
+
+ if ( ! gbox )
+ return (size_t)0;
+
+ loc = buf;
+
+ xmin = nextDown_f(gbox->xmin);
+ memcpy(loc, &xmin, sizeof(float));
+ loc += sizeof(float);
+
+ xmax = nextUp_f(gbox->xmax);
+ memcpy(loc, &xmax, sizeof(float));
+ loc += sizeof(float);
+
+ ymin = nextDown_f(gbox->ymin);
+ memcpy(loc, &ymin, sizeof(float));
+ loc += sizeof(float);
+
+ ymax = nextUp_f(gbox->ymax);
+ memcpy(loc, &ymax, sizeof(float));
+ loc += sizeof(float);
+
+ if( FLAGS_GET_GEODETIC(gbox->flags) || FLAGS_GET_Z(gbox->flags) )
+ {
+ float zmin, zmax;
+ zmin = nextDown_f(gbox->zmin);
+ memcpy(loc, &zmin, sizeof(float));
+ loc += sizeof(float);
+
+ zmax = nextUp_f(gbox->zmax);
+ memcpy(loc, &zmax, sizeof(float));
+ loc += sizeof(float);
+ }
+ if( FLAGS_GET_M(gbox->flags) )
+ {
+ float mmin, mmax;
+ mmin = nextDown_f(gbox->mmin);
+ memcpy(loc, &mmin, sizeof(float));
+ loc += sizeof(float);
+
+ mmax = nextUp_f(gbox->mmax);
+ memcpy(loc, &mmax, sizeof(float));
+ loc += sizeof(float);
+ }
+ return (size_t)(loc - buf);
+}
+
+/* Public function */
+
+GSERIALIZED* gserialized_from_lwgeom(LWGEOM *geom, int is_geodetic, size_t *size)
+{
+ size_t expected_size = 0;
+ size_t return_size = 0;
+ uchar *serialized = NULL;
+ uchar *ptr = NULL;
+ GSERIALIZED *g = NULL;
+ GBOX *gbox = NULL;
+ uchar flags = 0;
+
+ assert(geom);
+
+ flags = gflags(TYPE_HASZ(geom->type), TYPE_HASM(geom->type), is_geodetic);
+
+ /*
+ ** We need room for a bounding box in the serialized form.
+ ** Calculate the box and allocate enough size for it.
+ */
+ if( lwgeom_needs_bbox(geom) )
+ {
+ flags = FLAGS_SET_BBOX(flags, 1);
+ gbox = gbox_new(flags);
+ if( is_geodetic )
+ {
+ /* Geodetic input requires geodetic box calculation. */
+ if( lwgeom_calculate_gbox_geodetic(geom, gbox) == G_FAILURE )
+ {
+ flags = FLAGS_SET_BBOX(flags, 0);
+ lwfree(gbox);
+ gbox = NULL;
+ }
+ }
+ else
+ {
+ /* Cartesian box calculation. */
+ if( lwgeom_calculate_gbox(geom, gbox) == G_FAILURE )
+ {
+ flags = FLAGS_SET_BBOX(flags, 0);
+ lwfree(gbox);
+ gbox = NULL;
+ }
+ }
+ }
+
+ /* Set up the uchar buffer into which we are going to write the serialized geometry. */
+ expected_size = gserialized_from_lwgeom_size(geom, gbox);
+ serialized = lwalloc(expected_size);
+ ptr = serialized;
+
+ /* Move past srid and flags. */
+ ptr += 8;
+
+ /* Write in the serialized form of the gbox, if necessary. */
+ if( FLAGS_GET_BBOX(flags) )
+ ptr += gserialized_from_gbox(gbox, ptr);
+ if( gbox )
+ lwfree(gbox);
+
+ /* Write in the serialized form of the geometry. */
+ ptr += gserialized_from_lwgeom_any(geom, ptr);
+
+ /* Calculate size as returned by data processing functions. */
+ return_size = ptr - serialized;
+
+ if( expected_size != return_size ) /* Uh oh! */
+ {
+ lwerror("Return size (%d) not equal to expected size (%d)!", return_size, expected_size);
+ return NULL;
+ }
+
+ if( size ) /* Return the output size to the caller if necessary. */
+ *size = return_size;
+
+ g = (GSERIALIZED*)serialized;
+
+ /*
+ ** We are aping PgSQL code here, PostGIS code should use
+ ** VARSIZE to set this for real.
+ */
+ g->size = return_size << 2;
+
+ if( geom->SRID == 0 || geom->SRID == (uint32)(-1) ) /* Zero is the no-SRID value now. */
+ gserialized_set_srid(g, 0);
+ else
+ gserialized_set_srid(g, geom->SRID);
+
+ g->flags = flags;
+
+ return g;
+}
+
+/***********************************************************************
+* De-serialize GSERIALIZED into an LWGEOM.
+*/
+
+static LWGEOM* lwgeom_from_gserialized_buffer(uchar *data_ptr, uchar g_flags, size_t *g_size);
+
+static LWPOINT* lwpoint_from_gserialized_buffer(uchar *data_ptr, uchar g_flags, size_t *g_size)
+{
+ static uint32 type = POINTTYPE;
+ uchar *start_ptr = data_ptr;
+ LWPOINT *point;
+ uint32 npoints = 0;
+
+ assert(data_ptr);
+
+ point = (LWPOINT*)lwalloc(sizeof(LWPOINT));
+ point->SRID = -1; /* Default */
+ point->bbox = NULL;
+ point->type = lwgeom_makeType_full(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), 0, type, 0);
+
+ data_ptr += 4; /* Skip past the type. */
+ npoints = lw_get_uint32(data_ptr); /* Zero => empty geometry */
+ data_ptr += 4; /* Skip past the npoints. */
+
+ if( npoints > 0 )
+ point->point = pointArray_construct(data_ptr, FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), 1);
+ else
+ point->point = ptarray_construct(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), 0); /* Empty point */
+
+ data_ptr += npoints * FLAGS_NDIMS(g_flags) * sizeof(double);
+
+ if( g_size )
+ *g_size = data_ptr - start_ptr;
+
+ return point;
+}
+
+static LWLINE* lwline_from_gserialized_buffer(uchar *data_ptr, uchar g_flags, size_t *g_size)
+{
+ static uint32 type = LINETYPE;
+ uchar *start_ptr = data_ptr;
+ LWLINE *line;
+ uint32 npoints = 0;
+
+ assert(data_ptr);
+
+ line = (LWLINE*)lwalloc(sizeof(LWLINE));
+ line->SRID = -1; /* Default */
+ line->bbox = NULL;
+ line->type = lwgeom_makeType_full(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), 0, type, 0);
+
+ data_ptr += 4; /* Skip past the type. */
+ npoints = lw_get_uint32(data_ptr); /* Zero => empty geometry */
+ data_ptr += 4; /* Skip past the npoints. */
+
+ if ( npoints > 0 )
+ line->points = pointArray_construct(data_ptr, FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), npoints);
+ else
+ line->points = ptarray_construct(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), 0); /* Empty linestring */
+
+ data_ptr += FLAGS_NDIMS(g_flags) * npoints * sizeof(double);
+
+ if( g_size )
+ *g_size = data_ptr - start_ptr;
+
+ return line;
+}
+
+static LWPOLY* lwpoly_from_gserialized_buffer(uchar *data_ptr, uchar g_flags, size_t *g_size)
+{
+ static uint32 type = POLYGONTYPE;
+ uchar *start_ptr = data_ptr;
+ LWPOLY *poly;
+ uchar *ordinate_ptr;
+ uint32 nrings = 0;
+ int i = 0;
+
+ assert(data_ptr);
+
+ poly = (LWPOLY*)lwalloc(sizeof(LWPOLY));
+ poly->SRID = -1; /* Default */
+ poly->bbox = NULL;
+ poly->type = lwgeom_makeType_full(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), 0, type, 0);
+
+ data_ptr += 4; /* Skip past the polygontype. */
+ nrings = lw_get_uint32(data_ptr); /* Zero => empty geometry */
+ poly->nrings = nrings;
+ LWDEBUGF(4, "nrings = %d", nrings);
+ data_ptr += 4; /* Skip past the nrings. */
+
+ ordinate_ptr = data_ptr; /* Start the ordinate pointer. */
+ if( nrings > 0)
+ {
+ poly->rings = (POINTARRAY**)lwalloc( sizeof(POINTARRAY*) * nrings );
+ ordinate_ptr += nrings * 4; /* Move past all the npoints values. */
+ if( nrings % 2 ) /* If there is padding, move past that too. */
+ ordinate_ptr += 4;
+ }
+ else /* Empty polygon */
+ {
+ poly->rings = NULL;
+ }
+
+ for ( i = 0; i < nrings; i++ )
+ {
+ uint32 npoints = 0;
+
+ /* Read in the number of points. */
+ npoints = lw_get_uint32(data_ptr);
+ data_ptr += 4;
+
+ /* Make a point array for the ring, and move the ordinate pointer past the ring ordinates. */
+ poly->rings[i] = pointArray_construct(ordinate_ptr, FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), npoints);
+ ordinate_ptr += sizeof(double) * FLAGS_NDIMS(g_flags) * npoints;
+ }
+
+ if( g_size )
+ *g_size = ordinate_ptr - start_ptr;
+
+ return poly;
+}
+
+static LWCIRCSTRING* lwcircstring_from_gserialized_buffer(uchar *data_ptr, uchar g_flags, size_t *g_size)
+{
+ static uint32 type = CIRCSTRINGTYPE;
+ uchar *start_ptr = data_ptr;
+ LWCIRCSTRING *circstring;
+ uint32 npoints = 0;
+
+ assert(data_ptr);
+
+ circstring = (LWCIRCSTRING*)lwalloc(sizeof(LWCIRCSTRING));
+ circstring->SRID = -1; /* Default */
+ circstring->bbox = NULL;
+ circstring->type = lwgeom_makeType_full(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), 0, type, 0);
+
+ data_ptr += 4; /* Skip past the circstringtype. */
+ npoints = lw_get_uint32(data_ptr); /* Zero => empty geometry */
+ data_ptr += 4; /* Skip past the npoints. */
+
+ if ( npoints > 0 )
+ circstring->points = pointArray_construct(data_ptr, FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), npoints);
+ else
+ circstring->points = ptarray_construct(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), 0); /* Empty circularstring */
+
+ data_ptr += FLAGS_NDIMS(g_flags) * npoints * sizeof(double);
+
+ if( g_size )
+ *g_size = data_ptr - start_ptr;
+
+ return circstring;
+}
+
+static int lwcollection_from_gserialized_allowed_types(int collectiontype, int subtype)
+{
+ if( collectiontype == COLLECTIONTYPE )
+ return G_TRUE;
+ if( collectiontype == MULTIPOINTTYPE &&
+ subtype == POINTTYPE )
+ return G_TRUE;
+ if( collectiontype == MULTILINETYPE &&
+ subtype == LINETYPE )
+ return G_TRUE;
+ if( collectiontype == MULTIPOLYGONTYPE &&
+ subtype == POLYGONTYPE )
+ return G_TRUE;
+ if( collectiontype == COMPOUNDTYPE &&
+ (subtype == LINETYPE || subtype == CIRCSTRINGTYPE) )
+ return G_TRUE;
+ if( collectiontype == CURVEPOLYTYPE &&
+ (subtype == CIRCSTRINGTYPE || subtype == LINETYPE || subtype == COMPOUNDTYPE) )
+ return G_TRUE;
+ if( collectiontype == MULTICURVETYPE &&
+ (subtype == CIRCSTRINGTYPE || subtype == LINETYPE || subtype == COMPOUNDTYPE) )
+ return G_TRUE;
+ if( collectiontype == MULTISURFACETYPE &&
+ (subtype == POLYGONTYPE || subtype == CURVEPOLYTYPE) )
+ return G_TRUE;
+
+ /* Must be a bad combination! */
+ return G_FALSE;
+}
+
+static LWCOLLECTION* lwcollection_from_gserialized_buffer(uchar *data_ptr, uchar g_flags, size_t *g_size)
+{
+ uint32 type;
+ uchar *start_ptr = data_ptr;
+ LWCOLLECTION *collection;
+ uint32 ngeoms = 0;
+ int i = 0;
+
+ assert(data_ptr);
+
+ type = lw_get_uint32(data_ptr);
+ data_ptr += 4; /* Skip past the type. */
+
+ collection = (LWCOLLECTION*)lwalloc(sizeof(LWCOLLECTION));
+ collection->SRID = -1; /* Default */
+ collection->bbox = NULL;
+ collection->type = lwgeom_makeType_full(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), 0, type, 0);
+
+ ngeoms = lw_get_uint32(data_ptr);
+ collection->ngeoms = ngeoms; /* Zero => empty geometry */
+ data_ptr += 4; /* Skip past the ngeoms. */
+
+ if( ngeoms > 0 )
+ collection->geoms = lwalloc(sizeof(LWGEOM*) * ngeoms);
+ else
+ collection->geoms = NULL;
+
+ for( i = 0; i < ngeoms; i++ )
+ {
+ uint32 subtype = lw_get_uint32(data_ptr);
+ size_t subsize = 0;
+
+ if( ! lwcollection_from_gserialized_allowed_types(type, subtype) )
+ {
+ lwerror("Invalid subtype (%s) for collection type (%s)", lwgeom_typename(subtype), lwgeom_typename(type));
+ lwfree(collection);
+ return NULL;
+ }
+ collection->geoms[i] = lwgeom_from_gserialized_buffer(data_ptr, g_flags, &subsize);
+ data_ptr += subsize;
+ }
+
+ if( g_size )
+ *g_size = data_ptr - start_ptr;
+
+ return collection;
+}
+
+
+LWGEOM* lwgeom_from_gserialized_buffer(uchar *data_ptr, uchar g_flags, size_t *g_size)
+{
+ uint32 type;
+
+ assert(data_ptr);
+
+ type = lw_get_uint32(data_ptr);
+
+ LWDEBUGF(2, "Got type %d (%s)", type, lwgeom_typename(type));
+
+ switch (type)
+ {
+ case POINTTYPE:
+ return (LWGEOM *)lwpoint_from_gserialized_buffer(data_ptr, g_flags, g_size);
+ case LINETYPE:
+ return (LWGEOM *)lwline_from_gserialized_buffer(data_ptr, g_flags, g_size);
+ case CIRCSTRINGTYPE:
+ return (LWGEOM *)lwcircstring_from_gserialized_buffer(data_ptr, g_flags, g_size);
+ case POLYGONTYPE:
+ return (LWGEOM *)lwpoly_from_gserialized_buffer(data_ptr, g_flags, g_size);
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTIPOLYGONTYPE:
+ case COMPOUNDTYPE:
+ case CURVEPOLYTYPE:
+ case MULTICURVETYPE:
+ case MULTISURFACETYPE:
+ case COLLECTIONTYPE:
+ return (LWGEOM *)lwcollection_from_gserialized_buffer(data_ptr, g_flags, g_size);
+ default:
+ lwerror("Unknown geometry type: %d", type);
+ return NULL;
+ }
+}
+
+
+LWGEOM* lwgeom_from_gserialized(GSERIALIZED *g)
+{
+ uchar g_flags = 0;
+ uchar has_srid = 0;
+ uchar *data_ptr = NULL;
+ uint32 g_srid = 0;
+ uint32 g_type = 0;
+ LWGEOM *lwgeom = NULL;
+ size_t g_size = 0;
+
+ assert(g);
+
+ g_srid = gserialized_get_srid(g);
+ LWDEBUGF(4, "Got srid %d", g_srid);
+ g_flags = g->flags;
+ if( g_srid > 0 )
+ has_srid = 1;
+ g_type = gserialized_get_type(g);
+ LWDEBUGF(4, "Got type %d", g_type);
+
+ data_ptr = (uchar*)g->data;
+ if( FLAGS_GET_BBOX(g_flags) )
+ data_ptr += gbox_serialized_size(g_flags);
+
+ lwgeom = lwgeom_from_gserialized_buffer(data_ptr, g_flags, &g_size);
+
+ if( ! lwgeom ) return NULL; /* Ooops! */
+
+ lwgeom->type = lwgeom_makeType_full(FLAGS_GET_Z(g_flags), FLAGS_GET_M(g_flags), has_srid, g_type, FLAGS_GET_BBOX(g_flags));
+
+ if( FLAGS_GET_BBOX(g_flags) )
+ {
+ double *dptr = (double*)g->data;
+ BOX2DFLOAT4 *bbox = (BOX2DFLOAT4*)lwalloc(sizeof(BOX2DFLOAT4));
+ bbox->xmin = (float)(*dptr); dptr++;
+ bbox->xmax = (float)(*dptr); dptr++;
+ bbox->ymin = (float)(*dptr); dptr++;
+ bbox->ymax = (float)(*dptr); dptr++;
+ lwgeom->bbox = bbox;
+ }
+ else
+ {
+ lwgeom->bbox = NULL;
+ }
+
+ if( has_srid )
+ lwgeom->SRID = g_srid;
+ else
+ lwgeom->SRID = -1;
+
+ return lwgeom;
+}
+
+/***********************************************************************
+* Calculate geocentric bounding box from geodetic coordinates
+* of GSERIALIZED. To be used in index calculations to get the box
+* of smaller features on the fly, and in feature creation, to
+* calculate the box that will be added to the feature.
+*/
+
+static int gserialized_calculate_gbox_geocentric_from_any(uchar *data_ptr, size_t *g_size, GBOX *gbox);
+
+static int gserialized_calculate_gbox_geocentric_from_point(uchar *data_ptr, size_t *g_size, GBOX *gbox)
+{
+ uchar *start_ptr = data_ptr;
+ int npoints = 0;
+ POINTARRAY *pa;
+
+ assert(data_ptr);
+
+ data_ptr += 4; /* Move past type integer. */
+ npoints = lw_get_uint32(data_ptr);
+ data_ptr += 4; /* Move past npoints. */
+
+ if ( npoints == 0 ) /* Empty point */
+ {
+ if(g_size) *g_size = data_ptr - start_ptr;
+ return G_FAILURE;
+ }
+
+ pa = pointArray_construct(data_ptr, FLAGS_GET_Z(gbox->flags), FLAGS_GET_M(gbox->flags), npoints);
+
+ if( ptarray_calculate_gbox_geodetic(pa, gbox) == G_FAILURE )
+ return G_FAILURE;
+
+ /* Move past all the double ordinates. */
+ data_ptr += sizeof(double) * FLAGS_NDIMS(gbox->flags);
+
+ if(g_size)
+ *g_size = data_ptr - start_ptr;
+
+ lwfree(pa);
+
+ return G_SUCCESS;
+}
+
+static int gserialized_calculate_gbox_geocentric_from_line(uchar *data_ptr, size_t *g_size, GBOX *gbox)
+{
+ uchar *start_ptr = data_ptr;
+ int npoints = 0;
+ POINTARRAY *pa;
+
+ assert(data_ptr);
+
+ data_ptr += 4; /* Move past type integer. */
+ npoints = lw_get_uint32(data_ptr);
+ data_ptr += 4; /* Move past npoints. */
+
+ if ( npoints == 0 ) /* Empty linestring */
+ {
+ if(g_size) *g_size = data_ptr - start_ptr;
+ return G_FAILURE;
+ }
+
+ pa = pointArray_construct(data_ptr, FLAGS_GET_Z(gbox->flags), FLAGS_GET_M(gbox->flags), npoints);
+
+ if( ptarray_calculate_gbox_geodetic(pa, gbox) == G_FAILURE )
+ return G_FAILURE;
+
+ /* Move past all the double ordinates. */
+ data_ptr += sizeof(double) * FLAGS_NDIMS(gbox->flags) * npoints;
+
+ if(g_size)
+ *g_size = data_ptr - start_ptr;
+
+ lwfree(pa);
+
+ return G_SUCCESS;
+}
+
+static int gserialized_calculate_gbox_geocentric_from_polygon(uchar *data_ptr, size_t *g_size, GBOX *gbox)
+{
+ uchar *start_ptr = data_ptr;
+ int npoints0 = 0; /* Points in exterior ring. */
+ int npoints = 0; /* Points in all rings. */
+ int nrings = 0;
+ POINTARRAY *pa;
+ int i;
+
+ assert(data_ptr);
+
+ data_ptr += 4; /* Move past type integer. */
+
+ nrings = lw_get_uint32(data_ptr);
+ data_ptr += 4; /* Move past nrings. */
+
+ if( nrings <= 0 )
+ {
+ if(g_size) *g_size = data_ptr - start_ptr;
+ return G_FAILURE; /* Empty polygon */
+ }
+
+ npoints0 = lw_get_uint32(data_ptr); /* NPoints in first (exterior) ring. */
+
+ for( i = 0; i < nrings; i++ )
+ {
+ npoints += lw_get_uint32(data_ptr);
+ data_ptr += 4; /* Move past this npoints value. */
+ }
+
+ if ( nrings % 2 ) /* Move past optional padding. */
+ data_ptr += 4;
+
+ pa = pointArray_construct(data_ptr, FLAGS_GET_Z(gbox->flags), FLAGS_GET_M(gbox->flags), npoints);
+
+ /* Bounds of exterior ring is bounds of whole polygon. */
+ if( ptarray_calculate_gbox_geodetic(pa, gbox) == G_FAILURE )
+ return G_FAILURE;
+
+ /* Move past all the double ordinates. */
+ data_ptr += sizeof(double) * FLAGS_NDIMS(gbox->flags) * npoints;
+
+ if(g_size)
+ *g_size = data_ptr - start_ptr;
+
+ lwfree(pa);
+
+ return G_SUCCESS;
+}
+
+static int gserialized_calculate_gbox_geocentric_from_collection(uchar *data_ptr, size_t *g_size, GBOX *gbox)
+{
+ uchar *start_ptr = data_ptr;
+ int ngeoms = 0;
+ int i;
+ int first = G_TRUE;
+ int result = G_FAILURE;
+
+ assert(data_ptr);
+
+ data_ptr += 4; /* Move past type integer. */
+ ngeoms = lw_get_uint32(data_ptr);
+ data_ptr += 4; /* Move past ngeoms. */
+
+ if( ngeoms <= 0 ) return G_FAILURE; /* Empty collection */
+
+ for( i = 0; i < ngeoms; i++ )
+ {
+ size_t subgeom_size = 0;
+ GBOX subbox;
+ subbox.flags = gbox->flags;
+ if ( gserialized_calculate_gbox_geocentric_from_any(data_ptr, &subgeom_size, &subbox) != G_FAILURE )
+ {
+ if ( first )
+ {
+ gbox_duplicate(&subbox, gbox);
+ first = G_FALSE;
+ }
+ else
+ {
+ gbox_merge(gbox, &subbox);
+ }
+ result = G_SUCCESS;
+ }
+ data_ptr += subgeom_size;
+ }
+
+ if(g_size)
+ *g_size = data_ptr - start_ptr;
+
+ return result;
+}
+
+static int gserialized_calculate_gbox_geocentric_from_any(uchar *data_ptr, size_t *g_size, GBOX *gbox)
+{
+
+ uint32 type;
+
+ assert(data_ptr);
+
+ type = lw_get_uint32(data_ptr);
+
+ LWDEBUGF(2, "Got type %d (%s)", type, lwgeom_typename(type));
+ LWDEBUGF(3, "Got gbox pointer (%p)", gbox);
+
+ switch (type)
+ {
+ case POINTTYPE:
+ return gserialized_calculate_gbox_geocentric_from_point(data_ptr, g_size, gbox);
+ case LINETYPE:
+ return gserialized_calculate_gbox_geocentric_from_line(data_ptr, g_size, gbox);
+ case POLYGONTYPE:
+ return gserialized_calculate_gbox_geocentric_from_polygon(data_ptr, g_size, gbox);
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTIPOLYGONTYPE:
+ case COMPOUNDTYPE:
+ case CURVEPOLYTYPE:
+ case MULTICURVETYPE:
+ case MULTISURFACETYPE:
+ case COLLECTIONTYPE:
+ return gserialized_calculate_gbox_geocentric_from_collection(data_ptr, g_size, gbox);
+ default:
+ lwerror("Unsupported geometry type: %d", type);
+ return G_FAILURE;
+ }
+}
+
+int gserialized_calculate_gbox_geocentric_p(GSERIALIZED *g, GBOX *g_box)
+{
+ uchar *data_ptr = NULL;
+ size_t g_size = 0;
+ int result = G_SUCCESS;
+
+ assert(g);
+
+ /* This function only works for geodetics. */
+ if( ! FLAGS_GET_GEODETIC(g->flags) )
+ {
+ lwerror("Function only accepts geodetic inputs.");
+ return G_FAILURE;
+ }
+
+ LWDEBUGF(4, "Got input %s", gserialized_to_string(g));
+
+ data_ptr = (uchar*)g->data;
+ g_box->flags = g->flags;
+
+ /* If the serialized form already has a box, skip past it. */
+ if( FLAGS_GET_BBOX(g->flags) )
+ {
+ int ndims = FLAGS_GET_GEODETIC(g->flags) ? 3 : FLAGS_NDIMS(g->flags);
+ data_ptr += 2 * ndims * sizeof(float); /* Copy the bounding box and return. */
+ LWDEBUG(3,"Serialized form has box already, skipping past...");
+ }
+
+ LWDEBUG(3,"Calculating box...");
+ /* Calculate the bounding box from the geometry. */
+ result = gserialized_calculate_gbox_geocentric_from_any(data_ptr, &g_size, g_box);
+
+ if( result == G_FAILURE )
+ {
+ LWDEBUG(3,"Unable to calculate geocentric bounding box.");
+ return G_FAILURE; /* Ooops! */
+ }
+
+ LWDEBUGF(3,"Returning box: %s", gbox_to_string(g_box));
+ return result;
+}
+
+GBOX* gserialized_calculate_gbox_geocentric(GSERIALIZED *g)
+{
+ GBOX g_box;
+ int result = G_SUCCESS;
+
+ result = gserialized_calculate_gbox_geocentric_p(g, &g_box);
+
+ if( result == G_FAILURE )
+ return NULL; /* Ooops! */
+
+ LWDEBUGF(3,"Returning box: %s", gbox_to_string(&g_box));
+ return gbox_copy(&g_box);
+}
+
+
+GSERIALIZED* gserialized_copy(GSERIALIZED *g)
+{
+ GSERIALIZED *g_out = NULL;
+ assert(g);
+ g_out = (GSERIALIZED*)lwalloc(SIZE_GET(g->size));
+ memcpy((uchar*)g_out,(uchar*)g,SIZE_GET(g->size));
+ return g_out;
+}
+
+char* gserialized_to_string(GSERIALIZED *g)
+{
+ LWGEOM_UNPARSER_RESULT lwg_unparser_result;
+ LWGEOM *lwgeom = lwgeom_from_gserialized(g);
+ uchar *serialized_lwgeom;
+ int result;
+
+ assert(g);
+
+ if( ! lwgeom )
+ {
+ lwerror("Unable to create lwgeom from gserialized");
+ return NULL;
+ }
+
+ serialized_lwgeom = lwgeom_serialize(lwgeom);
+ lwgeom_release(lwgeom);
+
+ result = serialized_lwgeom_to_ewkt(&lwg_unparser_result, serialized_lwgeom, PARSER_CHECK_NONE);
+ lwfree(serialized_lwgeom);
+
+ return lwg_unparser_result.wkoutput;
+
+}
+
--- /dev/null
+
+GSERIALIZED FORM
+=================
+
+The new serialized form, used by GEOGRAPHY, attempts to learn from the lessons of the SERIALIZED_LWGEOM, making slightly different trade-offs between alignment and overall compactness.
+
+* It is understood that GSERIALIZED is be used primarily (only) by PostGIS. Other users of the geometry library (for example, the loaders and dumpers) will serialize to WKB or EWKB or various text representations. Therefore, GSERIALIZED includes the uint32 "size" field at the top used by PostgreSQL for varlena types.
+* Like SERIALIZED_LWGEOM, GSERIALIZED is built to be read and written recursively, so that nested collections of collections (of ...) are possible without restrictions on depth.
+* GSERIALIZED preserves double alignment of ordinate arrays. This will allow coordinate access without memcpy.
+* GSERIALIZED includes a mandatory SRID, in recognition of the fact that most production use of PostGIS does actually use an SRID. In SERIALIZED_LWGEOM the SRID is optional.
+* GSERIALIZED places the dimensionality information, the SRID information and the bounding boxes at the front of the structure, and all sub-components inherit from those parent values.
+* GSERIALIZED retains the idea of optional bounding boxes, so that small features do not carry the extra storage overhead of a largely redundant bounding box.
+
+STRUCTURE
+---------
+
+Most of the internals of GSERIALIZED are anonymous. One thing that differs from SERIALIZED_LWGEOM is that the geometry type is no longer directly accessible from the structure (it was inside the one byte "type" at the top of the SERIALIZED_LWGEOM). To access the type in GSERIALIZED, you first have to figure out whether there is a an bounding box, how many dimensions that bounding box has, and move the data pointer appropriately before dereferencing. The gserialized_get_type(GSERIALIZED *s) function carries out this operation.
+
+typedef struct
+{
+ uint32 size; /* For PgSQL use, use VAR* macros to manipulate. */
+ uchar srid[3]; /* 24 bits of SRID */
+ uchar flags; /* HasZ, HasM, HasBBox, IsGeodetic */
+ uchar data[1]; /* See gserialized.txt */
+} GSERIALIZED;
+
+The standard header is as follows (using <> to denote 4-byte fields and [] to denote 8-byte fields):
+
+<size> size /* Used by PgSQL */
+<srid /* 3 bytes */
+ flags> /* 1 byte */
+<bbox-xmin> /* bounding box is optional */
+<bbox-xmax> /* number of dimensions is variable and indicated in the flags */
+<bbox-ymin>
+<bbox-ymax>
+
+After the header comes the recursively searchable geometry representations. Each type is double aligned, so any combination is double aligned, and we start after a double aligned standard header, so we are golden:
+
+<pointype>
+<npoints> /* 0 if empty, 1 otherwise */
+[double]
+[double]
+[double]
+
+<linestringtype>
+<npoints> /* 0 if empty, N otherwise */
+[double]
+...
+[double]
+
+<polygontype>
+<nrings> /* 0 if empty, N otherwise */
+<npointsring1>
+<npointsring2>
+<?pad?> /* pad if nrings % 2 > 0 */
+[double]
+...
+[double]
+
+<collectiontype>
+<ngeoms> /* 0 if empty, N otherwise */
+[geom]
+...
+[geom]
+
+<circularstringtype>
+<npoints> /* 0 if empty, N otherwise */
+[double]
+...
+[double]
+
+<compoundstringtype>
+<ngeoms> /* 0 if empty, N otherwise */
+[geom]
+...
+[geom]
+
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "libgeom.h"
+
+uchar gflags(int hasz, int hasm, int geodetic)
+{
+ unsigned char flags = 0;
+ if ( hasz )
+ flags = FLAGS_SET_Z(flags, 1);
+ if ( hasm )
+ flags = FLAGS_SET_M(flags, 1);
+ if ( geodetic )
+ flags = FLAGS_SET_GEODETIC(flags, 1);
+ return flags;
+}
+
+/**
+* Calculate type integer and dimensional flags from string input.
+* Case insensitive, and insensitive to spaces at front and back.
+* Type == 0 in the case of the string "GEOMETRY" or "GEOGRAPHY".
+* Return G_SUCCESS for success.
+*/
+int geometry_type_from_string(char *str, int *type, int *z, int *m)
+{
+ regex_t rx_point;
+ regex_t rx_linestring;
+ regex_t rx_polygon;
+ regex_t rx_geometrycollection;
+ regex_t rx_geometry;
+ size_t rx_nmatch = 5;
+ regmatch_t rx_matchptr[5];
+
+ assert(str);
+ assert(type);
+ assert(z);
+ assert(m);
+
+ /* Initialize. */
+ *type = 0;
+ *z = 0;
+ *m = 0;
+
+ regcomp(&rx_point, "^ *(st_)?(multi)?point(z)?(m)? *$", REG_ICASE | REG_EXTENDED);
+ regcomp(&rx_linestring, "^ *(st_)?(multi)?linestring(z)?(m)? *$", REG_ICASE | REG_EXTENDED);
+ regcomp(&rx_polygon, "^ *(st_)?(multi)?polygon(z)?(m)? *$", REG_ICASE | REG_EXTENDED);
+ regcomp(&rx_geometrycollection, "^ *(st_)?geometrycollection(z)?(m)? *$", REG_ICASE | REG_EXTENDED);
+ regcomp(&rx_geometry, "^ *(st_)?geometry(z)?(m)? *$", REG_ICASE | REG_EXTENDED);
+
+ if( ! regexec(&rx_point, str, rx_nmatch, rx_matchptr, 0) )
+ {
+ *type = POINTTYPE;
+ if(rx_matchptr[2].rm_so != -1) /* MULTI */
+ *type = MULTIPOINTTYPE;
+
+ if(rx_matchptr[3].rm_so != -1) /* Z */
+ *z = 1;
+
+ if(rx_matchptr[4].rm_so != -1) /* M */
+ *m = 1;
+
+ return G_SUCCESS;
+ }
+ if( ! regexec(&rx_linestring, str, rx_nmatch, rx_matchptr, 0) )
+ {
+ *type = LINETYPE;
+ if(rx_matchptr[2].rm_so != -1) /* MULTI */
+ *type = MULTILINETYPE;
+
+ if(rx_matchptr[3].rm_so != -1) /* Z */
+ *z = 1;
+
+ if(rx_matchptr[4].rm_so != -1) /* M */
+ *m = 1;
+
+ return G_SUCCESS;
+ }
+ if( ! regexec(&rx_polygon, str, rx_nmatch, rx_matchptr, 0) )
+ {
+ *type = POLYGONTYPE;
+ if(rx_matchptr[2].rm_so != -1) /* MULTI */
+ *type = MULTIPOLYGONTYPE;
+
+ if(rx_matchptr[3].rm_so != -1) /* Z */
+ *z = 1;
+
+ if(rx_matchptr[4].rm_so != -1) /* M */
+ *m = 1;
+
+ return G_SUCCESS;
+ }
+ if( ! regexec(&rx_geometrycollection, str, rx_nmatch, rx_matchptr, 0) )
+ {
+ *type = COLLECTIONTYPE;
+
+ if(rx_matchptr[2].rm_so != -1) /* Z */
+ *z = 1;
+
+ if(rx_matchptr[3].rm_so != -1) /* M */
+ *m = 1;
+
+ return G_SUCCESS;
+ }
+ if( ! regexec(&rx_geometry, str, rx_nmatch, rx_matchptr, 0) )
+ {
+ *type = 0; /* Generic geometry type. */
+
+ if(rx_matchptr[2].rm_so != -1) /* Z */
+ *z = 1;
+
+ if(rx_matchptr[3].rm_so != -1) /* M */
+ *m = 1;
+
+ return G_SUCCESS;
+ }
+
+ return G_FAILURE;
+
+}
+
+
+
+
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ * Copyright 2001-2006 Refractions Research Inc.
+ * Copyright 2007-2008 Mark Cave-Ayland
+ * Copyright 2008 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "liblwgeom.h"
+#include <string.h>
+#include <assert.h>
+#include <regex.h>
+
+/**
+* PI
+*/
+#define PI 3.1415926535897932384626433832795
+
+/**
+* Return types for functions with status returns.
+*/
+
+#define G_FAILURE 0
+#define G_SUCCESS 1
+#define G_TRUE 1
+#define G_FALSE 0
+
+/**
+* Maximum allowed SRID value.
+* Currently we are using 20 bits (1048575) of storage for SRID.
+*/
+#define SRID_MAXIMUM 999999
+
+/**
+* Macros for manipulating the 'flags' byte. A uchar used as follows:
+* ---RGBMZ
+* Three unused bits, followed by ReadOnly, Geodetic, HasBBox, HasM and HasZ flags.
+*/
+#define FLAGS_GET_Z(flags) ((flags) & 0x01)
+#define FLAGS_GET_M(flags) (((flags) & 0x02)>>1)
+#define FLAGS_GET_BBOX(flags) (((flags) & 0x4)>>2)
+#define FLAGS_GET_GEODETIC(flags) (((flags) & 0x08)>>3)
+#define FLAGS_GET_READONLY(flags) (((flags) & 0x10)>>4)
+#define FLAGS_SET_Z(flags, value) ((value) ? ((flags) | 0x01) : ((flags) & 0xFE))
+#define FLAGS_SET_M(flags, value) ((value) ? ((flags) | 0x02) : ((flags) & 0xFD))
+#define FLAGS_SET_BBOX(flags, value) ((value) ? ((flags) | 0x04) : ((flags) & 0xFB))
+#define FLAGS_SET_GEODETIC(flags, value) ((value) ? ((flags) | 0x08) : ((flags) & 0xF7))
+#define FLAGS_SET_READONLY(flags, value) ((value) ? ((flags) | 0x10) : ((flags) & 0xEF))
+#define FLAGS_NDIMS(flags) (2 + FLAGS_GET_Z(flags) + FLAGS_GET_M(flags))
+
+/**
+* Macro for reading the size from the GSERIALIZED size attribute.
+* Cribbed from PgSQL, top 30 bits are size. Use VARSIZE() when working
+* internally with PgSQL.
+*/
+#define SIZE_GET(varsize) (((varsize) >> 2) & 0x3FFFFFFF)
+#define SIZE_SET(varsize, size) (((varsize) & 0x00000003)|(((size) & 0x3FFFFFFF) << 2 ))
+
+/**
+* Macros for manipulating the 'typemod' int. An int32 used as follows:
+* Plus/minus = Top bit.
+* Spare bits = Next 3 bits.
+* SRID = Next 20 bits.
+* TYPE = Next 6 bits.
+* ZM Flags = Bottom 2 bits.
+*/
+#define TYPMOD_GET_SRID(typmod) ((typmod & 0x0FFFFF00)>>8)
+#define TYPMOD_SET_SRID(typmod, srid) ((typmod & 0x000000FF) | ((srid & 0x000FFFFF)<<8))
+#define TYPMOD_GET_TYPE(typmod) ((typmod & 0x000000FC)>>2)
+#define TYPMOD_SET_TYPE(typmod, type) ((typmod & 0xFFFFFF03) | ((type & 0x0000003F)<<2))
+#define TYPMOD_GET_Z(typmod) ((typmod & 0x00000002)>>1)
+#define TYPMOD_SET_Z(typmod) (typmod | 0x00000002)
+#define TYPMOD_GET_M(typmod) (typmod & 0x00000001)
+#define TYPMOD_SET_M(typmod) (typmod | 0x00000001)
+#define TYPMOD_GET_NDIMS(typmod) (2+TYPMOD_GET_Z(typmod)+TYPMOD_GET_M(typmod))
+
+
+
+/**
+* GBOX structure.
+* Include the flags, so we don't have to constantly pass them
+* into the functions.
+*/
+typedef struct
+{
+ uchar flags;
+ double xmin;
+ double xmax;
+ double ymin;
+ double ymax;
+ double zmin;
+ double zmax;
+ double mmin;
+ double mmax;
+} GBOX;
+
+
+
+/**
+* Coordinate structure. This will be created on demand from GPTARRAY when
+* needed by algorithms and whatnot.
+*/
+typedef struct
+{
+ uchar flags;
+ double *ordinates;
+} GCOORDINATE;
+
+
+
+/* Start with space for 16 points */
+#define G_PT_ARRAY_DEFAULT_POINTS 16
+
+/**
+* GPTARRAY structure. Holder of our ordinates. Dynamically sized
+* container for coordinate arrays.
+*/
+typedef struct
+{
+ uchar flags;
+ size_t capacity;
+ uint32 npoints;
+ double *ordinates;
+} GPTARRAY;
+
+
+/**
+* GSERIALIZED
+*/
+typedef struct
+{
+ uint32 size; /* For PgSQL use only, use VAR* macros to manipulate. */
+ uchar srid[3]; /* 24 bits of SRID */
+ uchar flags; /* HasZ, HasM, HasBBox, IsGeodetic, IsReadOnly */
+ uchar data[1]; /* See gserialized.txt */
+} GSERIALIZED;
+
+/**
+* G_GEOMETRY
+*
+* flags = bit flags, see FLAGS_HAS_* defines
+* type = enumerated OGC geometry type number
+*/
+typedef struct
+{
+ uchar flags;
+ uint32 type;
+ GBOX *bbox; /* NULL == unneeded */
+ uint32 srid; /* 0 == unknown */
+ void *data;
+} G_GEOMETRY;
+
+/* POINTTYPE */
+typedef struct
+{
+ uchar flags;
+ uint32 type;
+ GBOX *bbox; /* NULL == unneeded */
+ uint32 srid; /* 0 == unknown */
+ GPTARRAY *point;
+} G_POINT;
+
+/* LINETYPE */
+/* CIRCSTRINGTYPE */
+typedef struct
+{
+ uchar flags;
+ uint32 type;
+ GBOX *bbox; /* NULL == unneeded */
+ uint32 srid; /* 0 == unknown */
+ GPTARRAY *points;
+} G_LINESTRING;
+
+typedef G_LINESTRING G_CIRCULARSTRING;
+
+/* POLYGONTYPE */
+typedef struct
+{
+ uchar flags;
+ uint32 type;
+ GBOX *bbox; /* NULL == unneeded */
+ uint32 srid; /* 0 == unknown */
+ size_t capacity; /* How much space is allocated for *rings? */
+ uint32 nrings;
+ GPTARRAY **rings; /* rings[0] is the exterior ring */
+} G_POLYGON;
+
+/* MULTIPOINTTYPE */
+/* MULTILINETYPE */
+/* MULTIPOINTTYPE */
+/* COMPOUNDTYPE */
+/* CURVEPOLYTYPE */
+/* MULTICURVETYPE */
+/* MULTISURFACETYPE */
+/* COLLECTIONTYPE */
+typedef struct
+{
+ uchar flags;
+ uint32 type;
+ GBOX *bbox; /* NULL == unneeded */
+ uint32 srid; /* 0 == unknown */
+ size_t capacity; /* How much space is allocated for *geoms? */
+ uint32 ngeoms;
+ G_GEOMETRY **geoms;
+} G_COLLECTION;
+
+/*
+** All the collection types share the same physical structure as the
+** generic geometry collection. We add type aliases so we can be more
+** explicit in our functions later.
+*/
+typedef G_COLLECTION G_MPOINT;
+typedef G_COLLECTION G_MLINESTRING;
+typedef G_COLLECTION G_MPOLYGON;
+typedef G_COLLECTION G_MSURFACE;
+typedef G_COLLECTION G_MCURVE;
+typedef G_COLLECTION G_COMPOUNDCURVE;
+typedef G_COLLECTION G_CURVEPOLYGON;
+
+/*
+* Utility casts from GEOMETRY to concrete type.
+* Return NULL if cast is illegal.
+extern G_POINT* g_geometry_as_point(G_GEOMETRY *g);
+extern G_LINESTRING* g_geometry_as_linestring(G_GEOMETRY *g);
+extern G_POLYGON* g_geometry_as_polygon(G_GEOMETRY *g);
+extern G_MPOINT* g_geometry_as_mpoint(G_GEOMETRY *g);
+extern G_MLINESTRING* g_geometry_as_mlinestring(G_GEOMETRY *g);
+extern G_MPOLYGON* g_geometry_as_mpolygon(G_GEOMETRY *g);
+*/
+
+/*
+* Utility casts from concrete type to GEOMETRY.
+* Return NULL if cast is illegal.
+extern G_GEOMETRY* g_point_as_geometry(G_POINT *g);
+extern G_GEOMETRY* g_linestring_as_geometry(G_LINESTRING *g);
+extern G_GEOMETRY* g_polygon_as_geometry(G_POLYGON *g);
+extern G_GEOMETRY* g_mpoint_as_geometry(G_MPOINT *g);
+extern G_GEOMETRY* g_mlinestring_as_geometry(G_MLINESTRING *g);
+extern G_GEOMETRY* g_mpolygon_as_geometry(G_POLYGON *g);
+extern G_GEOMETRY* g_collection_as_geometry(G_COLLECTION *g);
+*/
+
+/***********************************************************************
+* Coordinate creator, set, get and destroy.
+*/
+extern GCOORDINATE* gcoord_new(int ndims);
+extern GCOORDINATE* gcoord_new_with_flags(uchar flags);
+extern GCOORDINATE* gcoord_new_with_flags_and_ordinates(uchar flags, double *ordinates);
+extern GCOORDINATE* gcoord_copy(GCOORDINATE *coord);
+extern void gcoord_free(GCOORDINATE *coord);
+extern void gcoord_set_x(GCOORDINATE *coord, double x);
+extern void gcoord_set_y(GCOORDINATE *coord, double y);
+extern void gcoord_set_z(GCOORDINATE *coord, double z);
+extern void gcoord_set_m(GCOORDINATE *coord, double m);
+extern void gcoord_set_ordinates(GCOORDINATE *coord, double *ordinates);
+extern void gcoord_set_ordinate(GCOORDINATE *coord, double ordinate, int index);
+extern double gcoord_get_x(GCOORDINATE *coord);
+extern double gcoord_get_y(GCOORDINATE *coord);
+extern double gcoord_get_z(GCOORDINATE *coord);
+extern double gcoord_get_m(GCOORDINATE *coord);
+
+/***********************************************************************
+* Point arrays, set, get and destroy.
+*/
+extern GPTARRAY* gptarray_new(uchar flags);
+extern GPTARRAY* gptarray_new_with_size(uchar flags, int npoints);
+extern GPTARRAY* gptarray_new_with_ordinates(uchar flags, int npoints, double *ordinates);
+extern GPTARRAY* gptarray_copy(GPTARRAY *ptarray);
+extern void gptarray_free(GPTARRAY *ptarray);
+extern void gptarray_add_coord(GPTARRAY *ptarray, GCOORDINATE *coord);
+extern GCOORDINATE* gptarray_get_coord_ro(GPTARRAY *ptarray, int i);
+extern GCOORDINATE* gptarray_get_coord_new(GPTARRAY *ptarray, int i);
+extern void gptarray_set_coord(GPTARRAY *ptarray, int i, GCOORDINATE *coord);
+extern void gptarray_set_x(GPTARRAY *ptarray, int i, double x);
+extern void gptarray_set_y(GPTARRAY *ptarray, int i, double y);
+extern void gptarray_set_z(GPTARRAY *ptarray, int i, double z);
+extern void gptarray_set_m(GPTARRAY *ptarray, int i, double m);
+extern double gptarray_get_x(GPTARRAY *ptarray, int i);
+extern double gptarray_get_y(GPTARRAY *ptarray, int i);
+extern double gptarray_get_z(GPTARRAY *ptarray, int i);
+extern double gptarray_get_m(GPTARRAY *ptarray, int i);
+
+
+/***********************************************************************
+** Linestrings
+*/
+extern G_LINESTRING* glinestring_new_from_gptarray(GPTARRAY *ptarray);
+extern G_LINESTRING* glinestring_new(uchar flags);
+
+
+/***********************************************************************
+** Utility functions for flag byte and srid_flag integer.
+*/
+
+/**
+* Construct a new flags char.
+*/
+extern uchar gflags(int hasz, int hasm, int geodetic);
+
+/**
+* Extract the geometry type from the serialized form (it hides in
+* the anonymous data area, so this is a handy function).
+*/
+extern uint32 gserialized_get_type(GSERIALIZED *g);
+
+/**
+* Extract the SRID from the serialized form (it is packed into
+* three bytes so this is a handy function).
+*/
+extern uint32 gserialized_get_srid(GSERIALIZED *g);
+
+/**
+* Write the SRID into the serialized form (it is packed into
+* three bytes so this is a handy function).
+*/
+extern void gserialized_set_srid(GSERIALIZED *g, uint32 srid);
+
+
+/***********************************************************************
+** Functions for managing serialized forms and bounding boxes.
+*/
+
+/**
+* Calculate the geocentric bounding box directly from the serialized
+* form of the geodetic coordinates. Only accepts serialized geographies
+* flagged as geodetic. Caller is responsible for disposing of the GBOX.
+*/
+extern GBOX* gserialized_calculate_gbox_geocentric(GSERIALIZED *g);
+
+/**
+* Calculate the geocentric bounding box directly from the serialized
+* form of the geodetic coordinates. Only accepts serialized geographies
+* flagged as geodetic.
+*/
+int gserialized_calculate_gbox_geocentric_p(GSERIALIZED *g, GBOX *g_box);
+
+/**
+* Return a WKT representation of the gserialized geometry.
+* Caller is responsible for disposing of the char*.
+*/
+extern char* gserialized_to_string(GSERIALIZED *g);
+
+/**
+* Return a copy of the input serialized geometry.
+*/
+extern GSERIALIZED* gserialized_copy(GSERIALIZED *g);
+
+/**
+* Check that coordinates of LWGEOM are all within the geodetic range.
+*/
+extern int lwgeom_check_geodetic(const LWGEOM *geom);
+
+/**
+* Calculate the geodetic bounding box for an LWGEOM. Z/M coordinates are
+* ignored for this calculation. Pass in non-null, geodetic bounding box for function
+* to fill out. LWGEOM must have been built from a GSERIALIZED to provide
+* double aligned point arrays.
+*/
+extern int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox);
+
+/**
+* Calculate the 2-4D bounding box of a geometry. Z/M coordinates are honored
+* for this calculation, though for curves they are not included in calculations
+* of curvature.
+*/
+extern int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox);
+
+/**
+* New function to read doubles directly from the double* coordinate array
+* of an aligned lwgeom #POINTARRAY (built by de-serializing a #GSERIALIZED).
+*/
+extern int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point);
+
+/**
+* @brief Check whether or not a lwgeom is big enough to warrant a bounding box.
+*
+* Check whether or not a lwgeom is big enough to warrant a bounding box
+* when stored in the serialized form on disk. Currently only points are
+* considered small enough to not require a bounding box, because the
+* index operations can generate a large number of box-retrieval operations
+* when scanning keys.
+*/
+extern int lwgeom_needs_bbox(LWGEOM *geom);
+
+/**
+* Count the total number of vertices in any #LWGEOM.
+*/
+extern int lwgeom_count_vertices(LWGEOM *geom);
+
+/**
+* Calculate box and add values to gbox. Return #G_SUCCESS on success.
+*/
+extern int ptarray_calculate_gbox_geodetic(POINTARRAY *pa, GBOX *gbox);
+
+/**
+* Create a new gbox with the dimensionality indicated by the flags. Caller
+* is responsible for freeing.
+*/
+extern GBOX* gbox_new(uchar flags);
+
+/**
+* Update the merged #GBOX to be large enough to include itself and the new box.
+*/
+extern int gbox_merge(GBOX *merged_box, GBOX *new_box);
+
+/**
+* Allocate a string representation of the #GBOX, based on dimensionality of flags.
+*/
+extern char* gbox_to_string(GBOX *gbox);
+
+/**
+* Return a copy of the #GBOX, based on dimensionality of flags.
+*/
+extern GBOX* gbox_copy(GBOX *gbox);
+
+/**
+* Return #G_TRUE if the #GBOX overlaps, #G_FALSE otherwise.
+*/
+extern int gbox_overlaps(GBOX *g1, GBOX *g2);
+
+/**
+* Copy the values of original #GBOX into duplicate.
+*/
+extern void gbox_duplicate(GBOX *original, GBOX *duplicate);
+
+/**
+* Return the number of bytes necessary to hold a #GBOX of this dimension in
+* serialized form.
+*/
+extern size_t gbox_serialized_size(uchar flags);
+
+/**
+* Utility function to get type number from string. For example, a string 'POINTZ'
+* would return type of 1 and z of 1 and m of 0. Valid
+*/
+extern int geometry_type_from_string(char *str, int *type, int *z, int *m);
+
+/**
+* Calculate required memory segment to contain a serialized form of the LWGEOM.
+* Primarily used internally by the serialization code. Exposed to allow the cunit
+* tests to exercise it.
+*/
+extern size_t gserialized_from_lwgeom_size(LWGEOM *geom, GBOX *gbox);
+
+/**
+* Allocate a new #GSERIALIZED from an #LWGEOM. For all non-point types, a bounding
+* box will be calculated and embedded in the serialization. The geodetic flag is used
+* to control the box calculation (cartesian or geocentric). If set, the size pointer
+* will contain the size of the final output, which is useful for setting the PgSQL
+* VARSIZE information.
+*/
+extern GSERIALIZED* gserialized_from_lwgeom(LWGEOM *geom, int is_geodetic, size_t *size);
+
+/**
+* Allocate a new #LWGEOM from a #GSERIALIZED. The resulting #LWGEOM will have coordinates
+* that are double aligned and suitable for direct reading using getPoint2d_p_ro
+*/
+extern LWGEOM* lwgeom_from_gserialized(GSERIALIZED *g);
+
+/**
+* Serialize/deserialize from/to #G_GEOMETRY into #GSERIALIZED
+extern size_t gserialized_from_ggeometry_size(G_GEOMETRY *geom);
+extern GSERIALIZED* gserialized_from_ggeometry(G_GEOMETRY *geom);
+extern G_GEOMETRY* ggeometry_from_gserialized(GSERIALIZED *g);
+*/
+
#define NO_Z_VALUE NO_VALUE
#define NO_M_VALUE NO_VALUE
+#ifndef MAXFLOAT
+#define MAXFLOAT 3.402823466e+38F
+#endif
+
#ifndef C_H
typedef unsigned int uint32;
extern BOX3D *ptarray_compute_box3d(const POINTARRAY *pa);
extern int ptarray_compute_box3d_p(const POINTARRAY *pa, BOX3D *out);
+
+
/*
* size of point represeneted in the POINTARRAY
* 16 for 2d, 24 for 3d, 32 for 4d
-#ifndef MAXFLOAT
-#define MAXFLOAT 3.402823466e+38F
-#endif
-
/*
* Construct a new LWCIRCSTRING. points will *NOT* be copied
* use SRID=-1 for unknown SRID (will have 8bit type's S = 0)
return size;
}
+
BOX3D *
lwcircle_compute_box3d(POINT4D *p1, POINT4D *p2, POINT4D *p3)
{
LWDEBUGF(4, "circularstring %d x=(%.16f,%.16f) y=(%.16f,%.16f) z=(%.16f,%.16f)", i/2, box->xmin, box->xmax, box->ymin, box->ymax, box->zmin, box->zmax);
}
-
return box;
}
+
int
lwcircstring_compute_box2d_p(LWCIRCSTRING *curve, BOX2DFLOAT4 *result)
{
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include <math.h>
+#include "libgeom.h"
+
+/*
+* This function can only be used on LWGEOM that is built on top of
+* GSERIALIZED, otherwise alignment errors will ensue.
+*/
+int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point)
+{
+ uchar *pa_ptr = NULL;
+ assert(pa);
+ assert(n >= 0);
+ assert(n < pa->npoints);
+
+ pa_ptr = getPoint_internal(pa, n);
+ //printf( "pa_ptr[0]: %g\n", *((double*)pa_ptr));
+ *point = (POINT2D*)pa_ptr;
+
+ return G_SUCCESS;
+}
+
+int ptarray_calculate_gbox_geodetic(POINTARRAY *pa, GBOX *gbox)
+{
+ //static double wgs84_radius = 6371008.7714;
+ int i;
+/* int pt_spacing; */
+ double *sphericalcoords;
+ POINT2D *pt;
+
+ assert(gbox);
+ assert(pa);
+
+/* pt_spacing = 2;
+ if( TYPE_HASZ(pa->dims) ) pt_spacing++;
+ if( TYPE_HASM(pa->dims) ) pt_spacing++;
+*/
+
+ if ( pa->npoints == 0 ) return G_FAILURE;
+
+ sphericalcoords = (double*)(pa->serialized_pointlist);
+
+ for( i = 0; i < pa->npoints; i++ )
+ {
+ getPoint2d_p_ro(pa, i, &pt);
+ /* Convert to radians */
+ /* register double theta = sphericalcoords[pt_spacing*i] * PI / 180.0; */
+ /* register double phi = sphericalcoords[pt_spacing*i + 1] * PI / 180; */
+ register double theta = pt->x * PI / 180.0;
+ register double phi = pt->y * PI / 180;
+ /* Convert to geocentric */
+ register double x = cos(theta);
+ register double y = sin(theta);
+ register double z = sin(phi);
+ /* Initialize the box */
+ if( i == 0 )
+ {
+ gbox->xmin = gbox->xmax = x;
+ gbox->ymin = gbox->ymax = y;
+ gbox->zmin = gbox->zmax = z;
+ gbox->mmin = gbox->mmax = 0.0;
+ }
+ /* Expand the box where necessary */
+ if( x < gbox->xmin ) gbox->xmin = x;
+ if( y < gbox->ymin ) gbox->ymin = y;
+ if( z < gbox->zmin ) gbox->zmin = z;
+ if( x > gbox->xmax ) gbox->xmax = x;
+ if( y > gbox->ymax ) gbox->ymax = y;
+ if( z > gbox->zmax ) gbox->zmax = z;
+ }
+
+ return G_SUCCESS;
+
+}
+
+static int lwpoint_calculate_gbox_geodetic(LWPOINT *point, GBOX *gbox)
+{
+ assert(point);
+ if( ptarray_calculate_gbox_geodetic(point->point, gbox) == G_FAILURE )
+ return G_FAILURE;
+ return G_SUCCESS;
+}
+
+static int lwline_calculate_gbox_geodetic(LWLINE *line, GBOX *gbox)
+{
+ assert(line);
+ if( ptarray_calculate_gbox_geodetic(line->points, gbox) == G_FAILURE )
+ return G_FAILURE;
+ return G_SUCCESS;
+}
+
+static int lwpolygon_calculate_gbox_geodetic(LWPOLY *poly, GBOX *gbox)
+{
+ GBOX ringbox;
+ int i;
+ int first = G_TRUE;
+ assert(poly);
+ if( poly->nrings == 0 )
+ return G_FAILURE;
+ ringbox.flags = gbox->flags;
+ for( i = 0; i < poly->nrings; i++ )
+ {
+ if( ptarray_calculate_gbox_geodetic(poly->rings[i], &ringbox) == G_FAILURE )
+ return G_FAILURE;
+ if( first )
+ {
+ gbox_duplicate(gbox, &ringbox);
+ first = G_FALSE;
+ }
+ else
+ {
+ gbox_merge(gbox, &ringbox);
+ }
+ }
+ return G_SUCCESS;
+}
+
+static int lwcollection_calculate_gbox_geodetic(LWCOLLECTION *coll, GBOX *gbox)
+{
+ GBOX subbox;
+ int i;
+ int result = G_FAILURE;
+ int first = G_TRUE;
+ assert(coll);
+ if( coll->ngeoms == 0 )
+ return G_FAILURE;
+
+ subbox.flags = gbox->flags;
+
+ for( i = 0; i < coll->ngeoms; i++ )
+ {
+ if( lwgeom_calculate_gbox_geodetic((LWGEOM*)(coll->geoms[i]), &subbox) == G_FAILURE )
+ {
+ continue;
+ }
+ else
+ {
+ if( first )
+ {
+ gbox_duplicate(gbox, &subbox);
+ first = G_FALSE;
+ }
+ else
+ {
+ gbox_merge(gbox, &subbox);
+ }
+ result = G_SUCCESS;
+ }
+ }
+ return result;
+}
+
+int lwgeom_calculate_gbox_geodetic(const LWGEOM *geom, GBOX *gbox)
+{
+ int result = G_FAILURE;
+ LWDEBUGF(4, "got type %d", TYPE_GETTYPE(geom->type));
+ if ( ! FLAGS_GET_GEODETIC(gbox->flags) )
+ {
+ lwerror("lwgeom_get_gbox_geodetic: non-geodetic gbox provided");
+ }
+ switch (TYPE_GETTYPE(geom->type))
+ {
+ case POINTTYPE:
+ result = lwpoint_calculate_gbox_geodetic((LWPOINT*)geom, gbox);
+ break;
+ case LINETYPE:
+ result = lwline_calculate_gbox_geodetic((LWLINE *)geom, gbox);
+ break;
+ case POLYGONTYPE:
+ result = lwpolygon_calculate_gbox_geodetic((LWPOLY *)geom, gbox);
+ break;
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTIPOLYGONTYPE:
+ case COLLECTIONTYPE:
+ result = lwcollection_calculate_gbox_geodetic((LWCOLLECTION *)geom, gbox);
+ break;
+ default:
+ lwerror("unsupported input geometry type: %d", TYPE_GETTYPE(geom->type));
+ break;
+ }
+ return result;
+}
+
+
+
+static int ptarray_check_geodetic(POINTARRAY *pa)
+{
+ int t;
+ POINT2D *pt;
+
+ assert(pa);
+
+ for (t=0; t<pa->npoints; t++)
+ {
+ getPoint2d_p_ro(pa, t, &pt);
+ //printf( "%d (%g, %g)\n", t, pt->x, pt->y);
+ if ( pt->x < -180.0 || pt->y < -90.0 || pt->x > 180.0 || pt->y > 90.0 )
+ return G_FALSE;
+ }
+
+ return G_TRUE;
+}
+
+static int lwpoint_check_geodetic(LWPOINT *point)
+{
+ assert(point);
+ return ptarray_check_geodetic(point->point);
+}
+
+static int lwline_check_geodetic(LWLINE *line)
+{
+ assert(line);
+ return ptarray_check_geodetic(line->points);
+}
+
+static int lwpoly_check_geodetic(LWPOLY *poly)
+{
+ int i = 0;
+ assert(poly);
+
+ for( i = 0; i < poly->nrings; i++ )
+ {
+ if( ptarray_check_geodetic(poly->rings[i]) == G_FALSE )
+ return G_FALSE;
+ }
+ return G_TRUE;
+}
+
+static int lwcollection_check_geodetic(LWCOLLECTION *col)
+{
+ int i = 0;
+ assert(col);
+
+ for( i = 0; i < col->ngeoms; i++ )
+ {
+ if( lwgeom_check_geodetic(col->geoms[i]) == G_FALSE )
+ return G_FALSE;
+ }
+ return G_TRUE;
+}
+
+int lwgeom_check_geodetic(const LWGEOM *geom)
+{
+ switch (TYPE_GETTYPE(geom->type))
+ {
+ case POINTTYPE:
+ return lwpoint_check_geodetic((LWPOINT *)geom);
+ case LINETYPE:
+ return lwline_check_geodetic((LWLINE *)geom);
+ case POLYGONTYPE:
+ return lwpoly_check_geodetic((LWPOLY *)geom);
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTIPOLYGONTYPE:
+ case COLLECTIONTYPE:
+ return lwcollection_check_geodetic((LWCOLLECTION *)geom);
+ default:
+ lwerror("unsupported input geometry type: %d", TYPE_GETTYPE(geom->type));
+ }
+ return G_FALSE;
+}
+
+
+/*
+** Count points in an LWGEOM.
+*/
+
+static int lwcollection_count_vertices(LWCOLLECTION *col)
+{
+ int i = 0;
+ int v = 0; /* vertices */
+ assert(col);
+ for ( i = 0; i < col->ngeoms; i++ )
+ {
+ v += lwgeom_count_vertices(col->geoms[i]);
+ }
+ return v;
+}
+
+static int lwpolygon_count_vertices(LWPOLY *poly)
+{
+ int i = 0;
+ int v = 0; /* vertices */
+ assert(poly);
+ for ( i = 0; i < poly->nrings; i ++ )
+ {
+ v += poly->rings[i]->npoints;
+ }
+ return v;
+}
+
+static int lwline_count_vertices(LWLINE *line)
+{
+ assert(line);
+ if ( ! line->points )
+ return 0;
+ return line->points->npoints;
+}
+
+static int lwpoint_count_vertices(LWPOINT *point)
+{
+ assert(point);
+ if ( ! point->point )
+ return 0;
+ return 1;
+}
+
+int lwgeom_count_vertices(LWGEOM *geom)
+{
+ int result = 0;
+ LWDEBUGF(4, "got type %d", TYPE_GETTYPE(geom->type));
+ switch (TYPE_GETTYPE(geom->type))
+ {
+ case POINTTYPE:
+ result = lwpoint_count_vertices((LWPOINT *)geom);;
+ break;
+ case LINETYPE:
+ result = lwline_count_vertices((LWLINE *)geom);
+ break;
+ case POLYGONTYPE:
+ result = lwpolygon_count_vertices((LWPOLY *)geom);
+ break;
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTIPOLYGONTYPE:
+ case COLLECTIONTYPE:
+ result = lwcollection_count_vertices((LWCOLLECTION *)geom);
+ break;
+ default:
+ lwerror("unsupported input geometry type: %d", TYPE_GETTYPE(geom->type));
+ break;
+ }
+ LWDEBUGF(3, "counted %d vertices", result);
+ return result;
+}
+
+int lwgeom_needs_bbox(LWGEOM *geom)
+{
+ assert(geom);
+ if( TYPE_GETTYPE(geom->type) == POINTTYPE )
+ {
+ return G_FALSE;
+ }
+ return G_TRUE;
+}
+
+
+
+
+
return NULL;
}
+
int
lwgeom_compute_box2d_p(LWGEOM *lwgeom, BOX2DFLOAT4 *buf)
{
return ret;
}
+
/* find length of this deserialized line */
size_t
lwline_serialize_size(LWLINE *line)
return ptarray_compute_box3d(point->point);
}
+
+
/*
* Convenience functions to hide the POINTARRAY
* TODO: obsolete this
size_t pasize;
uint32 npoints;
+ LWDEBUGF(4, "TYPE_GETZM(poly->type) == %d", TYPE_GETZM(poly->type));
+ LWDEBUGF(4, "TYPE_GETZM(pa->dims) == %d", TYPE_GETZM(pa->dims));
+
if ( TYPE_GETZM(poly->type) != TYPE_GETZM(pa->dims) )
lwerror("Dimensions mismatch in lwpoly");
return result;
}
+
/* find length of this serialized polygon */
size_t
lwgeom_size_poly(const uchar *serialized_poly)
{
"Unknown",
"Point",
- "Line",
+ "LineString",
"Polygon",
"MultiPoint",
- "MultiLine",
+ "MultiLineString",
"MultiPolygon",
"GeometryCollection",
"CircularString",
TYPE_SETZM(dims, hasz?1:0, hasm?1:0);
size = TYPE_NDIMS(dims)*npoints*sizeof(double);
- ptlist = (uchar *)lwalloc(size);
+ if( size )
+ ptlist = (uchar *)lwalloc(size);
+ else
+ ptlist = NULL;
+
pa = lwalloc(sizeof(POINTARRAY));
pa->dims = dims;
pa->serialized_pointlist = ptlist;
# **********************************************************************
-# * $Id: Makefile.in
+# * $Id$
# *
# * PostGIS - Spatial Types for PostgreSQL
# * http://postgis.refractions.net
MODULE_big=postgis-@POSTGIS_MAJOR_VERSION@.@POSTGIS_MINOR_VERSION@
# Files to be copied to the contrib/ directory
-DATA_built=postgis.sql uninstall_postgis.sql postgis_upgrade.sql
+DATA_built=postgis.sql uninstall_postgis.sql postgis_upgrade.sql geography.sql
DATA=../spatial_ref_sys.sql
# SQL objects (files requiring C pre-processing)
-SQL_OBJS=postgis.sql.in uninstall_postgis.sql.in
+SQL_OBJS=postgis.sql.in uninstall_postgis.sql.in geography.sql.in
# PostgreSQL objects
PG_OBJS=lwgeom_pg.o \
lwgeom_functions_lrs.o \
long_xact.o \
lwgeom_sqlmm.o \
- lwgeom_rtree.o
+ lwgeom_rtree.o \
+ geography_inout.o \
+ geography_gist.o
# Objects to build using PGXS
OBJS=$(PG_OBJS)
# Borrow the $libdir substitution from PGXS but customise by adding the version number
%.sql: %.sql.in
- sed 's,MODULE_PATHNAME,$$libdir/$*-@POSTGIS_MAJOR_VERSION@.@POSTGIS_MINOR_VERSION@,g' $< >$@
+ sed 's,MODULE_PATHNAME,$$libdir/postgis-@POSTGIS_MAJOR_VERSION@.@POSTGIS_MINOR_VERSION@,g' $< >$@
postgis_upgrade.sql: postgis.sql
$(PERL) ../utils/postgis_proc_upgrade.pl $< > $@
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+/* Useful functions for all GEOGRAPHY handlers. */
+GSERIALIZED* geography_serialize(LWGEOM *lwgeom);
+void geography_valid_typmod(LWGEOM *lwgeom, int32 typmod);
+void geography_valid_type(uchar type);
--- /dev/null
+---------------------------------------------------------------------------
+-- $Id$
+--
+-- PostGIS - Spatial Types for PostgreSQL
+-- Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+--
+-- This is free software; you can redistribute and/or modify it under
+-- the terms of the GNU General Public Licence. See the COPYING file.
+--
+---------------------------------------------------------------------------
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_typmod_in(cstring[])
+ RETURNS integer
+ AS 'MODULE_PATHNAME','geography_typmod_in'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_typmod_out(integer)
+ RETURNS cstring
+ AS 'MODULE_PATHNAME','geography_typmod_out'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_in(cstring, oid, integer)
+ RETURNS geography
+ AS 'MODULE_PATHNAME','geography_in'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_out(geography)
+ RETURNS cstring
+ AS 'MODULE_PATHNAME','geography_out'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE TYPE geography (
+ internallength = variable,
+ input = geography_in,
+ output = geography_out,
+ typmod_in = geography_typmod_in,
+ typmod_out = geography_typmod_out,
+ storage = main,
+ alignment = double
+);
+
+--
+-- GIDX type is used by the GiST index bindings.
+-- In/out functions are stubs, as all access should be internal.
+---
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION gidx_in(cstring)
+ RETURNS gidx
+ AS 'MODULE_PATHNAME','gidx_in'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION gidx_out(gidx)
+ RETURNS cstring
+ AS 'MODULE_PATHNAME','gidx_out'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE TYPE gidx (
+ internallength = variable,
+ input = gidx_in,
+ output = gidx_out,
+ storage = plain,
+ alignment = double
+);
+
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography(geography, integer, boolean)
+ RETURNS geography
+ AS 'MODULE_PATHNAME','geography_enforce_typmod'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE CAST (geography AS geography) WITH FUNCTION geography(geography, integer, boolean) AS IMPLICIT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION ST_AsText(geography)
+ RETURNS text
+ AS 'MODULE_PATHNAME','geography_as_text'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION ST_GeographyFromText(text)
+ RETURNS geography
+ AS 'MODULE_PATHNAME','geography_from_text'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION ST_AsBinary(geography)
+ RETURNS bytea
+ AS 'MODULE_PATHNAME','geography_as_binary'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION ST_GeographyFromBinary(bytea)
+ RETURNS geography
+ AS 'MODULE_PATHNAME','geography_from_binary'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_typmod_dims(integer)
+ RETURNS integer
+ AS 'MODULE_PATHNAME','geography_typmod_dims'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_typmod_srid(integer)
+ RETURNS integer
+ AS 'MODULE_PATHNAME','geography_typmod_srid'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_typmod_type(integer)
+ RETURNS text
+ AS 'MODULE_PATHNAME','geography_typmod_type'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE VIEW geography_columns AS
+ SELECT
+ current_database() AS f_table_catalog,
+ n.nspname AS f_table_schema,
+ c.relname AS f_table_name,
+ a.attname AS f_geography_column,
+ geography_typmod_dims(a.atttypmod) AS coord_dimension,
+ geography_typmod_srid(a.atttypmod) AS srid,
+ geography_typmod_type(a.atttypmod) AS type
+ FROM
+ pg_class c,
+ pg_attribute a,
+ pg_type t,
+ pg_namespace n
+ WHERE c.relkind IN('r','v')
+ AND t.typname = 'geography'
+ AND a.attisdropped = false
+ AND a.atttypid = t.oid
+ AND a.attrelid = c.oid
+ AND c.relnamespace = n.oid;
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography(geometry)
+ RETURNS geography
+ AS 'MODULE_PATHNAME','geography_from_geometry'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE CAST (geometry AS geography) WITH FUNCTION geography(geometry) AS IMPLICIT;
+
+-- ---------- ---------- ---------- ---------- ---------- ---------- ----------
+-- GiST Support Functions
+-- ---------- ---------- ---------- ---------- ---------- ---------- ----------
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_gist_consistent(internal,geometry,int4)
+ RETURNS bool
+ AS 'MODULE_PATHNAME' ,'geography_gist_consistent'
+ LANGUAGE 'C';
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_gist_compress(internal)
+ RETURNS internal
+ AS 'MODULE_PATHNAME','geography_gist_compress'
+ LANGUAGE 'C';
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_gist_penalty(internal,internal,internal)
+ RETURNS internal
+ AS 'MODULE_PATHNAME' ,'geography_gist_penalty'
+ LANGUAGE 'C';
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_gist_picksplit(internal, internal)
+ RETURNS internal
+ AS 'MODULE_PATHNAME' ,'geography_gist_picksplit'
+ LANGUAGE 'C';
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_gist_union(bytea, internal)
+ RETURNS internal
+ AS 'MODULE_PATHNAME' ,'geography_gist_union'
+ LANGUAGE 'C';
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_gist_same(box2d, box2d, internal)
+ RETURNS internal
+ AS 'MODULE_PATHNAME' ,'geography_gist_same'
+ LANGUAGE 'C';
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_gist_decompress(internal)
+ RETURNS internal
+ AS 'MODULE_PATHNAME' ,'geography_gist_decompress'
+ LANGUAGE 'C';
+
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION geography_overlaps(geography, geography)
+ RETURNS boolean
+ AS 'MODULE_PATHNAME' ,'geography_overlaps'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- Availability: 1.5.0
+CREATE OPERATOR && (
+ LEFTARG = geography, RIGHTARG = geography, PROCEDURE = geography_overlaps,
+ COMMUTATOR = '&&'
+-- ,RESTRICT = geography_gist_selectivity,
+-- JOIN = geography_gist_join_selectivity
+);
+
+
+-- Availability: 1.5.0
+CREATE OPERATOR CLASS gist_geography_ops
+ DEFAULT FOR TYPE geography USING GIST AS
+ STORAGE gidx,
+ OPERATOR 3 && ,
+-- OPERATOR 6 ~= ,
+-- OPERATOR 7 ~ ,
+-- OPERATOR 8 @ ,
+ FUNCTION 1 geography_gist_consistent (internal, geometry, int4),
+ FUNCTION 2 geography_gist_union (bytea, internal),
+ FUNCTION 3 geography_gist_compress (internal),
+ FUNCTION 4 geography_gist_decompress (internal),
+ FUNCTION 5 geography_gist_penalty (internal, internal, internal),
+ FUNCTION 6 geography_gist_picksplit (internal, internal),
+ FUNCTION 7 geography_gist_same (box2d, box2d, internal);
+
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+/*
+** R-Tree Bibliography
+**
+** [1] A. Guttman. R-tree: a dynamic index structure for spatial searching.
+** Proceedings of the ACM SIGMOD Conference, pp 47-57, June 1984.
+** [2] C.-H. Ang and T. C. Tan. New linear node splitting algorithm for
+** R-Trees. Advances in Spatial Databases - 5th International Symposium,
+** 1997
+** [3] N. Beckmann, H.-P. Kriegel, R. Schneider, B. Seeger. The R*tree: an
+** efficient and robust access method for points and rectangles.
+** Proceedings of the ACM SIGMOD Conference. June 1990.
+*/
+
+#include "postgres.h"
+#include "access/gist.h" /* For GiST */
+#include "access/itup.h"
+#include "access/skey.h"
+
+#include "../postgis_config.h"
+
+#include "libgeom.h" /* For standard geometry types. */
+#include "lwgeom_pg.h" /* For debugging macros. */
+#include "geography.h" /* For utility functions. */
+
+/*********************************************************************************
+** GIDX structure.
+**
+** This is an n-dimensional (practically, the
+** implementation is 2-4 dimensions) box used for index keys. The
+** coordinates are anonymous, so we can have any number of dimensions.
+** The sizeof a GIDX is 1 + 2 * ndims * sizeof(float).
+** The order of values in a GIDX is
+** xmin,xmax,ymin,ymax,zmin,zmax...
+*/
+typedef struct
+{
+ int32 varsize;
+ float c[1];
+} GIDX;
+
+
+/*
+** For some GiST support functions, it is more efficient to allocate
+** memory for our GIDX off the stack and cast that memory into a GIDX.
+** But, GIDX is variable length, what to do? We'll bake in the assumption
+** that no GIDX is more than 4-dimensional for now, and ensure that much
+** space is available.
+** 4 bytes varsize + 4 dimensions * 2 ordinates * 4 bytes float size = 36 bytes
+*/
+#define GIDX_MAX_SIZE 36
+
+/*
+** When is a node split not so good? If more than 90% of the entries
+** end up in one of the children.
+*/
+#define LIMIT_RATIO 0.1
+
+
+/*
+** GiST prototypes
+*/
+Datum geography_gist_consistent(PG_FUNCTION_ARGS);
+Datum geography_gist_compress(PG_FUNCTION_ARGS);
+Datum geography_gist_decompress(PG_FUNCTION_ARGS);
+Datum geography_gist_penalty(PG_FUNCTION_ARGS);
+Datum geography_gist_picksplit(PG_FUNCTION_ARGS);
+Datum geography_gist_union(PG_FUNCTION_ARGS);
+Datum geography_gist_same(PG_FUNCTION_ARGS);
+
+/*
+** Index key type stub prototypes
+*/
+Datum gidx_out(PG_FUNCTION_ARGS);
+Datum gidx_in(PG_FUNCTION_ARGS);
+
+/*
+** Operator prototypes
+*/
+Datum geography_overlaps(PG_FUNCTION_ARGS);
+
+
+/*********************************************************************************
+** GIDX support functions.
+**
+** We put the GIDX support here rather than libgeom because it is a specialized
+** type only used for indexing purposes. It also makes use of some PgSQL
+** infrastructure like the VARSIZE() macros.
+*/
+
+/*
+** Inline some of the most common and simplest utility functions.
+*/
+
+/* Returns number of dimensions for this GIDX */
+#define gidx_ndims(gidx) ((VARSIZE((gidx)) - VARHDRSZ) / (2 * sizeof(float)))
+/* Minimum accessor. */
+#define gidx_get_min(gidx, dimension) ((gidx)->c[2*(dimension)])
+/* Maximum accessor. */
+#define gidx_get_max(gidx, dimension) ((gidx)->c[2*(dimension)+1])
+/* Minimum setter. */
+#define gidx_set_min(gidx, dimension, value) ((gidx)->c[2*(dimension)] = (value))
+/* Maximum setter. */
+#define gidx_set_max(gidx, dimension, value) ((gidx)->c[2*(dimension)+1] = (value))
+/* Returns the size required to store a GIDX of requested dimension */
+#define gidx_size(dimensions) (sizeof(int32) + 2*(dimensions)*sizeof(float))
+
+
+/* Allocates a new GIDX on the heap of the requested dimensionality */
+static GIDX* gidx_new(int ndims)
+{
+ size_t size = gidx_size(ndims);
+ GIDX *g = (GIDX*)palloc(size);
+ POSTGIS_DEBUGF(5,"created new gidx of %d dimensions, size %d", ndims, (int)size);
+ SET_VARSIZE(g, size);
+ return g;
+}
+
+/* Generate human readable form for GIDX. */
+#if POSTGIS_DEBUG_LEVEL > 0
+static char* gidx_to_string(GIDX *a)
+{
+ char *str, *rv;
+ int i, ndims;
+
+ if( a == NULL )
+ return pstrdup("<NULLPTR>");
+
+ str = (char*)palloc(128);
+ rv = str;
+ ndims = gidx_ndims(a);
+
+ str += sprintf(str, "GIDX(");
+ for ( i = 0; i < ndims; i++ )
+ str += sprintf(str, " %.12g", gidx_get_min(a,i));
+ str += sprintf(str, ",");
+ for ( i = 0; i < ndims; i++ )
+ str += sprintf(str, " %.12g", gidx_get_max(a,i));
+ str += sprintf(str, " )");
+
+ return rv;
+}
+#endif
+
+/* Allocate a new copy of GIDX */
+static GIDX* gidx_copy(GIDX *b)
+{
+ GIDX *c = (GIDX*)palloc(VARSIZE(b));
+ POSTGIS_DEBUGF(5, "copied gidx (%p) to gidx (%p)", b, c);
+ memcpy((void*)c, (void*)b, VARSIZE(b));
+ return c;
+}
+
+/* Ensure all minimums are below maximums. */
+static inline void gidx_validate(GIDX *b)
+{
+ int i;
+ Assert(b);
+ POSTGIS_DEBUGF(5,"validating gidx (%s)", gidx_to_string(b));
+ for( i = 0; i < gidx_ndims(b); i++ )
+ {
+ if( gidx_get_min(b,i) > gidx_get_max(b,i) )
+ {
+ float tmp;
+ tmp = gidx_get_min(b,i);
+ gidx_set_min(b,i,gidx_get_max(b,i));
+ gidx_set_max(b,i,tmp);
+ }
+ }
+ return;
+}
+
+/* Enlarge b_union to contain b_new. If b_new contains more
+ dimensions than b_union, expand b_union to contain those dimensions. */
+static void gidx_merge(GIDX **b_union, GIDX *b_new)
+{
+ int i, dims_union, dims_new;
+ Assert(b_union);
+ Assert(*b_union);
+ Assert(b_new);
+
+ dims_union = gidx_ndims(*b_union);
+ dims_new = gidx_ndims(b_new);
+
+ POSTGIS_DEBUGF(4, "merging gidx (%s) into gidx (%s)", gidx_to_string(b_new), gidx_to_string(*b_union));
+
+ if( dims_new > dims_union )
+ {
+ POSTGIS_DEBUGF(5, "reallocating b_union from %d dims to %d dims", dims_union, dims_new);
+ *b_union = (GIDX*)repalloc(*b_union, gidx_size(dims_new));
+ SET_VARSIZE(*b_union, VARSIZE(b_new));
+ dims_union = dims_new;
+ }
+
+ for( i = 0; i < dims_new; i++ )
+ {
+ /* Adjust minimums */
+ gidx_set_min(*b_union, i, Min(gidx_get_min(*b_union,i),gidx_get_min(b_new,i)));
+ /* Adjust maximums */
+ gidx_set_max(*b_union, i, Max(gidx_get_max(*b_union,i),gidx_get_max(b_new,i)));
+ }
+
+ POSTGIS_DEBUGF(5, "merge complete (%s)", gidx_to_string(*b_union));
+ return;
+}
+
+/* Calculate the volume (in n-d units) of the GIDX */
+static float gidx_volume(GIDX *a)
+{
+ float result;
+ int i;
+ if ( a == NULL )
+ {
+ elog(ERROR, "gidx_volume received a null argument");
+ return 0.0;
+ }
+ result = gidx_get_max(a,0) - gidx_get_min(a,0);
+ for( i = 1; i < gidx_ndims(a); i++ )
+ result *= (gidx_get_max(a,i) - gidx_get_min(a,i));
+ POSTGIS_DEBUGF(5, "calculated volume of %s as %.12g", gidx_to_string(a), result);
+ return result;
+}
+
+/* Ensure the first argument has the higher dimensionality. */
+static void gidx_dimensionality_check(GIDX **a, GIDX **b)
+{
+ if( gidx_ndims(*a) < gidx_ndims(*b) )
+ {
+ GIDX *tmp = *b;
+ *b = *a;
+ *a = tmp;
+ }
+}
+
+/* Calculate the volume of the union of the boxes. Avoids creating an intermediate box. */
+static float gidx_union_volume(GIDX *a, GIDX *b)
+{
+ float result;
+ int i;
+ int ndims_a, ndims_b;
+
+ POSTGIS_DEBUG(5,"entered function");
+
+ if ( a == NULL && b == NULL )
+ {
+ elog(ERROR, "gidx_union_volume received two null arguments");
+ return 0.0;
+ }
+ if ( a == NULL )
+ return gidx_volume(b);
+
+ if ( b == NULL )
+ return gidx_volume(a);
+
+ /* Ensure 'a' has the most dimensions. */
+ gidx_dimensionality_check(&a, &b);
+
+ ndims_a = gidx_ndims(a);
+ ndims_b = gidx_ndims(b);
+
+ /* Initialize with maximal length of first dimension. */
+ result = Max(gidx_get_max(a,0),gidx_get_max(b,0)) - Min(gidx_get_min(a,0),gidx_get_min(b,0));
+
+ /* Multiply by maximal length of remaining dimensions. */
+ for( i = 1; i < ndims_b; i++ )
+ {
+ result *= (Max(gidx_get_max(a,i),gidx_get_max(b,i)) - Min(gidx_get_min(a,i),gidx_get_min(b,i)));
+ }
+
+ /* Add in dimensions of higher dimensional box. */
+ for( i = ndims_b; i < ndims_a; i++ )
+ {
+ result *= (gidx_get_max(a,i) - gidx_get_min(a,i));
+ }
+
+ POSTGIS_DEBUGF(5, "volume( %s union %s ) = %.12g", gidx_to_string(a), gidx_to_string(b), result);
+
+ return result;
+}
+
+/* Calculate the volume of the intersection of the boxes. */
+static float gidx_inter_volume(GIDX *a, GIDX *b)
+{
+ int i;
+ float result;
+
+ POSTGIS_DEBUG(5,"entered function");
+
+ if ( a == NULL || b == NULL )
+ {
+ elog(ERROR, "gidx_inter_volume received a null argument");
+ return 0.0;
+ }
+
+ /* Ensure 'a' has the most dimensions. */
+ gidx_dimensionality_check(&a, &b);
+
+ /* Initialize with minimal length of first dimension. */
+ result = Min(gidx_get_max(a,0),gidx_get_max(b,0)) - Max(gidx_get_min(a,0),gidx_get_min(b,0));
+
+ /* If they are disjoint (max < min) then return zero. */
+ if ( result < 0.0 ) return 0.0;
+
+ /* Continue for remaining dimensions. */
+ for( i = 1; i < gidx_ndims(b); i++ )
+ {
+ float width = Min(gidx_get_max(a,i),gidx_get_max(b,i)) - Max(gidx_get_min(a,i),gidx_get_min(b,i));
+ if ( width < 0.0 ) return 0.0;
+ /* Multiply by minimal length of remaining dimensions. */
+ result *= width;
+ }
+ POSTGIS_DEBUGF(5, "volume( %s intersection %s ) = %.12g", gidx_to_string(a), gidx_to_string(b), result);
+ return result;
+}
+
+/* Convert a double-based GBOX into a float-based GIDX,
+ ensuring the float box is larger than the double box */
+static int gidx_from_gbox_p(GBOX *box, GIDX *a)
+{
+ int ndims;
+
+ if ( box == NULL )
+ {
+ elog(ERROR, "gidx_from_gbox received a null argument");
+ return G_FAILURE;
+ }
+
+ ndims = (FLAGS_GET_GEODETIC(box->flags) ? 3 : FLAGS_NDIMS(box->flags));
+ SET_VARSIZE(a, VARHDRSZ + ndims * 2 * sizeof(float));
+
+ gidx_set_min(a,0,nextDown_f(box->xmin));
+ gidx_set_max(a,0,nextUp_f(box->xmax));
+ gidx_set_min(a,1,nextDown_f(box->ymin));
+ gidx_set_max(a,1,nextUp_f(box->ymax));
+
+ /* Geodetic indexes are always 3d, geocentric x/y/z */
+ if ( FLAGS_GET_GEODETIC(box->flags) )
+ {
+ gidx_set_min(a,2,nextDown_f(box->zmin));
+ gidx_set_max(a,2,nextUp_f(box->zmax));
+ }
+ else
+ {
+ /* Cartesian with Z implies Z is third dimension */
+ if ( FLAGS_GET_Z(box->flags) )
+ {
+ gidx_set_min(a,2,nextDown_f(box->zmin));
+ gidx_set_max(a,2,nextUp_f(box->zmax));
+ if ( FLAGS_GET_M(box->flags) )
+ {
+ gidx_set_min(a,3,nextDown_f(box->mmin));
+ gidx_set_max(a,3,nextUp_f(box->mmax));
+ }
+ }
+ /* Unless there's no Z, in which case M is third dimension */
+ else if ( FLAGS_GET_M(box->flags) )
+ {
+ gidx_set_min(a,2,nextDown_f(box->mmin));
+ gidx_set_max(a,2,nextUp_f(box->mmax));
+ }
+ }
+
+ POSTGIS_DEBUGF(5, "converted %s to %s", gbox_to_string(box), gidx_to_string(a));
+
+ return G_SUCCESS;
+}
+
+static GIDX* gidx_from_gbox(GBOX *box)
+{
+ int ndims;
+ GIDX *a;
+
+ if ( box == NULL )
+ {
+ elog(ERROR, "gidx_from_gbox received a null argument");
+ return NULL;
+ }
+ ndims = (FLAGS_GET_GEODETIC(box->flags) ? 3 : FLAGS_NDIMS(box->flags));
+ a = gidx_new(ndims);
+ gidx_from_gbox_p(box, a);
+ return a;
+}
+
+
+/*
+** Overlapping GIDX box test.
+**
+** Box(A) Overlap Box(B) IFF (pt(a)LL < pt(B)UR) && (pt(b)LL < pt(a)UR)
+*/
+static bool gidx_overlaps(GIDX *a, GIDX *b)
+{
+ int i;
+ int ndims_b;
+ POSTGIS_DEBUG(5, "entered function");
+ if( (a == NULL) || (b == NULL) ) return FALSE;
+
+ /* Ensure 'a' has the most dimensions. */
+ gidx_dimensionality_check(&a, &b);
+
+ ndims_b = gidx_ndims(b);
+
+ /* compare within the dimensions of (b) */
+ for( i = 0; i < ndims_b; i++ )
+ {
+ if( gidx_get_min(a,i) > gidx_get_max(b,i) )
+ return FALSE;
+ if( gidx_get_min(b,i) > gidx_get_max(a,i) )
+ return FALSE;
+ }
+
+ /* compare to zero those dimensions in (a) absent in (b) */
+ for( i = ndims_b; i < gidx_ndims(a); i++ )
+ {
+ if( gidx_get_min(a,i) > 0.0 )
+ return FALSE;
+ if( gidx_get_max(a,i) < 0.0 )
+ return FALSE;
+ }
+ return TRUE;
+}
+
+/*
+** Containment GIDX test.
+**
+** Box(A) CONTAINS Box(B) IFF (pt(A)LL < pt(B)LL) && (pt(A)UR > pt(B)UR)
+*/
+static bool gidx_contains(GIDX *a, GIDX *b)
+{
+ int i, dims_a, dims_b;
+
+ POSTGIS_DEBUG(5, "entered function");
+
+ if( (a == NULL) || (b == NULL) ) return FALSE;
+
+ dims_a = gidx_ndims(a);
+ dims_b = gidx_ndims(b);
+
+ if( dims_a < dims_b )
+ {
+ /*
+ ** If (b) is of higher dimensionality than (a) it can only be contained
+ ** if those higher dimensions are zeroes.
+ */
+ for (i = dims_a; i < dims_b; i++)
+ {
+ if( gidx_get_min(b,i) != 0 )
+ return FALSE;
+ if( gidx_get_max(b,i) != 0 )
+ return FALSE;
+ }
+ }
+
+ /* Excess dimensions of (a), don't matter, it just has to contain (b) in (b)'s dimensions */
+ for (i = 0; i < Min(dims_a, dims_b); i++)
+ {
+ if ( gidx_get_min(a,i) > gidx_get_min(b,i) )
+ return FALSE;
+ if ( gidx_get_max(a,i) < gidx_get_max(b,i) )
+ return FALSE;
+ }
+
+ return TRUE;
+}
+/*
+** Equality GIDX test.
+**
+** Box(A) EQUALS Box(B) IFF (pt(A)LL == pt(B)LL) && (pt(A)UR == pt(B)UR)
+*/
+static bool gidx_equals(GIDX *a, GIDX *b)
+{
+ int i;
+
+ POSTGIS_DEBUG(5, "entered function");
+
+ if( (a == NULL) && (b == NULL) ) return TRUE;
+ if( (a == NULL) || (b == NULL) ) return FALSE;
+
+ /* Ensure 'a' has the most dimensions. */
+ gidx_dimensionality_check(&a, &b);
+
+ /* For all shared dimensions min(a) == min(b), max(a) == max(b) */
+ for (i = 0; i < gidx_ndims(b); i++)
+ {
+ if ( gidx_get_min(a,i) != gidx_get_min(b,i) )
+ return FALSE;
+ if ( gidx_get_max(a,i) != gidx_get_max(b,i) )
+ return FALSE;
+ }
+ /* For all unshared dimensions min(a) == 0.0, max(a) == 0.0 */
+ for (i = gidx_ndims(b); i < gidx_ndims(a); i++)
+ {
+ if ( gidx_get_min(a,i) != 0.0 )
+ return FALSE;
+ if ( gidx_get_max(a,i) != 0.0 )
+ return FALSE;
+ }
+ return TRUE;
+}
+
+/*
+** Peak into a geography (gserialized) datum to find the bounding box. If the
+** box is there, copy it out and return it. If not, calculate the box from the
+** full geography and return the box based on that. If no box is available,
+** return G_FAILURE, otherwise G_SUCCESS.
+*/
+static int geography_datum_gidx(Datum geography_datum, GIDX *gidx)
+{
+ GSERIALIZED *gpart;
+ int result = G_SUCCESS;
+
+ POSTGIS_DEBUG(4, "entered function");
+
+ /*
+ ** The most info we need is the 8 bytes of serialized header plus the 24 bytes
+ ** of floats necessary to hold the 6 floats of the geocentric index
+ ** bounding box, so 32 bytes.
+ */
+ gpart = (GSERIALIZED*)PG_DETOAST_DATUM_SLICE(geography_datum, 0, 32);
+
+ POSTGIS_DEBUGF(4, "got flags %d", gpart->flags);
+
+ if ( FLAGS_GET_BBOX(gpart->flags) && FLAGS_GET_GEODETIC(gpart->flags) )
+ {
+ const size_t size = 2 * 3 * sizeof(float);
+ POSTGIS_DEBUG(4, "copying box out of serialization");
+ memcpy(gidx->c, gpart->data, size);
+ SET_VARSIZE(gidx, VARHDRSZ + size);
+ }
+ else
+ {
+ GSERIALIZED *g = (GSERIALIZED*)PG_DETOAST_DATUM(geography_datum);
+ GBOX gbox;
+ POSTGIS_DEBUG(4, "calculating new box from scratch");
+ if( gserialized_calculate_gbox_geocentric_p(g, &gbox) == G_FAILURE )
+ {
+ POSTGIS_DEBUG(4, "calculated null bbox, returning null");
+ return G_FAILURE;
+ }
+ result = gidx_from_gbox_p(&gbox, gidx);
+ }
+ if( result == G_SUCCESS )
+ {
+ POSTGIS_DEBUGF(4, "got gidx %s", gidx_to_string(gidx));
+ }
+
+ return result;
+}
+
+/***********************************************************************
+* GiST Support Functions
+*/
+
+/*
+** GiST support function. Based on two geographies return true if
+** they overlap and false otherwise.
+*/
+PG_FUNCTION_INFO_V1(geography_overlaps);
+Datum geography_overlaps(PG_FUNCTION_ARGS)
+{
+ /* Put aside some stack memory and use it for GIDX pointers. */
+ char gboxmem1[GIDX_MAX_SIZE];
+ char gboxmem2[GIDX_MAX_SIZE];
+ GIDX *gbox1 = (GIDX*)gboxmem1;
+ GIDX *gbox2 = (GIDX*)gboxmem2;
+
+ geography_datum_gidx(PG_GETARG_DATUM(0), gbox1);
+ geography_datum_gidx(PG_GETARG_DATUM(1), gbox2);
+
+ if ( gidx_overlaps(gbox1, gbox2) == G_TRUE )
+ PG_RETURN_BOOL(TRUE);
+
+ PG_RETURN_BOOL(FALSE);
+}
+
+/*
+** GiST support function. Given a geography, return a "compressed"
+** version. In this case, we convert the geography into a geocentric
+** bounding box. If the geography already has the box embedded in it
+** we pull that out and hand it back.
+*/
+PG_FUNCTION_INFO_V1(geography_gist_compress);
+Datum geography_gist_compress(PG_FUNCTION_ARGS)
+{
+ GISTENTRY *entry_in = (GISTENTRY*)PG_GETARG_POINTER(0);
+ GISTENTRY *entry_out = NULL;
+ char gidxmem[GIDX_MAX_SIZE];
+ GIDX *bbox_out = (GIDX*)gidxmem;
+ int result = G_SUCCESS;
+
+ int i;
+
+ POSTGIS_DEBUG(4, "[GIST] 'compress' function called");
+
+ /*
+ ** Not a leaf key? There's nothing to do.
+ ** Return the input unchanged.
+ */
+ if( ! entry_in->leafkey )
+ {
+ POSTGIS_DEBUG(4, "[GIST] non-leafkey entry, returning input unaltered");
+ PG_RETURN_POINTER(entry_in);
+ }
+
+ POSTGIS_DEBUG(4, "[GIST] processing leafkey input");
+ entry_out = palloc(sizeof(GISTENTRY));
+
+ /*
+ ** Null key? Make a copy of the input entry and
+ ** return.
+ */
+ if ( DatumGetPointer(entry_in->key) == NULL )
+ {
+ POSTGIS_DEBUG(4, "[GIST] leafkey is null");
+ gistentryinit(*entry_out, (Datum) 0, entry_in->rel,
+ entry_in->page, entry_in->offset, FALSE);
+ POSTGIS_DEBUG(4, "[GIST] returning copy of input");
+ PG_RETURN_POINTER(entry_out);
+ }
+
+ /* Extract our index key from the GiST entry. */
+ result = geography_datum_gidx(entry_in->key, bbox_out);
+
+ /* Is the bounding box valid (non-empty, non-infinite)? If not, return input uncompressed. */
+ if ( result == G_FAILURE )
+ {
+ POSTGIS_DEBUG(4, "[GIST] empty geometry!");
+ PG_RETURN_POINTER(entry_in);
+ }
+
+ POSTGIS_DEBUGF(4, "[GIST] got entry_in->key: %s", gidx_to_string(bbox_out));
+
+ /* Check all the dimensions for finite values */
+ for( i = 0; i < gidx_ndims(bbox_out); i++ )
+ {
+ if( ! finite(gidx_get_max(bbox_out, i)) || ! finite(gidx_get_min(bbox_out, i)) )
+ {
+ POSTGIS_DEBUG(4, "[GIST] infinite geometry!");
+ PG_RETURN_POINTER(entry_in);
+ }
+ }
+
+ /* Enure bounding box has minimums below maximums. */
+ gidx_validate(bbox_out);
+
+ /* Prepare GISTENTRY for return. */
+ gistentryinit(*entry_out, PointerGetDatum(gidx_copy(bbox_out)),
+ entry_in->rel, entry_in->page, entry_in->offset, FALSE);
+
+ /* Return GISTENTRY. */
+ POSTGIS_DEBUG(4, "[GIST] 'compress' function complete");
+ PG_RETURN_POINTER(entry_out);
+}
+
+
+/*
+** GiST support function.
+** Decompress an entry. Unused for geography, so we return.
+*/
+PG_FUNCTION_INFO_V1(geography_gist_decompress);
+Datum geography_gist_decompress(PG_FUNCTION_ARGS)
+{
+ POSTGIS_DEBUG(5, "[GIST] 'decompress' function called");
+ /* We don't decompress. Just return the input. */
+ PG_RETURN_POINTER(PG_GETARG_POINTER(0));
+}
+
+/*
+** GiST support function. Called from geography_gist_consistent below.
+*/
+static inline bool geography_gist_consistent_leaf(GIDX *key, GIDX *query, StrategyNumber strategy)
+{
+ bool retval;
+
+ POSTGIS_DEBUGF(4, "[GIST] leaf consistent, strategy: %d", strategy);
+
+ switch (strategy)
+ {
+ case RTOverlapStrategyNumber:
+ retval = (bool) gidx_overlaps(key, query);
+ break;
+ case RTSameStrategyNumber:
+ retval = (bool) gidx_equals(key, query);
+ break;
+ case RTContainsStrategyNumber:
+ case RTOldContainsStrategyNumber:
+ retval = (bool) gidx_contains(key, query);
+ break;
+ case RTContainedByStrategyNumber:
+ case RTOldContainedByStrategyNumber:
+ retval = (bool) gidx_contains(query, key);
+ break;
+ default:
+ retval = FALSE;
+ }
+ return (retval);
+}
+
+/*
+** GiST support function. Called from geography_gist_consistent below.
+*/
+static inline bool geography_gist_consistent_internal(GIDX *key, GIDX *query, StrategyNumber strategy)
+{
+ bool retval;
+
+ POSTGIS_DEBUGF(4, "[GIST] internal consistent, strategy: %d", strategy);
+
+ switch (strategy)
+ {
+ case RTOverlapStrategyNumber:
+ retval = (bool) gidx_overlaps(key, query);
+ break;
+ case RTSameStrategyNumber:
+ case RTContainsStrategyNumber:
+ case RTOldContainsStrategyNumber:
+ retval = (bool) gidx_contains(key, query);
+ break;
+ case RTContainedByStrategyNumber:
+ case RTOldContainedByStrategyNumber:
+ retval = (bool) gidx_overlaps(key, query);
+ break;
+ default:
+ retval = FALSE;
+ }
+ return (retval);
+}
+
+/*
+** GiST support function. Take in a query and an entry and see what the
+** relationship is, based on the query strategy.
+*/
+PG_FUNCTION_INFO_V1(geography_gist_consistent);
+Datum geography_gist_consistent(PG_FUNCTION_ARGS)
+{
+ GISTENTRY *entry = (GISTENTRY*) PG_GETARG_POINTER(0);
+ StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
+ bool result;
+ char gidxmem[GIDX_MAX_SIZE];
+ GIDX *query_gbox_index = (GIDX*)gidxmem;
+
+#if POSTGIS_PGSQL_VERSION >= 84
+ /* PostgreSQL 8.4 and later require the RECHECK flag to be set here,
+ rather than being supplied as part of the operator class definition */
+ bool *recheck = (bool *) PG_GETARG_POINTER(4);
+
+ /* We set recheck to false to avoid repeatedly pulling every "possibly matched" geometry
+ out during index scans. For cases when the geometries are large, rechecking
+ can make things twice as slow. */
+ *recheck = false;
+#endif
+
+ POSTGIS_DEBUG(4, "[GIST] 'consistent' function called");
+
+ /* Quick sanity check on query argument. */
+ if ( DatumGetPointer(PG_GETARG_DATUM(1)) == NULL )
+ {
+ POSTGIS_DEBUG(4, "[GIST] null query pointer (!?!), returning false");
+ PG_RETURN_BOOL(FALSE); /* NULL query! This is screwy! */
+ }
+
+ /* Quick sanity check on entry key. */
+ if ( DatumGetPointer(entry->key) == NULL )
+ {
+ POSTGIS_DEBUG(4, "[GIST] null index entry, returning false");
+ PG_RETURN_BOOL(FALSE); /* NULL entry! */
+ }
+
+ /* Null box should never make this far. */
+ if ( geography_datum_gidx(PG_GETARG_DATUM(1), query_gbox_index) == G_FAILURE )
+ {
+ POSTGIS_DEBUG(4, "[GIST] null query_gbox_index!");
+ PG_RETURN_BOOL(FALSE);
+ }
+
+ /* Treat leaf node tests different from internal nodes */
+ if (GIST_LEAF(entry))
+ {
+ result = geography_gist_consistent_leaf(
+ (GIDX*)DatumGetPointer(entry->key),
+ query_gbox_index, strategy);
+ }
+ else
+ {
+ result = geography_gist_consistent_internal(
+ (GIDX*)DatumGetPointer(entry->key),
+ query_gbox_index, strategy);
+ }
+
+ PG_RETURN_BOOL(result);
+}
+
+
+/*
+** GiST support function. Calculate the "penalty" cost of adding this entry into an existing entry.
+** Calculate the change in volume of the old entry once the new entry is added.
+** TODO: Re-evaluate this in light of R*Tree penalty approaches.
+*/
+PG_FUNCTION_INFO_V1(geography_gist_penalty);
+Datum geography_gist_penalty(PG_FUNCTION_ARGS)
+{
+ GISTENTRY *origentry = (GISTENTRY*) PG_GETARG_POINTER(0);
+ GISTENTRY *newentry = (GISTENTRY*) PG_GETARG_POINTER(1);
+ float *result = (float*) PG_GETARG_POINTER(2);
+ GIDX *gbox_index_orig, *gbox_index_new;
+ float size_union, size_orig;
+
+ POSTGIS_DEBUG(4, "[GIST] 'penalty' function called");
+
+ gbox_index_orig = (GIDX*)DatumGetPointer(origentry->key);
+ gbox_index_new = (GIDX*)DatumGetPointer(newentry->key);
+
+ /* Drop out if we're dealing with null inputs. Shouldn't happen. */
+ if ( (gbox_index_orig == NULL) && (gbox_index_new == NULL) )
+ {
+ POSTGIS_DEBUG(4, "[GIST] both inputs NULL! returning penalty of zero");
+ *result = 0.0;
+ PG_RETURN_FLOAT8(*result);
+ }
+
+ /* Calculate the size difference of the boxes (volume difference in this case). */
+ size_union = gidx_union_volume(gbox_index_orig, gbox_index_new);
+ size_orig = gidx_volume(gbox_index_orig);
+ *result = size_union - size_orig;
+ POSTGIS_DEBUGF(4, "[GIST] union size (%f), original size (%f), penalty (%f)", size_union, size_orig, *result);
+ PG_RETURN_POINTER(result);
+}
+
+/*
+** GiST support function. Merge all the boxes in a page into one master box.
+*/
+PG_FUNCTION_INFO_V1(geography_gist_union);
+Datum geography_gist_union(PG_FUNCTION_ARGS)
+{
+ GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
+ int *sizep = (int *) PG_GETARG_POINTER(1); /* Size of the return value */
+ int numranges, i;
+ GIDX *box_cur, *box_union;
+
+ POSTGIS_DEBUG(4, "[GIST] 'union' function called");
+
+ numranges = entryvec->n;
+
+ box_cur = (GIDX*) DatumGetPointer(entryvec->vector[0].key);
+
+ box_union = gidx_copy(box_cur);
+
+ for ( i = 1; i < numranges; i++ )
+ {
+ box_cur = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
+ gidx_merge(&box_union, box_cur);
+ }
+
+ *sizep = VARSIZE(box_union);
+
+ POSTGIS_DEBUGF(4, "[GIST] union called, numranges(%i), pageunion %s", numranges, gidx_to_string(box_union));
+
+ PG_RETURN_POINTER(box_union);
+
+}
+
+/*
+** GiST support function. Test equality of keys.
+*/
+PG_FUNCTION_INFO_V1(geography_gist_same);
+Datum geography_gist_same(PG_FUNCTION_ARGS)
+{
+ GIDX *b1 = (GIDX*)PG_GETARG_POINTER(0);
+ GIDX *b2 = (GIDX*)PG_GETARG_POINTER(1);
+ bool *result = (bool*)PG_GETARG_POINTER(2);
+
+ POSTGIS_DEBUG(4, "[GIST] 'same' function called");
+
+ *result = gidx_equals(b1, b2);
+
+ PG_RETURN_POINTER(result);
+}
+
+
+/*
+** Utility function to add entries to the axis partition lists in the
+** picksplit function.
+*/
+static void geography_gist_picksplit_addlist(OffsetNumber *list, GIDX **box_union, GIDX *box_current, int *pos, int num)
+{
+ if( *pos )
+ gidx_merge(box_union, box_current);
+ else
+ memcpy((void*)(*box_union), (void*)box_current, VARSIZE(box_current));
+ list[*pos] = num;
+ (*pos)++;
+}
+
+/*
+** Utility function check whether the number of entries two halves of the
+** space constitute a "bad ratio" (poor balance).
+*/
+static int geography_gist_picksplit_badratio(int x, int y)
+{
+ POSTGIS_DEBUGF(4, "[GIST] checking split ratio (%d, %d)", x, y);
+ if( (y == 0) || (((float)x / (float)y) < LIMIT_RATIO) ||
+ (x == 0) || (((float)y / (float)x) < LIMIT_RATIO) )
+ return TRUE;
+
+ return FALSE;
+}
+
+/*
+** Where the picksplit algorithm cannot find any basis for splitting one way
+** or another, we simply split the overflowing node in half.
+*/
+static void geography_gist_picksplit_fallback(GistEntryVector *entryvec, GIST_SPLITVEC *v)
+{
+ OffsetNumber i, maxoff;
+ GIDX *unionL = NULL;
+ GIDX *unionR = NULL;
+ int nbytes;
+
+ POSTGIS_DEBUG(4, "[GIST] in fallback picksplit function");
+
+ maxoff = entryvec->n - 1;
+
+ nbytes = (maxoff + 2) * sizeof(OffsetNumber);
+ v->spl_left = (OffsetNumber*) palloc(nbytes);
+ v->spl_right = (OffsetNumber*) palloc(nbytes);
+ v->spl_nleft = v->spl_nright = 0;
+
+ for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
+ {
+ GIDX *cur = (GIDX*)DatumGetPointer(entryvec->vector[i].key);
+
+ if (i <= (maxoff - FirstOffsetNumber + 1) / 2)
+ {
+ v->spl_left[v->spl_nleft] = i;
+ if (unionL == NULL)
+ {
+ unionL = gidx_copy(cur);
+ }
+ else
+ {
+ gidx_merge(&unionL, cur);
+ }
+ v->spl_nleft++;
+ }
+ else
+ {
+ v->spl_right[v->spl_nright] = i;
+ if (unionR == NULL)
+ {
+ unionR = gidx_copy(cur);
+ }
+ else
+ {
+ gidx_merge(&unionR, cur);
+ }
+ v->spl_nright++;
+ }
+ }
+
+ if (v->spl_ldatum_exists)
+ gidx_merge(&unionL, (GIDX*)DatumGetPointer(v->spl_ldatum));
+
+ v->spl_ldatum = PointerGetDatum(unionL);
+
+ if (v->spl_rdatum_exists)
+ gidx_merge(&unionR, (GIDX*)DatumGetPointer(v->spl_rdatum));
+
+ v->spl_rdatum = PointerGetDatum(unionR);
+ v->spl_ldatum_exists = v->spl_rdatum_exists = false;
+}
+
+
+
+static void geography_gist_picksplit_constructsplit(GIST_SPLITVEC *v, OffsetNumber *list1, int nlist1, GIDX **union1, OffsetNumber *list2, int nlist2, GIDX **union2)
+{
+ bool firstToLeft = true;
+
+ POSTGIS_DEBUG(4, "[GIST] picksplit in constructsplit function");
+
+ if (v->spl_ldatum_exists || v->spl_rdatum_exists)
+ {
+ if (v->spl_ldatum_exists && v->spl_rdatum_exists)
+ {
+ GIDX *LRl = gidx_copy(*union1);
+ GIDX *LRr = gidx_copy(*union2);
+ GIDX *RLl = gidx_copy(*union2);
+ GIDX *RLr = gidx_copy(*union1);
+ double sizeLR, sizeRL;
+
+ gidx_merge(&LRl, (GIDX*)DatumGetPointer(v->spl_ldatum));
+ gidx_merge(&LRr, (GIDX*)DatumGetPointer(v->spl_rdatum));
+ gidx_merge(&RLl, (GIDX*)DatumGetPointer(v->spl_ldatum));
+ gidx_merge(&RLr, (GIDX*)DatumGetPointer(v->spl_rdatum));
+
+ sizeLR = gidx_inter_volume(LRl,LRr);
+ sizeRL = gidx_inter_volume(RLl,RLr);
+
+ if (sizeLR > sizeRL)
+ firstToLeft = false;
+
+ }
+ else
+ {
+ float p1, p2;
+ GISTENTRY oldUnion, addon;
+
+ gistentryinit(oldUnion, (v->spl_ldatum_exists) ? v->spl_ldatum : v->spl_rdatum,
+ NULL, NULL, InvalidOffsetNumber, FALSE);
+
+ gistentryinit(addon, PointerGetDatum(*union1), NULL, NULL, InvalidOffsetNumber, FALSE);
+ DirectFunctionCall3(geography_gist_penalty, PointerGetDatum(&oldUnion), PointerGetDatum(&addon), PointerGetDatum(&p1));
+ gistentryinit(addon, PointerGetDatum(*union2), NULL, NULL, InvalidOffsetNumber, FALSE);
+ DirectFunctionCall3(geography_gist_penalty, PointerGetDatum(&oldUnion), PointerGetDatum(&addon), PointerGetDatum(&p2));
+
+ if ((v->spl_ldatum_exists && p1 > p2) || (v->spl_rdatum_exists && p1 < p2))
+ firstToLeft = false;
+ }
+ }
+
+ if (firstToLeft)
+ {
+ v->spl_left = list1;
+ v->spl_right = list2;
+ v->spl_nleft = nlist1;
+ v->spl_nright = nlist2;
+ if (v->spl_ldatum_exists)
+ gidx_merge(union1, (GIDX*)DatumGetPointer(v->spl_ldatum));
+ v->spl_ldatum = PointerGetDatum(*union1);
+ if (v->spl_rdatum_exists)
+ gidx_merge(union2, (GIDX*)DatumGetPointer(v->spl_rdatum));
+ v->spl_rdatum = PointerGetDatum(*union2);
+ }
+ else
+ {
+ v->spl_left = list2;
+ v->spl_right = list1;
+ v->spl_nleft = nlist2;
+ v->spl_nright = nlist1;
+ if (v->spl_ldatum_exists)
+ gidx_merge(union2, (GIDX*)DatumGetPointer(v->spl_ldatum));
+ v->spl_ldatum = PointerGetDatum(*union2);
+ if (v->spl_rdatum_exists)
+ gidx_merge(union1, (GIDX*)DatumGetPointer(v->spl_rdatum));
+ v->spl_rdatum = PointerGetDatum(*union1);
+ }
+
+ v->spl_ldatum_exists = v->spl_rdatum_exists = false;
+}
+
+static bool geography_gist_picksplit_badratios(int *pos, int dims)
+{
+ int i;
+ for ( i = 0; i < dims; i++ )
+ {
+ if( geography_gist_picksplit_badratio(pos[2*i],pos[2*i+1]) == FALSE )
+ return FALSE;
+ }
+ return TRUE;
+}
+
+#define BELOW(d) (2*(d))
+#define ABOVE(d) ((2*(d))+1)
+
+/*
+** GiST support function. Split an overflowing node into two new nodes.
+** Uses linear algorithm from Ang & Tan [2], dividing node extent into
+** four quadrants, and splitting along the axis that most evenly distributes
+** entries between the new nodes.
+** TODO: Re-evaluate this in light of R*Tree picksplit approaches.
+*/
+PG_FUNCTION_INFO_V1(geography_gist_picksplit);
+Datum geography_gist_picksplit(PG_FUNCTION_ARGS)
+{
+
+ GistEntryVector *entryvec = (GistEntryVector*) PG_GETARG_POINTER(0);
+
+ GIST_SPLITVEC *v = (GIST_SPLITVEC*) PG_GETARG_POINTER(1);
+ OffsetNumber i;
+ /* One union box for each half of the space. */
+ GIDX **box_union;
+ /* One offset number list for each half of the space. */
+ OffsetNumber **list;
+ /* One position index for each half of the space. */
+ int *pos;
+ GIDX *box_pageunion;
+ GIDX *box_current;
+ int direction = -1;
+ bool all_entries_equal = true;
+ OffsetNumber max_offset;
+ int nbytes, ndims_pageunion, d;
+ int posmax = -1;
+
+ POSTGIS_DEBUG(4, "[GIST] 'picksplit' function called");
+
+ /*
+ ** First calculate the bounding box and maximum number of dimensions in this page.
+ */
+
+ max_offset = entryvec->n - 1;
+ box_current = (GIDX*) DatumGetPointer(entryvec->vector[FirstOffsetNumber].key);
+ POSTGIS_DEBUGF(5, "gidx minx %f", box_current->c[0]);
+ POSTGIS_DEBUGF(5, "gidx maxx %f", box_current->c[1]);
+ box_pageunion = gidx_copy(box_current);
+
+ POSTGIS_DEBUGF(4, "[GIST] box_current (%s)", gidx_to_string(box_current));
+
+ /* Calculate the containing box (box_pageunion) for the whole page we are going to split. */
+ for ( i = OffsetNumberNext(FirstOffsetNumber); i <= max_offset; i = OffsetNumberNext(i) )
+ {
+ box_current = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
+
+ if ( all_entries_equal == true && ! gidx_equals (box_pageunion, box_current) )
+ all_entries_equal = false;
+
+ gidx_merge( &box_pageunion, box_current );
+ }
+
+ POSTGIS_DEBUGF(4, "[GIST] box_pageunion (%s)", gidx_to_string(box_pageunion));
+
+ /* Every box in the page is the same! So, we split and just put half the boxes in each child. */
+ if ( all_entries_equal )
+ {
+ POSTGIS_DEBUG(4, "[GIST] picksplit finds all entries equal!");
+ geography_gist_picksplit_fallback(entryvec, v);
+ PG_RETURN_POINTER(v);
+ }
+
+ /* Initialize memory structures. */
+ nbytes = (max_offset + 2) * sizeof(OffsetNumber);
+ ndims_pageunion = gidx_ndims(box_pageunion);
+ pos = palloc(2*ndims_pageunion * sizeof(int));
+ list = palloc(2*ndims_pageunion * sizeof(OffsetNumber*));
+ box_union = palloc(2*ndims_pageunion * sizeof(GIDX*));
+ for( d = 0; d < ndims_pageunion; d++ )
+ {
+ list[BELOW(d)] = (OffsetNumber*) palloc(nbytes);
+ list[ABOVE(d)] = (OffsetNumber*) palloc(nbytes);
+ box_union[BELOW(d)] = gidx_new(ndims_pageunion);
+ box_union[ABOVE(d)] = gidx_new(ndims_pageunion);
+ pos[BELOW(d)] = 0;
+ pos[ABOVE(d)] = 0;
+ }
+
+ /*
+ ** Assign each entry in the node to the volume partitions it belongs to,
+ ** such as "above the x/y plane, left of the y/z plane, below the x/z plane".
+ ** Each entry thereby ends up in three of the six partitions.
+ */
+ POSTGIS_DEBUG(4, "[GIST] 'picksplit' calculating best split axis");
+ for ( i = FirstOffsetNumber; i <= max_offset; i = OffsetNumberNext(i) )
+ {
+ box_current = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
+
+ for ( d = 0; d < ndims_pageunion; d++ )
+ {
+ if( gidx_get_min(box_current,d)-gidx_get_min(box_pageunion,d) < gidx_get_max(box_pageunion,d)-gidx_get_max(box_current,d) )
+ {
+ geography_gist_picksplit_addlist(list[BELOW(d)], &(box_union[BELOW(d)]), box_current, &(pos[BELOW(d)]), i);
+ }
+ else
+ {
+ geography_gist_picksplit_addlist(list[ABOVE(d)], &(box_union[ABOVE(d)]), box_current, &(pos[ABOVE(d)]), i);
+ }
+
+ }
+
+ }
+
+ /*
+ ** "Bad disposition", too many entries fell into one octant of the space, so no matter which
+ ** plane we choose to split on, we're going to end up with a mostly full node. Where the
+ ** data is pretty homogeneous (lots of duplicates) entries that are equidistant from the
+ ** sides of the page union box can occasionally all end up in one place, leading
+ ** to this condition.
+ */
+ if ( geography_gist_picksplit_badratios(pos,ndims_pageunion) == TRUE )
+ {
+ /*
+ ** Instead we split on center points and see if we do better.
+ ** First calculate the average center point for each axis.
+ */
+ double *avgCenter = palloc(ndims_pageunion * sizeof(double));
+
+ for( d = 0; d < ndims_pageunion; d++ )
+ {
+ avgCenter[d] = 0.0;
+ }
+
+ POSTGIS_DEBUG(4, "[GIST] picksplit can't find good split axis, trying center point method");
+
+ for( i = FirstOffsetNumber; i <= max_offset; i = OffsetNumberNext(i) )
+ {
+ box_current = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
+ for( d = 0; d < ndims_pageunion; d++ )
+ {
+ avgCenter[d] += (gidx_get_max(box_current,d) + gidx_get_min(box_current,d)) / 2.0;
+ }
+ }
+ for( d = 0; d < ndims_pageunion; d++ )
+ {
+ avgCenter[d] /= max_offset;
+ pos[BELOW(d)] = pos[ABOVE(d)] = 0; /* Re-initialize our counters. */
+ POSTGIS_DEBUGF(4, "[GIST] picksplit average center point[%d] = %.12g", d, avgCenter[d]);
+ }
+
+ /* For each of our entries... */
+ for( i = FirstOffsetNumber; i <= max_offset; i = OffsetNumberNext(i) )
+ {
+ double center;
+ box_current = (GIDX*) DatumGetPointer(entryvec->vector[i].key);
+
+ for( d = 0; d < ndims_pageunion; d++ )
+ {
+ center = (gidx_get_min(box_current,d)+gidx_get_max(box_current,d))/2.0;
+ if( center < avgCenter[d] )
+ geography_gist_picksplit_addlist(list[BELOW(d)], &(box_union[BELOW(d)]), box_current, &(pos[BELOW(d)]), i);
+ else if( FPeq(center, avgCenter[d]) )
+ if( pos[BELOW(d)] > pos[ABOVE(d)] )
+ geography_gist_picksplit_addlist(list[ABOVE(d)], &(box_union[ABOVE(d)]), box_current, &(pos[ABOVE(d)]), i);
+ else
+ geography_gist_picksplit_addlist(list[BELOW(d)], &(box_union[BELOW(d)]), box_current, &(pos[BELOW(d)]), i);
+ else
+ geography_gist_picksplit_addlist(list[ABOVE(d)], &(box_union[ABOVE(d)]), box_current, &(pos[ABOVE(d)]), i);
+ }
+
+ }
+
+ /* Do we have a good disposition now? If not, screw it, just cut the node in half. */
+ if ( geography_gist_picksplit_badratios(pos,ndims_pageunion) == TRUE )
+ {
+ POSTGIS_DEBUG(4, "[GIST] picksplit still cannot find a good split! just cutting the node in half");
+ geography_gist_picksplit_fallback(entryvec, v);
+ PG_RETURN_POINTER(v);
+ }
+
+ }
+
+ /*
+ ** Now, what splitting plane gives us the most even ratio of
+ ** entries in our child pages? Since each split region has been apportioned entries
+ ** against the same number of total entries, the axis that has the smallest maximum
+ ** number of entries in its regions is the most evenly distributed.
+ ** TODO: what if the distributions are equal in two or more axes?
+ */
+ for( d = 0; d < ndims_pageunion; d++ )
+ {
+ int posd = Max(pos[ABOVE(d)],pos[BELOW(d)]);
+ if( posd > posmax )
+ direction = d;
+ }
+ if( direction == -1 || posmax == -1 )
+ {
+ /* ERROR OUT HERE */
+ }
+
+ POSTGIS_DEBUGF(4, "[GIST] 'picksplit' splitting on axis %d",direction);
+ geography_gist_picksplit_constructsplit(v, list[BELOW(direction)],
+ pos[BELOW(direction)],
+ &(box_union[BELOW(direction)]),
+ list[ABOVE(direction)],
+ pos[ABOVE(direction)],
+ &(box_union[ABOVE(direction)]) );
+
+ PG_RETURN_POINTER(v);
+
+}
+
+/*
+** The GIDX key must be defined as a PostgreSQL type, even though it is only
+** ever used internally. These no-op stubs are used to bind the type.
+*/
+PG_FUNCTION_INFO_V1(gidx_in);
+Datum gidx_in(PG_FUNCTION_ARGS)
+{
+// char *str = PG_GETARG_CSTRING(0);
+ PG_RETURN_POINTER(NULL);
+}
+
+PG_FUNCTION_INFO_V1(gidx_out);
+Datum gidx_out(PG_FUNCTION_ARGS)
+{
+ char *str = palloc(5);
+ SET_VARSIZE(str, 1);
+ str[4] = '\0';
+ PG_RETURN_POINTER(str);
+}
--- /dev/null
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "postgres.h"
+
+#include "../postgis_config.h"
+
+#include <math.h>
+#include <float.h>
+#include <string.h>
+#include <stdio.h>
+#include <errno.h>
+
+#include "utils/elog.h"
+#include "utils/array.h"
+#include "utils/builtins.h" /* for pg_atoi */
+#include "lib/stringinfo.h" /* For binary input */
+#include "catalog/pg_type.h" /* for CSTRINGOID */
+
+#include "libgeom.h" /* For standard geometry types. */
+#include "lwgeom_pg.h" /* For debugging macros. */
+#include "geography.h" /* For utility functions. */
+
+Datum geography_in(PG_FUNCTION_ARGS);
+Datum geography_out(PG_FUNCTION_ARGS);
+Datum geography_typmod_in(PG_FUNCTION_ARGS);
+Datum geography_typmod_out(PG_FUNCTION_ARGS);
+Datum geography_typmod_dims(PG_FUNCTION_ARGS);
+Datum geography_typmod_srid(PG_FUNCTION_ARGS);
+Datum geography_typmod_type(PG_FUNCTION_ARGS);
+
+Datum geography_enforce_typmod(PG_FUNCTION_ARGS);
+Datum geography_as_text(PG_FUNCTION_ARGS);
+Datum geography_from_text(PG_FUNCTION_ARGS);
+Datum geography_as_binary(PG_FUNCTION_ARGS);
+Datum geography_from_binary(PG_FUNCTION_ARGS);
+Datum geography_from_geometry(PG_FUNCTION_ARGS);
+
+/* Datum geography_gist_selectivity(PG_FUNCTION_ARGS); TBD */
+/* Datum geography_gist_join_selectivity(PG_FUNCTION_ARGS); TBD */
+/* Datum geography_send(PG_FUNCTION_ARGS); TBD */
+/* Datum geography_recv(PG_FUNCTION_ARGS); TBD */
+
+
+/**
+* Utility method to call the libgeom serialization and then set the
+* PgSQL varsize header appropriately with the serialized size.
+*/
+GSERIALIZED* geography_serialize(LWGEOM *lwgeom)
+{
+ static int is_geodetic = 1;
+ size_t ret_size = 0;
+ GSERIALIZED *g = NULL;
+
+ g = gserialized_from_lwgeom(lwgeom, is_geodetic, &ret_size);
+ if( ! g ) lwerror("Unable to serialize lwgeom.");
+ SET_VARSIZE(g, ret_size);
+ return g;
+}
+
+/**
+* The geography type only support POINT, LINESTRING, POLYGON, MULTI* variants
+* of same, and GEOMETRYCOLLECTION. If the input type is not one of those, shut
+* down the query.
+*/
+void geography_valid_type(uchar type)
+{
+ if( ! (
+ type == POINTTYPE ||
+ type == LINETYPE ||
+ type == POLYGONTYPE ||
+ type == MULTIPOINTTYPE ||
+ type == MULTILINETYPE ||
+ type == MULTIPOLYGONTYPE ||
+ type == COLLECTIONTYPE
+ ) )
+ {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Geography type does not support %s", lwgeom_typename(type) )));
+
+ }
+}
+
+/**
+* Check the consistency of the metadata we want to enforce in the typmod:
+* srid, type and dimensionality. If things are inconsistent, shut down the query.
+*/
+void geography_valid_typmod(LWGEOM *lwgeom, int32 typmod)
+{
+ int32 lwgeom_srid;
+ int32 lwgeom_type;
+ int32 lwgeom_z;
+ int32 lwgeom_m;
+ int32 typmod_srid = TYPMOD_GET_SRID(typmod);
+ int32 typmod_type = TYPMOD_GET_TYPE(typmod);
+ int32 typmod_z = TYPMOD_GET_Z(typmod);
+ int32 typmod_m = TYPMOD_GET_M(typmod);
+
+ assert(lwgeom);
+
+ lwgeom_type = TYPE_GETTYPE(lwgeom->type);
+ lwgeom_srid = lwgeom->SRID;
+ lwgeom_z = TYPE_HASZ(lwgeom->type);
+ lwgeom_m = TYPE_HASM(lwgeom->type);
+
+ POSTGIS_DEBUG(3, "Entered function");
+
+ /* No typmod (-1) => no preferences */
+ if (typmod < 0) return;
+
+ POSTGIS_DEBUGF(3, "Got lwgeom(type = %d, srid = %d, hasz = %d, hasm = %d)", lwgeom_type, lwgeom_srid, lwgeom_z, lwgeom_m);
+ POSTGIS_DEBUGF(3, "Got typmod(type = %d, srid = %d, hasz = %d, hasm = %d)", typmod_type, typmod_srid, typmod_z, typmod_m);
+
+ /* Typmod has a preference for SRID. */
+ if( typmod_srid > 0 && typmod_srid != lwgeom_srid)
+ {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Geography SRID (%d) does not match column SRID (%d)", lwgeom_srid, typmod_srid) ));
+ }
+
+ /* Typmod has a preference for geometry type. */
+ if( typmod_type > 0 &&
+ /* GEOMETRYCOLLECTION column can hold any kind of collection */
+ ((typmod_type == COLLECTIONTYPE && ! (lwgeom_type == COLLECTIONTYPE ||
+ lwgeom_type == MULTIPOLYGONTYPE ||
+ lwgeom_type == MULTIPOINTTYPE ||
+ lwgeom_type == MULTILINETYPE )) ||
+ /* Other types must be strictly equal. */
+ (typmod_type != lwgeom_type)) )
+ {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Geometry type (%s) does not match column type (%s)", lwgeom_typename(lwgeom_type), lwgeom_typename(typmod_type)) ));
+ }
+
+ /* Mismatched Z dimensionality. */
+ if( typmod_z && ! lwgeom_z )
+ {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Column has Z dimension but geometry does not" )));
+ }
+
+ /* Mismatched Z dimensionality (other way). */
+ if( lwgeom_z && ! typmod_z )
+ {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Geometry has Z dimension but column does not" )));
+ }
+
+ /* Mismatched M dimensionality. */
+ if( typmod_m && ! lwgeom_m )
+ {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Column has M dimension but geometry does not" )));
+ }
+
+ /* Mismatched M dimensionality (other way). */
+ if( lwgeom_m && ! typmod_m )
+ {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Geometry has M dimension but column does not" )));
+ }
+}
+
+/*
+** geography_in(cstring) returns *GSERIALIZED
+*/
+PG_FUNCTION_INFO_V1(geography_in);
+Datum geography_in(PG_FUNCTION_ARGS)
+{
+ char *geog_str = PG_GETARG_CSTRING(0);
+ /* Datum geog_oid = PG_GETARG_OID(1); Not needed. */
+ int32 geog_typmod = PG_GETARG_INT32(2);
+ LWGEOM_PARSER_RESULT lwg_parser_result;
+ LWGEOM *lwgeom = NULL;
+ GSERIALIZED *g_ser = NULL;
+ int result = 0;
+
+ /* Handles both HEXEWKB and EWKT */
+ result = serialized_lwgeom_from_ewkt(&lwg_parser_result, geog_str, PARSER_CHECK_ALL);
+ if (result)
+ PG_PARSER_ERROR(lwg_parser_result);
+
+ lwgeom = lwgeom_deserialize(lwg_parser_result.serialized_lwgeom);
+
+ geography_valid_type(TYPE_GETTYPE(lwgeom->type));
+
+ if( geog_typmod >= 0 )
+ {
+ geography_valid_typmod(lwgeom, geog_typmod);
+ POSTGIS_DEBUG(2,"typmod and geometry were consistent");
+ }
+ else
+ {
+ POSTGIS_DEBUG(2,"typmod was -1");
+ }
+
+ /*
+ ** Serialize our lwgeom and set the geodetic flag so subsequent
+ ** functions do the right thing.
+ */
+ g_ser = geography_serialize(lwgeom);
+ g_ser->flags = FLAGS_SET_GEODETIC(g_ser->flags, 1);
+
+ /*
+ ** Replace the unaligned lwgeom with a new aligned one based on GSERIALIZED.
+ */
+ lwgeom_release(lwgeom);
+ lwgeom = lwgeom_from_gserialized(g_ser);
+
+ /* Check if the geography has valid coordinate range. */
+ if( lwgeom_check_geodetic(lwgeom) == G_FALSE )
+ {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type" )));
+ }
+
+ PG_RETURN_POINTER(g_ser);
+}
+
+/*
+** geography_out(*GSERIALIZED) returns cstring
+*/
+PG_FUNCTION_INFO_V1(geography_out);
+Datum geography_out(PG_FUNCTION_ARGS)
+{
+ LWGEOM_UNPARSER_RESULT lwg_unparser_result;
+ LWGEOM *lwgeom = NULL;
+ GSERIALIZED *g = NULL;
+ int result = 0;
+
+ g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ lwgeom = lwgeom_from_gserialized(g);
+
+ result = serialized_lwgeom_to_hexwkb(&lwg_unparser_result, lwgeom_serialize(lwgeom), PARSER_CHECK_ALL, -1);
+ if (result)
+ PG_UNPARSER_ERROR(lwg_unparser_result);
+
+ PG_RETURN_CSTRING(lwg_unparser_result.wkoutput);
+}
+
+/*
+** geography_enforce_typmod(*GSERIALIZED) returns *GSERIALIZED
+*/
+PG_FUNCTION_INFO_V1(geography_enforce_typmod);
+Datum geography_enforce_typmod(PG_FUNCTION_ARGS)
+{
+ GSERIALIZED *arg = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ LWGEOM *lwgeom = NULL;
+ int32 typmod = PG_GETARG_INT32(1);
+ /* We don't need to have different behavior based on explicitness. */
+ /* bool isExplicit = PG_GETARG_BOOL(2); */
+
+ lwgeom = lwgeom_from_gserialized(arg);
+
+ /* Check if geometry typmod is consistent with the supplied one. */
+ geography_valid_typmod(lwgeom, typmod);
+
+ PG_RETURN_POINTER(geography_serialize(lwgeom));
+}
+
+/*
+** geography_typmod_in(cstring[]) returns int32
+**
+** Modified from ArrayGetIntegerTypmods in PostgreSQL 8.3
+*/
+PG_FUNCTION_INFO_V1(geography_typmod_in);
+Datum geography_typmod_in(PG_FUNCTION_ARGS)
+{
+
+ ArrayType *arr = (ArrayType *) DatumGetPointer(PG_GETARG_DATUM(0));
+ uint32 typmod = 0;
+ Datum *elem_values;
+ int n = 0;
+ int i = 0;
+
+ if (ARR_ELEMTYPE(arr) != CSTRINGOID)
+ ereport(ERROR,
+ (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
+ errmsg("typmod array must be type cstring[]")));
+
+ if (ARR_NDIM(arr) != 1)
+ ereport(ERROR,
+ (errcode(ERRCODE_ARRAY_SUBSCRIPT_ERROR),
+ errmsg("typmod array must be one-dimensional")));
+
+ if (ARR_HASNULL(arr))
+ ereport(ERROR,
+ (errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED),
+ errmsg("typmod array must not contain nulls")));
+
+ deconstruct_array(arr,
+ CSTRINGOID, -2, false, 'c', /* hardwire cstring representation details */
+ &elem_values, NULL, &n);
+
+ for (i = 0; i < n; i++)
+ {
+ if( i == 1 ) /* SRID */
+ {
+ int srid = pg_atoi(DatumGetCString(elem_values[i]), sizeof(int32), '\0');
+ if( srid > 0 )
+ {
+ POSTGIS_DEBUGF(1,"srid: %d", srid);
+ if( srid > SRID_MAXIMUM )
+ {
+ ereport(ERROR,
+ (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("SRID value may not exceed %d",
+ SRID_MAXIMUM)));
+ }
+ else
+ {
+ /* TODO: Check that the value provided is in fact a lonlat entry in spatial_ref_sys. */
+ /* For now, we only accept 4326. */
+ if( srid != 4326 )
+ {
+ ereport(ERROR,
+ (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Currently, only 4326 is accepted as an SRID for GEOGRAPHY")));
+ }
+ else
+ {
+ typmod = TYPMOD_SET_SRID(typmod, srid);
+ }
+ }
+ }
+ else
+ {
+ typmod = TYPMOD_SET_SRID(typmod, 0);
+ }
+ }
+ if( i == 0 ) /* TYPE */
+ {
+ char *s = DatumGetCString(elem_values[i]);
+ int type = 0;
+ int z = 0;
+ int m = 0;
+
+ if( geometry_type_from_string(s, &type, &z, &m) == G_FAILURE )
+ {
+ ereport(ERROR,
+ (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Invalid geometry type modifier: %s", s)));
+ }
+ else
+ {
+ typmod = TYPMOD_SET_TYPE(typmod, type);
+ if ( z )
+ typmod = TYPMOD_SET_Z(typmod);
+ if ( m )
+ typmod = TYPMOD_SET_M(typmod);
+ }
+ }
+ }
+
+ pfree(elem_values);
+
+ PG_RETURN_INT32(typmod);
+
+}
+
+
+/*
+** geography_typmod_out(int) returns cstring
+*/
+PG_FUNCTION_INFO_V1(geography_typmod_out);
+Datum geography_typmod_out(PG_FUNCTION_ARGS)
+{
+ char *s = (char*)palloc(64);
+ char *str = s;
+ uint32 typmod = PG_GETARG_INT32(0);
+ uint32 srid = TYPMOD_GET_SRID(typmod);
+ uint32 type = TYPMOD_GET_TYPE(typmod);
+ uint32 hasz = TYPMOD_GET_Z(typmod);
+ uint32 hasm = TYPMOD_GET_M(typmod);
+
+ POSTGIS_DEBUGF(3,"Got typmod(srid = %d, type = %d, hasz = %d, hasm = %d)", srid, type, hasz, hasm);
+
+ /* No SRID or type or dimensionality? Then no typmod at all. Return empty string. */
+ if( ! srid && ! type && ! hasz && ! hasz )
+ {
+ *str = '\0';
+ PG_RETURN_CSTRING(str);
+ }
+
+ /* Opening bracket. */
+ str += sprintf(str, "(");
+
+ /* Has type? */
+ if( type )
+ str += sprintf(str, "%s", lwgeom_typename(type));
+
+ /* Need dummy type to append Z/M to? */
+ if( !type & (hasz || hasz) )
+ str += sprintf(str, "Geometry");
+
+ /* Has Z? */
+ if( hasz )
+ str += sprintf(str, "%s", "Z");
+
+ /* Has M? */
+ if( hasm )
+ str += sprintf(str, "%s", "M");
+
+ /* Comma? */
+ if( srid && ( type || hasz || hasm ) )
+ str += sprintf(str, ",");
+
+ /* Has SRID? */
+ if( srid )
+ str += sprintf(str, "%d", srid);
+
+ /* Closing bracket. */
+ str += sprintf(str, ")");
+
+ PG_RETURN_CSTRING(s);
+
+}
+
+/*
+** geography_as_text(*GSERIALIZED) returns cstring
+*/
+PG_FUNCTION_INFO_V1(geography_as_text);
+Datum geography_as_text(PG_FUNCTION_ARGS)
+{
+ LWGEOM *lwgeom = NULL;
+ GSERIALIZED *g = NULL;
+ int result = 0;
+ char *semicolon_loc = NULL;
+ char *wkt = NULL;
+ uchar *lwgeom_serialized = NULL;
+ LWGEOM_UNPARSER_RESULT lwg_unparser_result;
+ size_t len = 0;
+
+ g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+
+ /* Convert to lwgeom so we can run the old functions */
+ lwgeom = lwgeom_from_gserialized(g);
+ lwgeom_serialized = lwgeom_serialize(lwgeom);
+
+ /* Generate WKT */
+ result = serialized_lwgeom_to_ewkt(&lwg_unparser_result, lwgeom_serialized, PARSER_CHECK_ALL);
+ if (result)
+ PG_UNPARSER_ERROR(lwg_unparser_result);
+
+ /* Strip SRID=NNNN;' off the front of the EWKT. */
+ semicolon_loc = strchr(lwg_unparser_result.wkoutput,';');
+ if (semicolon_loc == NULL)
+ semicolon_loc = lwg_unparser_result.wkoutput;
+ else
+ semicolon_loc = semicolon_loc + 1;
+
+ len = strlen(semicolon_loc) + VARHDRSZ;
+ wkt = palloc(len);
+ SET_VARSIZE(wkt, len);
+
+ memcpy(VARDATA(wkt), semicolon_loc, len - VARHDRSZ);
+
+ pfree(lwg_unparser_result.wkoutput);
+ pfree(lwgeom_serialized);
+ lwgeom_release(lwgeom);
+
+ PG_RETURN_POINTER(wkt);
+}
+
+/*
+** geography_from_text(*char) returns *GSERIALIZED
+**
+** Convert text (varlena) to cstring and then call geography_in().
+*/
+PG_FUNCTION_INFO_V1(geography_from_text);
+Datum geography_from_text(PG_FUNCTION_ARGS)
+{
+ text *wkt_text = PG_GETARG_TEXT_P(0);
+ size_t size = VARSIZE(wkt_text) - VARHDRSZ;
+ /* Extract the cstring from the varlena */
+ char *wkt = palloc(size + 1);
+ memcpy(wkt, VARDATA(wkt_text), size);
+ /* Null terminate it */
+ wkt[size] = '\0';
+
+ if ( size < 10 )
+ {
+ lwerror("Invalid OGC WKT (too short)");
+ PG_RETURN_NULL();
+ }
+
+ /* Pass the cstring to the input parser, and magic occurs! */
+ PG_RETURN_DATUM(DirectFunctionCall3(geography_in, PointerGetDatum(wkt), Int32GetDatum(0), Int32GetDatum(-1)));
+}
+
+/*
+** geography_as_binary(*GSERIALIZED) returns *char
+*/
+PG_FUNCTION_INFO_V1(geography_as_binary);
+Datum geography_as_binary(PG_FUNCTION_ARGS)
+{
+ LWGEOM_UNPARSER_RESULT lwg_unparser_result;
+ LWGEOM *lwgeom = NULL;
+ uchar *lwgeom_serialized = NULL;
+ size_t lwgeom_serialized_size = 0;
+ uchar *lwgeom_serialized_2d = NULL;
+ int result = 0;
+ char *wkb = NULL;
+ size_t wkb_size = 0;
+ GSERIALIZED *g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+
+ /* Drop SRID so that WKB does not contain SRID. */
+ gserialized_set_srid(g, 0);
+
+ /* Convert to lwgeom so we can run the old functions */
+ lwgeom = lwgeom_from_gserialized(g);
+ lwgeom_serialized_size = lwgeom_serialize_size(lwgeom);
+ lwgeom_serialized = lwgeom_serialize(lwgeom);
+
+ /* Force to 2D */
+ lwgeom_serialized_2d = lwalloc(lwgeom_serialized_size);
+ lwgeom_force2d_recursive(lwgeom_serialized, lwgeom_serialized_2d, &lwgeom_serialized_size);
+
+ /* Create WKB */
+ result = serialized_lwgeom_to_ewkb(&lwg_unparser_result, lwgeom_serialized_2d, PARSER_CHECK_ALL, NDR);
+ if (result)
+ PG_UNPARSER_ERROR(lwg_unparser_result);
+
+ /* Copy to varlena pointer */
+ wkb_size = lwg_unparser_result.size + VARHDRSZ;
+ wkb = palloc(wkb_size);
+ SET_VARSIZE(wkb, wkb_size);
+ memcpy(VARDATA(wkb), lwg_unparser_result.wkoutput, lwg_unparser_result.size);
+
+ /* Clean up */
+ pfree(lwg_unparser_result.wkoutput);
+ lwgeom_release(lwgeom);
+ lwfree(lwgeom_serialized);
+ lwfree(lwgeom_serialized_2d);
+
+ PG_RETURN_POINTER(wkb);
+}
+
+/*
+** geography_from_binary(*char) returns *GSERIALIZED
+*/
+PG_FUNCTION_INFO_V1(geography_from_binary);
+Datum geography_from_binary(PG_FUNCTION_ARGS)
+{
+ char *wkb_cstring = NULL;
+ char *wkb_hex = NULL;
+ char *wkb_bytea = (char*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ char *hexarg = palloc(4 + VARHDRSZ);
+ size_t wkb_hex_size;
+
+ /* Create our second argument to binary_encode */
+ SET_VARSIZE(hexarg, 4 + VARHDRSZ);
+ memcpy(VARDATA(hexarg), "hex", 4);
+
+ /* Convert the bytea into a hex representation cstring */
+ wkb_hex = DatumGetPointer(DirectFunctionCall2(binary_encode, PointerGetDatum(wkb_bytea), PointerGetDatum(hexarg)));
+ wkb_hex_size = VARSIZE(wkb_hex) - VARHDRSZ;
+ wkb_cstring = palloc(wkb_hex_size + 1);
+ memcpy(wkb_cstring, VARDATA(wkb_hex), wkb_hex_size);
+ wkb_cstring[wkb_hex_size] = '\0'; /* Null terminate the cstring */
+ pfree(hexarg);
+
+ /* Pass the cstring to the input parser, and magic occurs! */
+ PG_RETURN_DATUM(DirectFunctionCall3(geography_in, PointerGetDatum(wkb_cstring), Int32GetDatum(0), Int32GetDatum(-1)));
+}
+
+PG_FUNCTION_INFO_V1(geography_typmod_dims);
+Datum geography_typmod_dims(PG_FUNCTION_ARGS)
+{
+ int32 typmod = PG_GETARG_INT32(0);
+ int32 dims = 2;
+ if( typmod < 0 )
+ PG_RETURN_INT32(dims);
+ if( TYPMOD_GET_Z(typmod) )
+ dims++;
+ if( TYPMOD_GET_M(typmod) )
+ dims++;
+ PG_RETURN_INT32(dims);
+}
+
+PG_FUNCTION_INFO_V1(geography_typmod_srid);
+Datum geography_typmod_srid(PG_FUNCTION_ARGS)
+{
+ int32 typmod = PG_GETARG_INT32(0);
+ if( typmod < 0 )
+ PG_RETURN_INT32(0);
+ PG_RETURN_INT32(TYPMOD_GET_SRID(typmod));
+}
+
+PG_FUNCTION_INFO_V1(geography_typmod_type);
+Datum geography_typmod_type(PG_FUNCTION_ARGS)
+{
+ int32 typmod = PG_GETARG_INT32(0);
+ int32 type = TYPMOD_GET_TYPE(typmod);
+ char *s = (char*)palloc(64);
+ char *str = s;
+ int slen = 0;
+
+ /* Has type? */
+ if( typmod < 0 || type == 0 )
+ str += sprintf(str, "Geometry");
+ else
+ str += sprintf(str, "%s", lwgeom_typename(type));
+
+ /* Has Z? */
+ if( typmod >= 0 && TYPMOD_GET_Z(typmod) )
+ str += sprintf(str, "%s", "Z");
+
+ /* Has M? */
+ if( typmod >= 0 && TYPMOD_GET_M(typmod) )
+ str += sprintf(str, "%s", "M");
+
+ slen = strlen(s) + 1;
+ str = palloc(slen + VARHDRSZ);
+ SET_VARSIZE(str, slen + VARHDRSZ);
+ memcpy(VARDATA(str), s, slen);
+ pfree(s);
+ PG_RETURN_POINTER(str);
+}
+
+PG_FUNCTION_INFO_V1(geography_from_geometry);
+Datum geography_from_geometry(PG_FUNCTION_ARGS)
+{
+ LWGEOM *lwgeom = NULL;
+ GSERIALIZED *g_ser = NULL;
+ uchar *lwgeom_serialized = (uchar*)VARDATA(PG_DETOAST_DATUM(PG_GETARG_DATUM(0)));
+
+ geography_valid_type(TYPE_GETTYPE(lwgeom_serialized[0]));
+
+ lwgeom = lwgeom_deserialize(lwgeom_serialized);
+
+ /*
+ ** Serialize our lwgeom and set the geodetic flag so subsequent
+ ** functions do the right thing.
+ */
+ g_ser = geography_serialize(lwgeom);
+ g_ser->flags = FLAGS_SET_GEODETIC(g_ser->flags, 1);
+
+ /*
+ ** Replace the unaligned lwgeom with a new aligned one based on GSERIALIZED.
+ */
+ lwgeom_release(lwgeom);
+ lwgeom = lwgeom_from_gserialized(g_ser);
+
+ /* Check if the geography has valid coordinate range. */
+ if( lwgeom_check_geodetic(lwgeom) == G_FALSE )
+ {
+ ereport(ERROR, (
+ errcode(ERRCODE_INVALID_PARAMETER_VALUE),
+ errmsg("Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type" )));
+ }
+
+ PG_RETURN_POINTER(g_ser);
+
+}
+
bool *recheck = (bool *) PG_GETARG_POINTER(4);
/* We set recheck to false to avoid repeatedly pulling every "possibly matched" geometry
- out during index scans. For cases when the geometries are large, doing this
+ out during index scans. For cases when the geometries are large, rechecking
can make things twice as slow. */
*recheck = false;
#endif
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--- $Id: sqlmm.sql.in 2406 2006-11-02 13:56:52Z kneufeld $
+-- $Id$
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.refractions.net
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--- $Id: long_xact.sql.in 2803 2008-06-05 12:09:48Z mcayland $
+-- $Id$
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.refractions.net
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--- $Id: uninstall_postgis.sql.in.c 3993 2009-04-12 10:56:33Z mcayland $
+-- $Id$
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.refractions.net
-- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
--
--- $Id: sqlmm.sql.in 2406 2006-11-02 13:56:52Z kneufeld $
+-- $Id$
--
-- PostGIS - Spatial Types for PostgreSQL
-- http://postgis.refractions.net
projects / postgis / commitdiff