Changeset 1d03645f in git

Timestamp:

06/14/24 09:43:57 (11 months ago)

Author:

Paul Ramsey <pramsey@…>

Branches:

[('stable-3.4', '122ad6f692cdd6068b6fc8b2b9532a9204572fa0')]

Children:

['682776d3f96ab2400e662530b03f99a640fe0e19']

Parents:

['bb73a3fea9cf995da73ef63226c693849584c081']

Message:

Support geography estimated extents, references #5734

Files:

• NEWS (modified) (1 diff)
• liblwgeom/lwgeodetic.c (modified) (1 diff)
• liblwgeom/lwgeodetic.h (modified) (1 diff)
• postgis/gserialized_estimate.c (modified) (7 diffs)

NEWS

@@ -30 +30 @@
           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

liblwgeom/lwgeodetic.c

@@ -3580 +3580 @@
         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;
+}
+

liblwgeom/lwgeodetic.h

@@ -153 +153 @@
 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 */

postgis/gserialized_estimate.c

@@ -117 +117 @@
 #include "lwgeom_pg.h"       /* For debugging macros. */
 #include "gserialized_gist.h" /* For index common functions */
+#include "lwgeodetic.h"
 
 #include <math.h>
@@ -2285 +2286 @@
 
 
+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);
+}
+
 
 /**
@@ -2304 +2317 @@
         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() */
@@ -2347 +2362 @@
         }
 
+        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)
@@ -2364 +2382 @@
 
                 /* 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 */
@@ -2373 +2391 @@
 
                 /* 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);
@@ -2557 +2590 @@
         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