Calculating maximum distance within polygon in x-direction (east-west direction) in PostGIS?

Question

I am interested in the maximum width of a polygon e.g. a lake, in east-west direction. Bounding boxes will help only in simple polygons but not in complex concave polygons.

Answer 1

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)

Answer 2

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.

Answer 3

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.

Answer 4

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!

Answer 5

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).

Answer 6

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.

Answer 7

So, so far, my solution to your question looks like a new experimental spatial SQL function.

1. Create a function called ST_MaxDistanceWestEastIRRegularPolygon using pgAdmin:

CREATE OR REPLACE FUNCTION ST_MaxDistanceWestEastIRRegularPolygon(
geom GEOMETRY,
n integer)
RETURNS GEOMETRY AS
$BODY$
WITH
tbl_rigth AS (WITH 
    tbla AS (SELECT (ST_Dump(geom)).geom geom),
    tblb AS (SELECT ST_Boundary(ST_Envelope(geom)) geom FROM tbla),
    tblc AS (SELECT line1, line2 FROM (SELECT ST_MakeLine(ST_PointN(geom,1), ST_PointN(geom,2)) line1, 
             ST_MakeLine(ST_PointN(geom,4), ST_PointN(geom,3)) line2 FROM tblb) foo),
    tbld AS (SELECT generate_series (0, n) as steps),
    tble AS (SELECT steps AS stp1, ST_LineInterpolatePoint(line1, steps/(SELECT count(steps)::float-1 FROM tbld)) geom1 FROM tblc, tbld GROUP BY tbld.steps, geom1),
    tblf AS (SELECT steps AS stp2, ST_LineInterpolatePoint(line2, steps/(SELECT count(steps)::float-1 FROM tbld)) geom2 FROM tblc, tbld GROUP BY tbld.steps, geom2),
    tblg AS (SELECT ST_MakeLine(geom1, geom2) geom FROM tble JOIN tblf ON true AND stp1=stp2),
    tblh AS (SELECT (ST_Dump(ST_Intersection(a.geom, b.geom))).geom geom FROM tblg a JOIN tbla b ON true GROUP BY a.geom, b.geom),
    tbli AS (SELECT ST_Length(geography(geom)) lenght, geom FROM tblh),
    tblj AS (SELECT geom, MAX(lenght) max_length, degrees(ST_Azimuth(ST_StartPoint(geom), ST_EndPoint(geom))) azimuth FROM tbli GROUP BY geom ORDER BY max_length DESC LIMIT 1)
            SELECT geom FROM tblj WHERE azimuth IN (90,-90)),
tbl_lefth AS (WITH 
    tbla AS (SELECT (ST_Dump(geom)).geom geom),
    tblb AS (SELECT ST_Boundary(ST_Envelope(geom)) geom FROM tbla),
    tblc AS (SELECT line1, line2 FROM (SELECT ST_MakeLine(ST_PointN(geom,2), ST_PointN(geom,3)) line1, 
             ST_MakeLine(ST_PointN(geom,1), ST_PointN(geom,4)) line2 FROM tblb) foo),
    tbld AS (SELECT generate_series (0, n) as steps),
    tble AS (SELECT steps AS stp1, ST_LineInterpolatePoint(line1, steps/(SELECT count(steps)::float-1 FROM tbld)) geom1 FROM tblc, tbld GROUP BY tbld.steps, geom1),
    tblf AS (SELECT steps AS stp2, ST_LineInterpolatePoint(line2, steps/(SELECT count(steps)::float-1 FROM tbld)) geom2 FROM tblc, tbld GROUP BY tbld.steps, geom2),
    tblg AS (SELECT ST_MakeLine(geom1, geom2) geom FROM tble JOIN tblf ON true AND stp1=stp2),
    tblh AS (SELECT (ST_Dump(ST_Intersection(a.geom, b.geom))).geom geom FROM tblg a JOIN tbla b ON true GROUP BY a.geom, b.geom),
    tbli AS (SELECT ST_Length(geography(geom)) lenght, geom FROM tblh),
    tblj AS (SELECT geom, MAX(lenght) max_length, degrees(ST_Azimuth(ST_StartPoint(geom), ST_EndPoint(geom))) azimuth FROM tbli GROUP BY geom ORDER BY max_length DESC LIMIT 1)
            SELECT geom FROM tblj WHERE azimuth IN (90,-90))
            SELECT geom FROM tbl_rigth UNION SELECT geom FROM tbl_lefth
    $BODY$
LANGUAGE SQL

2. Run the SQL script:

WITH 
tbla AS (SELECT ST_MaxDistanceWestEastIRRegularPolygon(geom,300) geom FROM <table_name>) 
SELECT ST_Distance(ST_StartPoint(geom)::geography, ST_EndPoint(geom)::geography) distance FROM tbla

Explanation:

I specifically implemented a geometric vector solution, so the distance is defined this way.

Remember that increasing the point density will improve the result, but will have a negative impact on performance.

So far, the function does well with a small number of polygons, but with the increase in performance of the PC, this fact will not be an obstacle in the future :-)...

This construction allows to find not only the longest but also the shortest horizontal line and its length, as well as the longest and shortest vertical line, and it also works for undirected complex polygons (with function - ST_OrientedEnvelope())- but that is another question :-)...

Original spatial solutions...

Translated with www.DeepL.com/Translator (free version)