How to query latitudinal bounds on 95% of a polygon in PostGIS?

Question

I'm going to say thanks in advance. I'm very new to GIS and this forum is a tremendous help. My simple questions have been solved; here's one that might be a little trickier.

Given a polygon in a PostGIS geometry table, I want to find the northern and southern latitudes that would bound an arbitrary portion of the polygon's area (say 95%), with the remainder trimmed equally from the north and the south. I think what I want may be called a bounding box, though for a large polygon it might seem more like a conic section. I've noticed a number of functions relating to boxes in the documentation, but haven't figured out how they could be defined as I've described.

Answer 1

If you're willing to work with the geometry's bounding box, here's a possible solution:

First get the ymin, ymax and latitude span:

osm=# select st_ymax(st_envelope(way)), st_ymin(st_envelope(way)), st_ymax(st_envelope(way)) - st_ymin(st_envelope(way)) as ylength from planet_osm_polygon where name = 'Sector 1';
  st_ymax   |  st_ymin   | ylength  
------------+------------+----------
 5549655.89 | 5532946.68 | 16709.21
(1 row)

Compute ymax, ymin for the bbox covering 95% of the original bbox area (but not 95% of the geometry's area):

osm=# select 5549655.89 - ((16709.21 - 16709.21 * .95)/2) as ymax;
           ymax           
--------------------------
 5549238.1597500000000000
(1 row)

osm=# select 5532946.68 + ((16709.21 - 16709.21 * .95)/2) as ymin;
           ymin           
--------------------------
 5533364.4102500000000000

(1 row)

Make a box from the latitudes of the smaller box and the longitudes of the original box:

osm=# select st_makebox2d(st_point(2892122.63, 5533364.4102500000000000), st_point(2907085.65, 5549238.1597500000000000));
                 st_makebox2d                 
----------------------------------------------
 BOX(2892122.75 5533364.5,2907085.75 5549238)
(1 row)

Intersect the original geometry with the smaller bbox:

osm=# select st_intersection(way, st_setsrid(st_makebox2d(st_point(2892122.63, 5533364.4102500000000000), st_point(2907085.65, 5549238.1597500000000000)), 900913))  from planet_osm_polygon where name = 'Sector 1';

Here's how the geometries look like on the same data I've used:

Note that, as expected, the area is not 95% of the original (it's more like 99).

I've no idea if there's an algorithm out there to address the original requirement of controlling area span.

I think you could try something like slicing the polygon in, say 100 pieces (~100 latitude points spanning the latitude of the poly).

To compute the area span of a slice you need to compute the area of the intersection of the slice bbox and the original geometry.

Then you could select those slices whose areas sum close to the area span and aggregate them back to obtain your target poly.