From 256159f3b933dc160cec287b8e5ffcc464695801 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Fri, 6 Jul 2012 23:50:33 +0000 Subject: [PATCH] Handle best SRID finding for shapes the cross the poles and dateline. Also add some new larger zones in gnomic for larger shapes. (#1610) git-svn-id: http://svn.osgeo.org/postgis/trunk@10038 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/cunit/cu_geodetic.c | 46 +++++++++++- liblwgeom/g_box.c | 27 ++++--- liblwgeom/liblwgeom_internal.h | 26 ++++--- liblwgeom/lwgeodetic.c | 120 ++++++++++++++++++++++++++++++++ postgis/geography_measurement.c | 71 ++++++++++--------- regress/bestsrid_expected | 4 +- 6 files changed, 231 insertions(+), 63 deletions(-) diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index 6919cacc6..4cd708bef 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -980,7 +980,7 @@ static void test_spheroid_area(void) SPHEROID s; /* Init to WGS84 */ - spheroid_init(&s, 6378137.0, 6356752.314245179498); + spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS); gbox.flags = gflags(0, 0, 1); @@ -1035,6 +1035,49 @@ static void test_spheroid_area(void) lwgeom_free(lwg); } +static void test_gbox_utils(void) +{ + LWGEOM *lwg; + GBOX gbox; + double a1, a2; + SPHEROID s; + POINT2D pt; + + /* Init to WGS84 */ + spheroid_init(&s, WGS84_MAJOR_AXIS, WGS84_MINOR_AXIS); + + gbox.flags = gflags(0, 0, 1); + + /* One-degree square by equator */ + lwg = lwgeom_from_wkt("POLYGON((1 20,1 21,2 21,2 20,1 20))", LW_PARSER_CHECK_NONE); + lwgeom_calculate_gbox_geodetic(lwg, &gbox); + a1 = gbox_angular_width(&gbox); + a2 = gbox_angular_height(&gbox); + CU_ASSERT_DOUBLE_EQUAL(a1, 0.0177951, 0.0000001); + CU_ASSERT_DOUBLE_EQUAL(a2, 0.017764, 0.0000001); + lwgeom_free(lwg); + + /* One-degree square *across* dateline */ + lwg = lwgeom_from_wkt("POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", LW_PARSER_CHECK_NONE); + lwgeom_calculate_gbox_geodetic(lwg, &gbox); + a1 = gbox_angular_width(&gbox); + a2 = gbox_angular_height(&gbox); + //printf("a1=%g a2=%g\n", a1, a2); + CU_ASSERT_DOUBLE_EQUAL(a1, 0.0174613, 0.0000001); + CU_ASSERT_DOUBLE_EQUAL(a2, 0.0174553, 0.0000001); + lwgeom_free(lwg); + + /* One-degree square *across* dateline */ + lwg = lwgeom_from_wkt("POLYGON((178.5 2,178.5 1,-179.5 1,-179.5 2,178.5 2))", LW_PARSER_CHECK_NONE); + lwgeom_calculate_gbox_geodetic(lwg, &gbox); + gbox_centroid(&gbox, &pt); + //printf("POINT(%g %g)\n", pt.x, pt.y); + CU_ASSERT_DOUBLE_EQUAL(pt.x, 179.5, 0.0001); + CU_ASSERT_DOUBLE_EQUAL(pt.y, 1.50024, 0.0001); + lwgeom_free(lwg); + +} + /* ** Used by test harness to register the tests in this file. @@ -1057,6 +1100,7 @@ CU_TestInfo geodetic_tests[] = PG_TEST(test_spheroid_area), PG_TEST(test_lwpoly_covers_point2d), PG_TEST(test_ptarray_point_in_ring), + PG_TEST(test_gbox_utils), CU_TEST_INFO_NULL }; CU_SuiteInfo geodetic_suite = {"Geodetic Suite", NULL, NULL, geodetic_tests}; diff --git a/liblwgeom/g_box.c b/liblwgeom/g_box.c index eff95f7a4..10aec67ef 100644 --- a/liblwgeom/g_box.c +++ b/liblwgeom/g_box.c @@ -92,33 +92,32 @@ void gbox_expand(GBOX *g, double d) int gbox_union(const GBOX *g1, const GBOX *g2, GBOX *gout) { - if ( (g1 == NULL) && (g2 == NULL) ) + if ( ( ! g1 ) && ( ! g2 ) ) return LW_FALSE; - if (g1 == NULL) + if ( ! g1 ) { memcpy(gout, g2, sizeof(GBOX)); return LW_TRUE; } - if (g2 == NULL) + if ( ! g2 ) { memcpy(gout, g1, sizeof(GBOX)); return LW_TRUE; } + + gout->flags = g1->flags; - if (g1->xmin < g2->xmin) gout->xmin = g1->xmin; - else gout->xmin = g2->xmin; - - if (g1->ymin < g2->ymin) gout->ymin = g1->ymin; - else gout->ymin = g2->ymin; - - if (g1->xmax > g2->xmax) gout->xmax = g1->xmax; - else gout->xmax = g2->xmax; + gout->xmin = FP_MIN(g1->xmin, g2->xmin); + gout->xmax = FP_MAX(g1->xmax, g2->xmax); - if (g1->ymax > g2->ymax) gout->ymax = g1->ymax; - else gout->ymax = g2->ymax; + gout->ymin = FP_MIN(g1->ymin, g2->ymin); + gout->ymax = FP_MAX(g1->ymax, g2->ymax); + + gout->zmin = FP_MIN(g1->zmin, g2->zmin); + gout->zmax = FP_MAX(g1->zmax, g2->zmax); - return LW_TRUE; + return LW_TRUE; } int gbox_same(const GBOX *g1, const GBOX *g2) diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h index 85c7f5bec..e56cdb9e2 100644 --- a/liblwgeom/liblwgeom_internal.h +++ b/liblwgeom/liblwgeom_internal.h @@ -353,16 +353,16 @@ int lwpsurface_is_closed(const LWPSURFACE *psurface); int lwtin_is_closed(const LWTIN *tin); -/* - * Split a line by a point and push components to the provided multiline. - * If the point doesn't split the line, push nothing to the container. - * Returns 0 if the point is off the line. - * Returns 1 if the point is on the line boundary (endpoints). - * Return 2 if the point is on the interior of the line (only case in which - * a split happens). - * - * NOTE: the components pushed to the output vector have their SRID stripped - */ +/** +* Split a line by a point and push components to the provided multiline. +* If the point doesn't split the line, push nothing to the container. +* Returns 0 if the point is off the line. +* Returns 1 if the point is on the line boundary (endpoints). +* Return 2 if the point is on the interior of the line (only case in which +* a split happens). +* +* NOTE: the components pushed to the output vector have their SRID stripped +*/ int lwline_split_by_point_to(const LWLINE* ln, const LWPOINT* pt, LWMLINE* to); /** Ensure the collection can hold at least up to ngeoms geometries */ @@ -371,4 +371,10 @@ void lwcollection_reserve(LWCOLLECTION *col, int ngeoms); /** Check if subtype is allowed in collectiontype */ extern int lwcollection_allows_subtype(int collectiontype, int subtype); +/** GBOX utility functions to figure out coverage/location on the globe */ +double gbox_angular_height(const GBOX* gbox); +double gbox_angular_width(const GBOX* gbox); +int gbox_centroid(const GBOX* gbox, POINT2D* out); + + #endif /* _LIBLWGEOM_INTERNAL_H */ diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index ff0e2b711..0faa8ca2d 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -156,6 +156,126 @@ void geographic_point_init(double lon, double lat, GEOGRAPHIC_POINT *g) g->lon = longitude_radians_normalize(deg2rad(lon)); } +/** Returns the angular height (latitudinal span) of the box in radians */ +double +gbox_angular_height(const GBOX* gbox) +{ + double d[6]; + int i; + double zmin = MAXFLOAT; + double zmax = -1 * MAXFLOAT; + POINT3D pt; + + /* Take a copy of the box corners so we can treat them as a list */ + /* Elements are xmin, xmax, ymin, ymax, zmin, zmax */ + memcpy(d, &(gbox->xmin), 6*sizeof(double)); + + /* Generate all 8 corner vectors of the box */ + for ( i = 0; i < 8; i++ ) + { + pt.x = d[i / 4]; + pt.y = d[2 + (i % 4) / 2]; + pt.z = d[4 + (i % 2)]; + normalize(&pt); + if ( pt.z < zmin ) zmin = pt.z; + if ( pt.z > zmax ) zmax = pt.z; + } + return asin(zmax) - asin(zmin); +} + +/** Returns the angular width (longitudinal span) of the box in radians */ +double +gbox_angular_width(const GBOX* gbox) +{ + double d[6]; + int i, j; + POINT3D pt[3]; + double maxangle; + double magnitude; + + /* Take a copy of the box corners so we can treat them as a list */ + /* Elements are xmin, xmax, ymin, ymax, zmin, zmax */ + memcpy(d, &(gbox->xmin), 6*sizeof(double)); + + /* Start with the bottom corner */ + pt[0].x = gbox->xmin; + pt[0].y = gbox->ymin; + magnitude = sqrt(pt[0].x*pt[0].x + pt[0].y*pt[0].y); + pt[0].x /= magnitude; + pt[0].y /= magnitude; + + /* Generate all 8 corner vectors of the box */ + /* Find the vector furthest from our seed vector */ + for ( j = 0; j < 2; j++ ) + { + maxangle = -1 * MAXFLOAT; + for ( i = 0; i < 4; i++ ) + { + double angle, dotprod; + POINT3D pt_n; + + pt_n.x = d[i / 2]; + pt_n.y = d[2 + (i % 2)]; + magnitude = sqrt(pt_n.x*pt_n.x + pt_n.y*pt_n.y); + pt_n.x /= magnitude; + pt_n.y /= magnitude; + + dotprod = pt_n.x*pt[j].x + pt_n.y*pt[j].y; + angle = acos(dotprod > 1.0 ? 1.0 : dotprod); + if ( angle > maxangle ) + { + pt[j+1] = pt_n; + maxangle = angle; + } + } + } + + /* Return the distance between the two furthest vectors */ + return maxangle; +} + +/** Computes the average(ish) center of the box and returns success. */ +int +gbox_centroid(const GBOX* gbox, POINT2D* out) +{ + double d[6]; + GEOGRAPHIC_POINT g; + POINT3D pt; + int i; + + /* Take a copy of the box corners so we can treat them as a list */ + /* Elements are xmin, xmax, ymin, ymax, zmin, zmax */ + memcpy(d, &(gbox->xmin), 6*sizeof(double)); + + /* Zero out our return vector */ + pt.x = pt.y = pt.z = 0.0; + + for ( i = 0; i < 8; i++ ) + { + POINT3D pt_n; + + pt_n.x = d[i / 4]; + pt_n.y = d[2 + ((i % 4) / 2)]; + pt_n.z = d[4 + (i % 2)]; + normalize(&pt_n); + + pt.x += pt_n.x; + pt.y += pt_n.y; + pt.z += pt_n.z; + } + + pt.x /= 8.0; + pt.y /= 8.0; + pt.z /= 8.0; + normalize(&pt); + + cart2geog(&pt, &g); + out->x = longitude_degrees_normalize(rad2deg(g.lon)); + out->y = latitude_degrees_normalize(rad2deg(g.lat)); + + return LW_SUCCESS; +} + /** * Check to see if this geocentric gbox is wrapped around a pole. * Only makes sense if this gbox originated from a polygon, as it's assuming diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c index fec3e49e4..9cb512572 100644 --- a/postgis/geography_measurement.c +++ b/postgis/geography_measurement.c @@ -687,16 +687,14 @@ Datum geography_covers(PG_FUNCTION_ARGS) PG_FUNCTION_INFO_V1(geography_bestsrid); Datum geography_bestsrid(PG_FUNCTION_ARGS) { - LWGEOM *lwgeom1 = NULL; - LWGEOM *lwgeom2 = NULL; - GBOX gbox1; - GBOX gbox2; + GBOX gbox, gbox1, gbox2; GSERIALIZED *g1 = NULL; GSERIALIZED *g2 = NULL; int type1, type2; int empty1 = LW_FALSE; int empty2 = LW_FALSE; - double ycenter, xcenter, xwidth, ywidth; + double xwidth, ywidth; + POINT2D center; Datum d1 = PG_GETARG_DATUM(0); Datum d2 = PG_GETARG_DATUM(1); @@ -707,26 +705,23 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS) 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_cartesian(lwgeom1, &gbox1) == LW_FAILURE ) - elog(ERROR, "Error in geography_bestsrid calling lwgeom_calculate_gbox(lwgeom1, &gbox1)"); - + empty1 = gserialized_is_empty(g1); + /* Calculate a geocentric bounds for the objects */ + if ( ! empty1 && gserialized_get_gbox_p(g1, &gbox1) == LW_FAILURE ) + elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g1, &gbox1)"); + POSTGIS_DEBUGF(4, "calculated gbox = %s", gbox_to_string(&gbox1)); - /* If we have a unique second argument, fill in all the necessarily variables. */ + /* If we have a unique second argument, fill in all the necessary 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_cartesian(lwgeom2, &gbox2) == LW_FAILURE ) - elog(ERROR, "Error in geography_bestsrid calling lwgeom_calculate_gbox(lwgeom2, &gbox2)"); + empty2 = gserialized_is_empty(g2); + if ( ! empty1 && gserialized_get_gbox_p(g2, &gbox2) == LW_FAILURE ) + elog(ERROR, "Error in geography_bestsrid calling gserialized_get_gbox_p(g2, &gbox2)"); } /* ** If no unique second argument, copying the box for the first @@ -734,34 +729,39 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS) */ else { - gbox2 = gbox1; + gbox = 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. */ + /* One empty? We can use the other argument values as infill. Otherwise merge the boxen */ if ( empty1 ) - gbox1 = gbox2; + gbox = gbox2; + else if ( empty2 ) + gbox = gbox1; + else + gbox_union(&gbox1, &gbox2, &gbox); - if ( empty2 ) - gbox2 = gbox1; + gbox_centroid(&gbox, ¢er); - ycenter = (gbox1.ymin + gbox1.ymax + gbox2.ymin + gbox2.ymax) / 4.0; - xcenter = (gbox1.xmin + gbox1.xmax + gbox2.xmin + gbox2.xmax) / 4.0; + /* Width and height in degrees */ + xwidth = 180.0 * gbox_angular_width(&gbox) / M_PI; + ywidth = 180.0 * gbox_angular_height(&gbox) / M_PI; - xwidth = fabs(FP_MAX(gbox1.xmax, gbox2.xmax) - FP_MIN(gbox1.xmin, gbox2.xmin)); - ywidth = fabs(FP_MAX(gbox1.ymax, gbox2.ymax) - FP_MIN(gbox1.ymin, gbox2.ymin)); + POSTGIS_DEBUGF(2, "xwidth %g", xwidth); + POSTGIS_DEBUGF(2, "ywidth %g", ywidth); + POSTGIS_DEBUGF(2, "center POINT(%g %g)", center.x, center.y); /* Are these data arctic? Lambert Azimuthal Equal Area North. */ - if ( gbox1.ymin > 65.0 && gbox2.ymin > 65.0 ) + if ( center.y > 70.0 && ywidth < 45.0 ) { PG_RETURN_INT32(SRID_NORTH_LAMBERT); } /* Are these data antarctic? Lambert Azimuthal Equal Area South. */ - if ( gbox1.ymin < -65.0 && gbox2.ymin < -65.0 ) + if ( center.y < -70.0 && ywidth < 45.0 ) { PG_RETURN_INT32(SRID_SOUTH_LAMBERT); } @@ -774,13 +774,12 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS) */ if ( xwidth < 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((dzone + 180.0) / 6.0); + int zone = floor((center.x + 180.0) / 6.0); + if ( zone > 59 ) zone = 59; /* Are these data below the equator? UTM South. */ - if ( gbox1.ymax < 0.0 && gbox2.ymax < 0.0 ) + if ( center.y < 0.0 ) { PG_RETURN_INT32( SRID_SOUTH_UTM_START + zone ); } @@ -800,22 +799,22 @@ Datum geography_bestsrid(PG_FUNCTION_ARGS) if ( ywidth < 25.0 ) { int xzone = -1; - int yzone = 3 + floor(ycenter / 30.0); /* (range of 0-5) */ + int yzone = 3 + floor(center.y / 30.0); /* (range of 0-5) */ /* Equatorial band, 12 zones, 30 degrees wide */ if ( (yzone == 2 || yzone == 3) && xwidth < 30.0 ) { - xzone = 6 + floor(xcenter / 30.0); + xzone = 6 + floor(center.x / 30.0); } /* Temperate band, 8 zones, 45 degrees wide */ else if ( (yzone == 1 || yzone == 4) && xwidth < 45.0 ) { - xzone = 4 + floor(xcenter / 45.0); + xzone = 4 + floor(center.x / 45.0); } /* Arctic band, 4 zones, 90 degrees wide */ else if ( (yzone == 0 || yzone == 5) && xwidth < 90.0 ) { - xzone = 2 + floor(xcenter / 90.0); + xzone = 2 + floor(center.x / 90.0); } /* Did we fit into an appropriate xzone? */ diff --git a/regress/bestsrid_expected b/regress/bestsrid_expected index 256038f9f..f939e24af 100644 --- a/regress/bestsrid_expected +++ b/regress/bestsrid_expected @@ -60,7 +60,7 @@ 165|60|999058 171|60|999059 177|60|999060 --180|60|999001 +-180|60|999060 180|60|999060 -177|-60|999101 -171|-60|999102 @@ -122,6 +122,6 @@ 165|-60|999158 171|-60|999159 177|-60|999160 --180|-60|999101 +-180|-60|999160 180|-60|999160 world|999000 -- 2.40.0