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);
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.
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};
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)
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 */
/** 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 */
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
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);
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
*/
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);
}
*/
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 );
}
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? */
165|60|999058
171|60|999059
177|60|999060
--180|60|999001
+-180|60|999060
180|60|999060
-177|-60|999101
-171|-60|999102
165|-60|999158
171|-60|999159
177|-60|999160
--180|-60|999101
+-180|-60|999160
180|-60|999160
world|999000
projects / postgis / commitdiff