3

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.

Gregory
  • 1,661
  • 4
  • 20
  • 30
  • This problem is solved with an algorithm like the one described at http://gis.stackexchange.com/questions/5300/how-can-i-can-i-divide-a-polygon-into-specific-sizes-using-arcgis/5302#5302. Use a vertical sweep direction and a cylindrical equal-area projection. – whuber May 02 '12 at 19:50
  • @whuber Thanks for the link to your answer for a similar question. That approach looks great and I'm looking forward to trying it, unfortunately a (figurative) fire at work may keep me from it for as much couple weeks. I'll post back here when I've figured it out. – Gregory May 02 '12 at 20:12

1 Answers1

3

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:

QGIS screenshot

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.

diciu
  • 2,348
  • 16
  • 15
  • 1
    There is a good analysis and some nice ideas here, including a clear presentation of potential limitations. Notice, please, that to obtain what the OP asks for, you must be sure to work in a cylindrical projection (so that the horizontal slices are lines of latitude) and, for features near the poles or extending more than a few degrees north-south, it's important to use an equal-area projection for correct computation of areas. – whuber May 02 '12 at 19:54