[SCM] PostGIS branch stable-3.4 updated. 3.4.2-42-g1d03645f2
git at osgeo.org
git at osgeo.org
Fri Jun 14 09:44:14 PDT 2024
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, stable-3.4 has been updated
via 1d03645f2643e6d2b49861b3eecb4f957247a4e0 (commit)
from bb73a3fea9cf995da73ef63226c693849584c081 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 1d03645f2643e6d2b49861b3eecb4f957247a4e0
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Fri Jun 14 09:43:57 2024 -0700
Support geography estimated extents, references #5734
diff --git a/NEWS b/NEWS
index ad811061a..bded777e6 100644
--- a/NEWS
+++ b/NEWS
@@ -29,6 +29,7 @@ To take advantage of all SFCGAL featurs, SFCGAL 1.4.1+ is needed.
- #5653, Do not simplify away points when linestring doubles back on
itself (Paul Ramsey)
- #5720, Correctly mangle special column names in shp2pgsql (Paul Ramsey)
+ - #5734, Estimate geography extent more correctly (Paul Ramsey)
PostGIS 3.4.2
2024/02/08
diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 29eef8e95..176c5ed2c 100644
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -3579,3 +3579,135 @@ int ptarray_contains_point_sphere(const POINTARRAY *pa, const POINT2D *pt_outsid
return LW_FALSE;
}
+
+
+
+/*
+* Given a geodetic bounding volume, calculate a lon/lat bounding
+* box that should contain the original feature that gave rise to
+* the geodetic box, in plate-carre space (planar lon/lat).
+*/
+int gbox_geocentric_get_gbox_cartesian(const GBOX *gbox_geocentric, GBOX *gbox_planar)
+{
+ /* Normalized corners of the bounding volume */
+ POINT3D corners[8];
+ POINT3D cap_center = {0,0,0};
+ double furthest_angle = 0.0;
+ double cap_angle = 0.0;
+ int all_longitudes = LW_FALSE;
+ double lon0 = -M_PI, lon1 = M_PI;
+ double lat0, lat1;
+ GEOGRAPHIC_POINT cap_center_g;
+
+ if (!gbox_geocentric || !gbox_planar)
+ {
+ lwerror("Null pointer passed to %s", __func__);
+ return LW_FALSE;
+ }
+
+#define CORNER_SET(ii, xx, yy, zz) { \
+ corners[ii].x = gbox_geocentric->xx; \
+ corners[ii].y = gbox_geocentric->yy; \
+ corners[ii].z = gbox_geocentric->zz; \
+ }
+
+ /*
+ * First find a "centered" vector to serve as the mid-point
+ * of the input bounding volume.
+ */
+ CORNER_SET(0, xmin, ymin, zmin);
+ CORNER_SET(1, xmax, ymin, zmin);
+ CORNER_SET(2, xmin, ymax, zmin);
+ CORNER_SET(3, xmax, ymax, zmin);
+ CORNER_SET(4, xmin, ymin, zmax);
+ CORNER_SET(5, xmax, ymin, zmax);
+ CORNER_SET(6, xmin, ymax, zmax);
+ CORNER_SET(7, xmax, ymax, zmax);
+
+ /*
+ * Normalize the volume corners
+ * and normalize the final vector.
+ */
+ for (uint32_t i = 0; i < 8; i++)
+ {
+ normalize(&(corners[i]));
+ cap_center.x += corners[i].x;
+ cap_center.y += corners[i].y;
+ cap_center.z += corners[i].z;
+ }
+ normalize(&cap_center);
+
+ /*
+ * Find the volume corner that is furthest from the center,
+ * and calculate the angle between the center and the corner.
+ * Now we have a "cap" (center and angle)
+ */
+ for (uint32_t i = 0; i < 8; i++)
+ {
+ double angle = vector_angle(&cap_center, &(corners[i]));
+ if (angle > furthest_angle)
+ furthest_angle = angle;
+ }
+ cap_angle = furthest_angle;
+
+ /*
+ * Calculate the planar box that contains the cap.
+ * If the cap contains a pole, then we include all longitudes
+ */
+ cart2geog(&cap_center, &cap_center_g);
+
+ /* Check whether cap includes the south pole */
+ lat0 = cap_center_g.lat - cap_angle;
+ if (lat0 <= -M_PI_2)
+ {
+ lat0 = -M_PI_2;
+ all_longitudes = LW_TRUE;
+ }
+
+ /* Check whether cap includes the north pole */
+ lat1 = cap_center_g.lat + cap_angle;
+ if (lat1 >= M_PI_2)
+ {
+ lat1 = M_PI_2;
+ all_longitudes = LW_TRUE;
+ }
+
+ if (!all_longitudes)
+ {
+ // Compute the range of longitudes covered by the cap. We use the law
+ // of sines for spherical triangles. Consider the triangle ABC where
+ // A is the north pole, B is the center of the cap, and C is the point
+ // of tangency between the cap boundary and a line of longitude. Then
+ // C is a right angle, and letting a,b,c denote the sides opposite A,B,C,
+ // we have sin(a)/sin(A) = sin(c)/sin(C), or sin(A) = sin(a)/sin(c).
+ // Here "a" is the cap angle, and "c" is the colatitude (90 degrees
+ // minus the latitude). This formula also works for negative latitudes.
+ //
+ // The formula for sin(a) follows from the relationship h = 1 - cos(a).
+
+ double sin_a = sin(cap_angle);
+ double sin_c = cos(cap_center_g.lat);
+ if (sin_a <= sin_c)
+ {
+ double angle_A = asin(sin_a / sin_c);
+ lon0 = remainder(cap_center_g.lon - angle_A, 2 * M_PI);
+ lon1 = remainder(cap_center_g.lon + angle_A, 2 * M_PI);
+ }
+ else
+ {
+ lon0 = -M_PI;
+ lon1 = M_PI;
+ }
+ }
+
+ gbox_planar->xmin = rad2deg(lon0);
+ gbox_planar->ymin = rad2deg(lat0);
+ gbox_planar->xmax = rad2deg(lon1);
+ gbox_planar->ymax = rad2deg(lat1);
+ FLAGS_SET_GEODETIC(gbox_planar->flags, 0);
+ FLAGS_SET_Z(gbox_planar->flags, 0);
+ FLAGS_SET_M(gbox_planar->flags, 0);
+
+ return LW_TRUE;
+}
+
diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h
index 39053ef5f..7636e6ce6 100644
--- a/liblwgeom/lwgeodetic.h
+++ b/liblwgeom/lwgeodetic.h
@@ -152,6 +152,7 @@ double spheroid_distance(const GEOGRAPHIC_POINT *a, const GEOGRAPHIC_POINT *b, c
double spheroid_direction(const GEOGRAPHIC_POINT *r, const GEOGRAPHIC_POINT *s, const SPHEROID *spheroid);
int spheroid_project(const GEOGRAPHIC_POINT *r, const SPHEROID *spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g);
+int gbox_geocentric_get_gbox_cartesian(const GBOX *gbox_geocentric, GBOX *gbox_planar);
#endif /* _LWGEODETIC_H */
diff --git a/postgis/gserialized_estimate.c b/postgis/gserialized_estimate.c
index 40aa85d65..c36564380 100644
--- a/postgis/gserialized_estimate.c
+++ b/postgis/gserialized_estimate.c
@@ -116,6 +116,7 @@ dimensionality cases. (2D geometry) &&& (3D column), etc.
#include "liblwgeom.h"
#include "lwgeom_pg.h" /* For debugging macros. */
#include "gserialized_gist.h" /* For index common functions */
+#include "lwgeodetic.h"
#include <math.h>
#if HAVE_IEEEFP_H
@@ -2284,6 +2285,18 @@ Datum gserialized_gist_sel(PG_FUNCTION_ARGS)
}
+static Oid
+get_attype_by_name(Oid tbl, const char *colname)
+{
+ AttrNumber attnum = get_attnum(tbl, colname);
+ if (attnum == InvalidAttrNumber)
+ elog(ERROR, "cannot find column \"%s\" in table \"%s\"",
+ colname,
+ get_rel_name(tbl));
+
+ return get_atttype(tbl, attnum);
+}
+
/**
* Return the estimated extent of the table
@@ -2303,6 +2316,8 @@ Datum gserialized_estimated_extent(PG_FUNCTION_ARGS)
bool only_parent = false;
int key_type, att_num;
size_t sz;
+ Oid atttypid;
+ bool is_geography = false;
/* We need to initialize the internal cache to access it later via postgis_oid() */
postgis_initialize_cache();
@@ -2346,8 +2361,11 @@ Datum gserialized_estimated_extent(PG_FUNCTION_ARGS)
PG_RETURN_NULL();
}
- /* Read the extent from the head of the spatial index, if there is one */
+ atttypid = get_attype_by_name(tbl_oid, text_to_cstring(col));
+ if (atttypid == postgis_oid(GEOGRAPHYOID))
+ is_geography = true;
+ /* Read the extent from the head of the spatial index, if there is one */
idx_oid = table_get_spatial_index(tbl_oid, col, &key_type, &att_num);
if (idx_oid)
{
@@ -2363,7 +2381,7 @@ Datum gserialized_estimated_extent(PG_FUNCTION_ARGS)
/* Fall back to reading the stats, if no index is found */
/* Estimated extent only returns 2D bounds, so use mode 2 */
- nd_stats = pg_get_nd_stats_by_name(tbl_oid, col, 2, only_parent);
+ nd_stats = pg_get_nd_stats_by_name(tbl_oid, col, is_geography ? 3 : 2, only_parent);
/* Error out on no stats */
if ( ! nd_stats ) {
@@ -2372,17 +2390,32 @@ Datum gserialized_estimated_extent(PG_FUNCTION_ARGS)
}
/* Construct the box */
- gbox = palloc(sizeof(GBOX));
- FLAGS_SET_GEODETIC(gbox->flags, 0);
- FLAGS_SET_Z(gbox->flags, 0);
- FLAGS_SET_M(gbox->flags, 0);
+ gbox = gbox_new(0);
gbox->xmin = nd_stats->extent.min[0];
gbox->xmax = nd_stats->extent.max[0];
gbox->ymin = nd_stats->extent.min[1];
gbox->ymax = nd_stats->extent.max[1];
+ if (is_geography)
+ {
+ FLAGS_SET_Z(gbox->flags, 1);
+ gbox->zmin = nd_stats->extent.min[2];
+ gbox->zmax = nd_stats->extent.max[2];
+ }
+
pfree(nd_stats);
}
+ /* Convert geocentric geography box into a planar box */
+ /* that users understand */
+ if (is_geography)
+ {
+ GBOX *gbox_planar = gbox_new(0);
+ gbox_geocentric_get_gbox_cartesian(gbox, gbox_planar);
+ PG_RETURN_POINTER(gbox_planar);
+ }
+ else
+ PG_RETURN_POINTER(gbox);
+
PG_RETURN_POINTER(gbox);
}
@@ -2556,10 +2589,13 @@ spatial_index_read_extent(Oid idx_oid, int key_type, int att_num)
}
else if (key_type == STATISTIC_KIND_ND && bounds_gidx)
{
+ lwflags_t flags = 0;
if (gidx_is_unknown(bounds_gidx))
return NULL;
- gbox = gbox_new(0);
- gbox_from_gidx(bounds_gidx, gbox, 0);
+ FLAGS_SET_Z(flags, GIDX_NDIMS(bounds_gidx) > 2);
+ FLAGS_SET_M(flags, GIDX_NDIMS(bounds_gidx) > 3);
+ gbox = gbox_new(flags);
+ gbox_from_gidx(bounds_gidx, gbox, flags);
}
else
return NULL;
-----------------------------------------------------------------------
Summary of changes:
NEWS | 1 +
liblwgeom/lwgeodetic.c | 132 +++++++++++++++++++++++++++++++++++++++++
liblwgeom/lwgeodetic.h | 1 +
postgis/gserialized_estimate.c | 52 +++++++++++++---
4 files changed, 178 insertions(+), 8 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list