I am interested in the maximum width of a polygon, e.g. a lake, in the east-west direction. Bounding boxes will help only in simple polygons but not in complex concave polygons.
Hopefully your lake polygon has a number of points, you can create lines on these points with a azimuth (aspect, geo direction). Choose the length of the lines big enough (ST_MakePoint part), so that you can calculate the shortest line between the two most distant lines.
Here is an example:
The example shows the maximum width of the polygon. I choose ST_ShortestLine (red line) for this approach. ST_MakeLine would increase the value (blue line) and the the endpoint of the line (bottom left) would hit the blue line of the polygon. You have to calculate the distance with the centroids of the created (help)lines.
An idea for irregular or concave polygons for this approach. May be you have to intersect the polygon with a raster.
This likely requires some scripting in any GIS platform.
The most efficient method (asymptotically) is a vertical line sweep: it requires sorting the edges by their minimum y-coordinates and then processing the edges from bottom (minimum y) to top (maximum y), for a O(e * log(e)) algorithm when e edges are involved.
The procedure, although simple, is surprisingly tricky to get right in all cases. Polygons can be nasty: they can have dangles, slivers, holes, be disconnected, have duplicated vertices, runs of vertices along straight lines, and have undissolved boundaries between two adjacent components. Here's an example exhibiting many of these characteristics (and more):
We will specifically seek the maximum-length horizontal segment(s) lying entirely within the closure of the polygon. For instance, this eliminates the dangle between x=20 and x=40 emanating from the hole between x=10 and x=25. It is then straightforward to show that at least one of the maximum-length horizontal segments intersects at least one vertex. (If there are solutions intersecting no vertices, they will lie in the interior of some parallelogram bounded at the top and bottom by solutions which do intersect at least one vertex. This gives us a means to find all solutions.)
Accordingly, the line sweep must begin with the lowest vertices and then move upwards (that is, towards higher y values) to stop at each vertex. At each stop, we find any new edges emanating upwards from that elevation; eliminate any edges terminating from below at that elevation (this is one of the key ideas: it simplifies the algorithm and eliminates half of the potential processing); and carefully process any edges lying wholly at a constant elevation (the horizontal edges).
For example, consider the state when a level of y=10 is reached. From left to right, we find the following edges:
In this table, (x.min, y.min) are coordinates of the lower endpoint of the edge and (x.max, y.max) are coordinates of its upper endpoint. At this level (y=10), the first edge is intercepted within its interior, the second one is intercepted at its bottom, and so on. Some edges terminating at this level, such as from (10,0) to (10,10), are not included in the list.
To determine where the interior points are and the exterior ones are, imagine starting from the extreme left--which is outside the polygon, of course--and moving horizontally to the right. Each time we cross an edge which is not horizontal, we alternately switch from exterior to interior and back. (This is another key idea.) However, all the points within any horizontal edge are determined to be inside the polygon, no matter what. (The closure of a polygon always includes its edges.)
Continuing the example, here is the sorted list of x-coordinates where non-horizontal edges begin at or cross the y=10 line:
(Notice that x=40 is not in this list.) The values of the
Some of these intervals, such as the one from x=60 to x=63.3, do not intersect any vertices: a quick check against the x-coordinates of all vertices with y=10 eliminates such intervals.
During the scan we can monitor the lengths of these intervals, retaining the data concerning the maximum-length interval(s) found so far.
Notice some of the implications of this approach. A "v" shaped vertex, when encountered, is the origin of two edges. Therefore two switches occur when crossing it. Those switches cancel out. Any upside-down "v" is not even processed, because both of its edges are eliminated before starting the left-to-right scan. In both cases, such a vertex does not block off a horizontal segment.
More than two edges can share a vertex: this is illustrated at (10,0), (60,5), (25, 20), and--although it's hard to tell--at (20,10) and (40,10). (That's because the dangle goes (20,10) --> (40,10) --> (40,0) --> (40, -50) --> (40, 10) --> (20,10). Notice how the vertex at (40,0) is also in the interior of another edge...that's nasty.) This algorithm handles those situations just fine.
A tricky situation is illustrated at the very bottom: the x-coordinates of non-horizontal segments there are
This causes everything to the left of x=30 to be considered exterior, everything between 30 and 50 to be interior, and everything after 50 to be exterior again. The vertex at x=40 is never even considered in this algorithm.
Here is what the polygon looks like at the end of the scan. I show all the vertex-containing interior intervals in dark gray, any maximum-length intervals in red, and color the vertices according to their y-coordinates. The maximum interval is 64 units long.
The only geometric calculations involved are to compute where edges intersect horizontal lines: that is a simple linear interpolation. Calculations are also needed to determine which interior segments contain vertices: these are betweenness determinations, easily computed with a couple of inequalities. This simplicity makes the algorithm robust and appropriate both for integer and floating point coordinate representations.
If coordinates are geographic, then horizontal lines are really on circles of latitude. Their lengths are not difficult to compute: just multiply their Euclidean lengths by the cosine of their latitude (in a spherical model). Therefore this algorithm adapts nicely to geographic coordinates. (To handle wrapping around the +-180 meridian well, one might need first to find a curve from the south pole to the north pole that does not pass through the polygon. After re-expressing all x-coordinates as horizontal displacements relative to that curve, this algorithm will correctly find the maximum horizontal segment.)
The following is the
|show 1 more comment|
Here is a raster-based solution. It's fast (I did all the work from start to finish in 14 minutes), requires no scripting, takes only a few operations, and is reasonably accurate.
Start with a raster representation of the polygon. This one uses a grid of 550 rows and 1200 columns:
In this representation, the gray (inside) cells have the value 1 and all other cells are NoData.
Compute the flow accumulation in the west-to-east direction using the unit cell values for the weight grid (amount of "rainfall"):
Low accumulation is dark, increasing to highest accumulations in the bright yellow.
A zonal maximum (using the polygon for the grid and the flow accumulation for the values) identifies the cell(s) where the flow is greatest. To show these, I had to zoom in toward the bottom right:
The red cells mark the ends of the highest accumulations of flow: they are the rightmost endpoints of maximal-length interior segments of the polygon.
To find these segments, place all the weight at the red cells and run the flow backwards!
The red stripe near the bottom marks two rows of cells: within them lies the maximum-length horizontal segment. Use this representation as-is for further analysis or convert it to a polyline (or polygon) shape.
There is some discretization error made with a raster representation. It can be reduced by increasing the resolution, at some cost in computation time.
One really nice aspect of this approach is that typically we find extreme values of things as part of a larger workflow in which some objective needs to be attained: siting a pipeline or football field, creating ecological buffers, and so on. The process involves trade-offs. Thus, the very longest horizontal line might not be part of an optimal solution. We might care instead to know where almost longest lines would lie. This is simple: rather than selecting the zonal maximum flow, select all cells close to a zonal maximum. In this example, the zonal max equals 744 (the number of columns spanned by the longest interior segment). Instead, let's select all cells within 5% of that maximum:
Running the flow from east to west produces this collection of horizontal segments:
This is a map of locations where uninterrupted east-west extent is 95% or greater than the maximal east-west extent anywhere within the polygon.
I'm not sure that Fetzer's answer is what you want to do, but is so the st_box2d may do the job.
SS_Rebelious's idea N°1 will work in many cases but not for some concave polygons.
I think that you must create artificial lw-lines whom points follow edges when vertex-made lines cross the borders of the polygon if there is an east-west line possibility.
For this you could try to make a 4 nodes polygon where line length is high, make the polygon P* which is the preceding overlapping with you original polygon, and see if the min(y1) and max(y2) leave some x-line possibility. (where y1 is the set of point between top left cornet and top right corner and y2 the set of y between bottom left and bottom right corners of your 4 nodes polygon). This is not that easy I hope you will find psql tools to help you!
Ok. I've got another (better) idea (idea-№2). But I suppose that it is better to be realised as a python script, not as an SQL-querry. Again here is the common case, not only E-W.
You will need a bounding box for the polygon and an azimuth (A) as your measurement direction. Assume that the length of the BBox edges are LA and LB. Maximum possible distance (MD) within a polygon is:
I have an idea-№1 (Edit: for common case, not only E-W direction, and with some limitations that are described in comments). I will not provide the code, just a concept. The "x-direction" is actually an azimuth, which is calculated by ST_Azimuth. The proposed steps are:
|show 3 more comments|