# How to calculate the maximum distance within a polygon in x-direction (east-west direction) in PostGIS?

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.

-
I'm not familiar with postgis functionality. However, there may be a bounding box tool. The width of the bounding box would be the maximum distance in the E-W direction. –  Fezter Aug 30 '12 at 22:48
@Fetzter that's not correct: a counterexample, even for a simple complex polygon, is a thin rhombus extending from SW to NE. Its maximum east-west width can be an arbitrarily small fraction of the width of its bounding box. –  whuber Aug 31 '12 at 16:24
I have created an utility for this task based on this and this proposals. It is able to calculate max or min width of the polygon. Currently it works with a shp-files, but you may rewrite it to work with a PostGIS or just wait for some time till it evolve to a QGIS plugin that will work with PostGIS too. Detailed description and download link here. –  SS_Rebelious Sep 21 '12 at 15:29

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:

      x.min x.max y.min y.max
[1,]    10     0     0    30
[2,]    10    24    10    20
[3,]    20    24    10    20
[4,]    20    40    10    10
[5,]    40    20    10    10
[6,]    60     0     5    30
[7,]    60    60     5    30
[8,]    60    70     5    20
[9,]    60    70     5    15
[10,]    90   100    10    40


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:

x.array    6.7 10 20 48 60 63.3 65 90
interior     1  0  1  0  1    0  1  0


(Notice that x=40 is not in this list.) The values of the interior array mark the left endpoints of interior segments: 1 designates an interior interval, 0 an exterior interval. Thus, the first 1 indicates the interval from x=6.7 to x=10 is inside the polygon. The next 0 indicates the interval from x=10 to x=20 is outside the polygon. And so it proceeds: the array identifies four separate intervals as inside the polygon.

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

30, 50


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 R code implemented to perform the calculations and create the illustrations.

#
# Plotting functions.
#
points.polygon <- function(p, ...) {
points(p$v, ...) } plot.polygon <- function(p, ...) { apply(p$e, 1, function(e) lines(matrix(e[c("x.min", "x.max", "y.min", "y.max")], ncol=2), ...))
}
expand <- function(bb, e=1) {
a <- matrix(c(e, 0, 0, e), ncol=2)
origin <- apply(bb, 2, mean)
delta <-  origin %*% a - origin
t(apply(bb %*% a, 1, function(x) x - delta))
}
#
# Convert polygon to a better data structure.
#
# A polygon class has three attributes:
#   v is an array of vertex coordinates "x" and "y" sorted by increasing y;
#   e is an array of edges from (x.min, y.min) to (x.max, y.max) with y.max >= y.min, sorted by y.min;
#   bb is its rectangular extent (x0,y0), (x1,y1).
#
as.polygon <- function(p) {
#
# p is a list of linestrings, each represented as a sequence of 2-vectors
# with coordinates in columns "x" and "y".
#
f <- function(p) {
g <- function(i) {
v <- p[(i-1):i, ]
v[order(v[, "y"]), ]
}
sapply(2:nrow(p), g)
}
vertices <- do.call(rbind, p)
edges <- t(do.call(cbind, lapply(p, f)))
colnames(edges) <- c("x.min", "x.max", "y.min", "y.max")
#
# Sort by y.min.
#
vertices <- vertices[order(vertices[, "y"]), ]
vertices <- vertices[!duplicated(vertices), ]
edges <- edges[order(edges[, "y.min"]), ]

# Maintaining an extent is useful.
bb <- apply(vertices <- vertices[, c("x","y")], 2, function(z) c(min(z), max(z)))

# Package the output.
l <- list(v=vertices, e=edges, bb=bb); class(l) <- "polygon"
l
}
#
# Compute the maximal horizontal interior segments of a polygon.
#
fetch.x <- function(p) {
#
# Update moves the line from the previous level to a new, higher level, changing the
# state to represent all edges originating or strictly passing through level y.
#
update <- function(y) {
if (y > state$level) { state$level <<- y
#
# Remove edges below the new level from state$current. # current <- state$current
current <- current[current[, "y.max"] > y, ]
#
# Adjoin edges at this level.
#
i <- state$i while (i <= nrow(p$e) && p$e[i, "y.min"] <= y) { current <- rbind(current, p$e[i, ])
i <- i+1
}
state$i <<- i # # Sort the current edges by x-coordinate. # x.coord <- function(e, y) { if (e["y.max"] > e["y.min"]) { ((y - e["y.min"]) * e["x.max"] + (e["y.max"] - y) * e["x.min"]) / (e["y.max"] - e["y.min"]) } else { min(e["x.min"], e["x.max"]) } } if (length(current) > 0) { x.array <- apply(current, 1, function(e) x.coord(e, y)) i.x <- order(x.array) current <- current[i.x, ] x.array <- x.array[i.x] # # Scan and mark each interval as interior or exterior. # status <- FALSE interior <- numeric(length(x.array)) for (i in 1:length(x.array)) { if (current[i, "y.max"] == y) { interior[i] <- TRUE } else { status <- !status interior[i] <- status } } # # Simplify the data structure by retaining the last value of interior # within each group of common values of x.array. # interior <- sapply(split(interior, x.array), function(i) rev(i)[1]) x.array <- sapply(split(x.array, x.array), function(i) i[1]) print(y) print(current) print(rbind(x.array, interior)) markers <- c(1, diff(interior)) intervals <- x.array[markers != 0] # # Break into a list structure. # if (length(intervals) > 1) { if (length(intervals) %% 2 == 1) intervals <- intervals[-length(intervals)] blocks <- 1:length(intervals) - 1 blocks <- blocks - (blocks %% 2) intervals <- split(intervals, blocks) } else { intervals <- list() } } else { intervals <- list() } # # Update the state. # state$current <<- current
}
list(y=y, x=intervals)
} # Update()

process <- function(intervals, x, y) {
# intervals is a list of 2-vectors. Each represents the endpoints of
# an interior interval of a polygon.
# x is an array of x-coordinates of vertices.
#
# Retains only the intervals containing at least one vertex.
between <- function(i) {
1 == max(mapply(function(a,b) a && b, i[1] <= x, x <= i[2]))
}
is.good <- lapply(intervals$x, between) list(y=y, x=intervals$x[unlist(is.good)])
#intervals
}
#
# Group the vertices by common y-coordinate.
#
vertices.x <- split(p$v[, "x"], p$v[, "y"])
vertices.y <- lapply(split(p$v[, "y"], p$v[, "y"]), max)
#
# The "state" is a collection of segments and an index into edges.
# It will updated during the vertical line sweep.
#
state <- list(level=-Inf, current=c(), i=1, x=c(), interior=c())
#
# Sweep vertically from bottom to top, processing the intersection
# as we go.
#
mapply(function(x,y) process(update(y), x, y), vertices.x, vertices.y)
}

scale <- 10
p.raw = list(scale * cbind(x=c(0:10,7,6,0), y=c(3,0,0,-1,-1,-1,0,-0.5,0.75,1,4,1.5,0.5,3)),
scale *cbind(x=c(1,1,2.4,2,4,4,4,4,2,1), y=c(0,1,2,1,1,0,-0.5,1,1,0)),
scale *cbind(x=c(6,7,6,6), y=c(.5,2,3,.5)))

#p.raw = list(cbind(x=c(0,2,1,1/2,0), y=c(0,0,2,1,0)))
#p.raw = list(cbind(x=c(0, 35, 100, 65, 0), y=c(0, 50, 100, 50, 0)))

p <- as.polygon(p.raw)

results <- fetch.x(p)
#
# Find the longest.
#
dx <- matrix(unlist(results["x", ]), nrow=2)
length.max <- max(dx[2,] - dx[1,])
#
# Draw pictures.
#
segment.plot <- function(s, length.max, colors,  ...) {
lapply(s$x, function(x) { col <- ifelse (diff(x) >= length.max, colors[1], colors[2]) lines(x, rep(s$y,2), col=col, ...)
})
}
gray <- "#f0f0f0"
grayer <- "#d0d0d0"
plot(expand(p$bb, 1.1), type="n", xlab="x", ylab="y", main="After the Scan") sapply(1:length(p.raw), function(i) polygon(p.raw[[i]], col=c(gray, "White", grayer)[i])) apply(results, 2, function(s) segment.plot(s, length.max, colors=c("Red", "#b8b8a8"), lwd=4)) plot(p, col="Black", lty=3) points(p, pch=19, col=round(2 + 2*p$v[, "y"]/scale, 0))
points(p, cex=1.25)

-
Is there a theorem that proves that maximum length line inside non-convex polygon in any given direction intersects at least one vertex of this polygon? –  SS_Rebelious Sep 14 '12 at 11:57
@SS Yes, there is. Here's a sketch of a proof: if there is no intersection, then the endpoints of the segment both lie on the interiors of edges and the segment can be moved, at least a little, up and down. Its length is a linear function of the amount of displacement. Thus, it can have a maximum length only if the length does not change when being moved. This implies both (a) it is part of a parallelogram formed of maximal-length segments and (b) both the top and bottom edges of that parallelogram must meet a vertex, QED. –  whuber Sep 14 '12 at 14:43
And what is the name of this theorem? I'm struggling to find it. BTW, what about curved edges that has no vertex (I mean a theoretical approach)? A sketch of example of the figure I mean (a dumb-bell shaped polygon): "C=D". –  SS_Rebelious Sep 14 '12 at 15:05
@SS When the edges are curved, the theorem no longer holds. Techniques of differential geometry can be applied to obtain useful results. I learned these methods from Cheeger & Ebin's book, Comparison Theorems in Riemannian Geometry. However, most GISes will approximate curves by detailed polylines anyway, so the question (as a practical matter) is moot. –  whuber Sep 14 '12 at 15:40
could you specify the name of the theorem (and the page if possible)? I've got the book and I wasn't able to locate the needed theorem. –  SS_Rebelious Sep 14 '12 at 16:28

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!

-
This is on the right track. The longest E-W segment will be found among the intersections with the polygon's interior with horizontal lines passing through the vertices of the polygon. This requires code to loop over the vertices. There's an alternative (but equivalent) method available by following an artificial east-west flow across a raster representation of the polygon: the maximum flow length found in the polygon (which is one of its "zonal statistics") is the desired width. The raster solution is obtained in just 3 or 4 steps and requires no looping or scripting. –  whuber Aug 31 '12 at 16:37
@Aname, please add "№1" to the "SS_Rebelious's idea" to avoid misunderstanding: I've added another proposal. I cannot edit your answer myself because this edit is less then 6 characters. –  SS_Rebelious Aug 31 '12 at 20:21

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: MB = (LA^2 * LB^2)^(1/2), so seeking value (V) is not bigger than MB: V <= MB.

1. Starting from any vertex of the BBox create a line (LL) with the length MB and azimuth A.
2. Intersect line LL with the polygon to get intersection line (IL)
3. Check geometry of IL - if there only two points in IL line then calculate its length. If 4 or more - calculate segments and pick up the length of the longest one. Null (no intersection at all) - ignore.
4. Keep creating another LL lines moving from the starting point counter or clockwise towards the edges of the BBox untill you won't end up at starting point (you will do the whole loop across BBox).
5. Pick up the biggest IL length value (actually you dont't have to store all lengths, you can keep just 'so far' maximum value while looping) - it will be what you seek.
-
This sounds like a double loop over the vertices: that is inefficient enough that it should be avoided (except for very simplified polygons). –  whuber Aug 21 '13 at 20:56
@whuber, I don't see any extra loops here. There is just some meaningless processing of 2 sides of BB that will give nothing but nulls. But this processing was excluded in the script that I provided in the answer that was deleted here (it seems that it is a comment now, but I don't see it as a comment - only as a deleted answer) –  SS_Rebelious Aug 22 '13 at 9:44
(1) It is the third comment to the question. (2) You're right: upon reading your description very carefully, it appears to me now that you are finding the longest segment between (the four) vertices of the bounding box and vertices of the polygon. I don't see how this answers the question, though: the result definitely is not what the O.P. sought. –  whuber Aug 22 '13 at 11:51
@whuber, proposed algorithm finds the longest intersection of the polygon with the line representing the given direction. Apparently the result IS what was asked if the distance between intersection lines -> 0 or it passes all the vertexes (for not curved figures). –  SS_Rebelious Aug 22 '13 at 13:38

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:

1. Extract all vertices from polygon as points.
2. Create lines between every pair of points.
3. Select lines (lets call them lw-lines) that are within the original polygon (we don't need lines that will cross the borders of the polygon).
4. Find distances and azimuths for every lw-line.
5. Select longest distance from lw-lines where azimuth is equal to sought-for azimuth or lies in some interval (it might be that no azimuth will be exactly equal to sought-for azimuth).
-
This will fail to work even for some triangles, such as the one at vertices (0,0), (1000, 1000), and (501, 499). Its maximum east-west width is approximately 2; the azimuths are all around 45 degrees; and regardless, the shortest line segment between vertices is over 350 times longer than the east-west width. –  whuber Aug 31 '12 at 16:26
@whuber, you are right, it will fail for triangals, but for polygons, representing some nature features it may be useful. –  SS_Rebelious Aug 31 '12 at 17:45
It's hard to recommend a procedure that fails dramatically even for simple cases on the hope that it might sometimes get a correct answer! –  whuber Aug 31 '12 at 19:10
@whuber, so don't recommend it! ;-) I proposed this workaround because there were no answers to this question. Notice that you may post your own better answer. BTW, if you will place some points on the triangle edges, my proposal will work ;-) –  SS_Rebelious Aug 31 '12 at 19:19
I have suggested several approaches. A raster one is at gis.stackexchange.com/questions/32552/… and elaborated at forums.esri.com/Thread.asp?c=93&f=982&t=107703&mc=3. Another one--not quite as applicable, but with interesting uses--is at gis.stackexchange.com/questions/23664/… (the radon transform). This is illustrated at stats.stackexchange.com/a/33102. –  whuber Aug 31 '12 at 20:49

Have a look at my question and the answer from Evil Genius.

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.

-