From: Paul Ramsey Date: Mon, 19 Oct 2009 05:05:47 +0000 (+0000) Subject: Add _ST_BestSRID(Geography) utility function to support ST_Buffer(geography, radius... X-Git-Tag: 1.5.0b1~356 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=e40e70733581df6cd322d7789decc8929d61c7b9;p=postgis Add _ST_BestSRID(Geography) utility function to support ST_Buffer(geography, radius) hack that casts back and forth to geometry. git-svn-id: http://svn.osgeo.org/postgis/trunk@4664 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c index 70df83119..28bfe3c2c 100644 --- a/postgis/geography.sql.in.c +++ b/postgis/geography.sql.in.c @@ -498,6 +498,26 @@ CREATE OR REPLACE FUNCTION ST_CoveredBy(geography, geography) AS 'SELECT $1 && $2 AND _ST_Covers($2, $1)' LANGUAGE 'SQL' IMMUTABLE STRICT; +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION _ST_BestSRID(geography, geography) + RETURNS integer + AS 'MODULE_PATHNAME','geography_bestsrid' + LANGUAGE 'C' IMMUTABLE STRICT; + +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION _ST_BestSRID(geography) + RETURNS integer + AS 'SELECT _ST_BestSRID($1,$1)' + LANGUAGE 'SQL' IMMUTABLE STRICT; + +-- Availability: 1.5.0 +CREATE OR REPLACE FUNCTION ST_Buffer(geography, float8) + RETURNS geography + AS 'SELECT geography(ST_Transform(ST_Buffer(ST_Transform(geometry($1), _ST_BestSRID($1)), $2), 4326))' + LANGUAGE 'SQL' IMMUTABLE STRICT; + + + -- ---------- ---------- ---------- ---------- ---------- ---------- ---------- COMMIT; diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index 2a2d2e806..600fe3fac 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -29,6 +29,7 @@ Datum geography_length_sphere(PG_FUNCTION_ARGS); Datum geography_expand(PG_FUNCTION_ARGS); Datum geography_point_outside(PG_FUNCTION_ARGS); Datum geography_covers(PG_FUNCTION_ARGS); +Datum geography_bestsrid(PG_FUNCTION_ARGS); /* ** geography_distance_sphere(GSERIALIZED *g1, GSERIALIZED *g2, double tolerance) @@ -250,6 +251,13 @@ Datum geography_point_outside(PG_FUNCTION_ARGS) } +/* +** geography_covers(GSERIALIZED *g, GSERIALIZED *g) returns boolean +** Only works for (multi)points and (multi)polygons currently. +** Attempts a simple point-in-polygon test on the polygon and point. +** Current algorithm does not distinguish between points on edge +** and points within. +*/ PG_FUNCTION_INFO_V1(geography_covers); Datum geography_covers(PG_FUNCTION_ARGS) { @@ -299,3 +307,118 @@ Datum geography_covers(PG_FUNCTION_ARGS) PG_RETURN_BOOL(result); } + +/* +** geography_bestsrid(GSERIALIZED *g, GSERIALIZED *g) returns int +** Utility function. Returns negative SRID numbers that match to the +** numbers handled in code by the transform(lwgeom, srid) function. +** UTM, polar stereographic and mercator as fallback. To be used +** in wrapping existing geometry functions in SQL to provide access +** to them in the geography module. +*/ +PG_FUNCTION_INFO_V1(geography_bestsrid); +Datum geography_bestsrid(PG_FUNCTION_ARGS) +{ + LWGEOM *lwgeom1 = NULL; + LWGEOM *lwgeom2 = NULL; + GBOX gbox1; + GBOX gbox2; + GSERIALIZED *g1 = NULL; + GSERIALIZED *g2 = NULL; + int type1, type2; + int empty1 = LW_FALSE; + int empty2 = LW_FALSE; + int twoargs = LW_FALSE; + + Datum d1 = PG_GETARG_DATUM(0); + Datum d2 = PG_GETARG_DATUM(1); + + /* Get our geometry objects loaded into memory. */ + g1 = (GSERIALIZED*)PG_DETOAST_DATUM(d1); + /* Synchronize our box types */ + gbox1.flags = g1->flags; + /* Read our types */ + type1 = gserialized_get_type(g1); + /* Construct our working geometries */ + lwgeom1 = lwgeom_from_gserialized(g1); + /* Calculate if the geometry is empty. */ + empty1 = lwgeom_is_empty(lwgeom1); + /* Calculate a naive cartesian bounds for the objects */ + if( ! empty1 && lwgeom_calculate_gbox(lwgeom1, &gbox1) == G_FAILURE ) + elog(ERROR, "Error in geography_bestsrid calling lwgeom_calculate_gbox(lwgeom1, &gbox1)"); + + /* If we have a unique second argument, fill in all the necessarily variables. */ + if( d1 != d2 ) + { + g2 = (GSERIALIZED*)PG_DETOAST_DATUM(d2); + type2 = gserialized_get_type(g2); + gbox2.flags = g2->flags; + lwgeom2 = lwgeom_from_gserialized(g2); + empty2 = lwgeom_is_empty(lwgeom2); + if( ! empty2 && lwgeom_calculate_gbox(lwgeom2, &gbox2) == G_FAILURE ) + elog(ERROR, "Error in geography_bestsrid calling lwgeom_calculate_gbox(lwgeom2, &gbox2)"); + } + /* + ** If no unique second argument, copying the box for the first + ** argument will give us the right answer for all subsequent tests. + */ + else + { + gbox2 = gbox1; + } + + /* Both empty? We don't have an answer. */ + if ( empty1 && empty2 ) + PG_RETURN_NULL(); + + /* One empty? We can use the other argument values as infill. */ + if ( empty1 ) + gbox1 = gbox2; + + if ( empty2 ) + gbox2 = gbox1; + + + /* Are these data arctic? Polar stereographic north. */ + if( gbox1.ymin > 60.0 && gbox2.ymin > 60.0 ) + { + PG_RETURN_INT32(-3995); + } + + /* Are these data antarctic? Polar stereographic south. */ + if( gbox1.ymin < -60.0 && gbox2.ymin < -60.0 ) + { + PG_RETURN_INT32(-3995); + } + + /* + ** Can we fit these data into one UTM zone? We will assume we can push things as + ** far as a half zone past a zone boundary. Note we have no handling for the + ** date line in here. + */ + if( fabs(FP_MAX(gbox1.xmax, gbox2.xmax) - FP_MIN(gbox1.xmin, gbox2.xmin)) < 6.0 ) + { + /* Cheap hack to pick a zone. Average of the box center points. */ + double dzone = (gbox1.xmin + gbox1.xmax + gbox2.xmin + gbox2.xmax) / 4.0; + int zone = floor(1.0 + ((dzone + 180.0) / 6.0)); + + /* Are these data below the equator? UTM South. */ + if( gbox1.ymax < 0.0 && gbox2.ymax < 0.0 ) + { + PG_RETURN_INT32( -32700 - zone ); + } + /* Are these data above the equator? UTM North. */ + else + { + PG_RETURN_INT32( -32600 - zone ); + } + } + + /* + ** Running out of options... fall-back to Mercator and hope for the best. + */ + PG_RETURN_INT32(-3395); + +} + + diff --git a/postgis/lwgeom_transform.c b/postgis/lwgeom_transform.c index 8fcc10373..579d16f97 100644 --- a/postgis/lwgeom_transform.c +++ b/postgis/lwgeom_transform.c @@ -428,7 +428,12 @@ static char* GetProj4String(int srid) { strncpy(proj_str, "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", maxproj4len ); } - + /* World Mercator */ + if( id == 3395 ) + { + strncpy(proj_str, "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", maxproj4len ); + } + POSTGIS_DEBUGF(3, "returning on SRID=%d: %s", srid, proj_str); return proj_str; }