Intro

In a previous post, we've seen how PostGIS and Openstreetmap can be used to leverage geographical data.

A common task in GIS(in particular land use planning) is land subdivision. This involves taking different shapes, polygons usually, and cutting them into smaller polygons.

This post will focus on describing an algorithm for partitioning a land polygon into parcels of a given area. Parcels will be cut off from the initial polygon until no more parcels can be formed.

This is part2 of a series:

Algorithm description

We have a polygon P and the nearest road R. We get the bounding box B for P and all our searches for cutting points/lines will be confined to B. We compute the extreme points on the bounding box and we label them with the cardinal directions they represent.

Our goal is to cut a corner C (also called subdivision) from P such that it contains the nearest boundary point to a road. The cut will be done using two lines, one horizontal and one vertical. We want C to be of a given area A.

Now in order to find the corner to cut, we look at the extreme points and check which one of them is closest to a road.

We'll use sweeping-lines to find the parcel we need. Any sweeping-lines mentioned will be moving away from the corner (in other words, away from the closest road point).

In what follows, we assume the north-west corner needs to be cut.

We place an inset (a horizontal line) that will be located sqrt(A) to the south (relative to the north edge). The inset is positioned there because we anticipate the target area to be in the form of a square.

If the area above the inset (the one we aim for) is larger than our target, we split the polygon, take the upper half and use another sweeping line that goes from west to east, to find another cutting line that allows us to find a parcel with the required area.

If the area above the inset is insufficient (below the target area), we search for a better position for it, using binary search, along the north-south direction.

Additional details: The way the cut search works, using the inset, is such that we avoid getting thin horizontal strips when our initial polygon is a square/rectangle (and it is expected to be a square in the vast majority of cases).

Details about corner cases (other than NW which was covered above):

  • NE corner: horizontal goes north->south and vertical goes east->west
  • SE corner: horizontal goes south->north and vertical goes east->west
  • SW corner: horizontal goes south->north and vertical goes west->east

So the sweep lines always move away from the corner.

After the parcel with target area was found, it will be cut off from the polygon, the original polygon in the GIS database will be updated and the new parcel will be inserted (separate from the original polygon that we split).

Demo

OSM was used to extract test data. Specifically, in the run below, data for the Herastrau Park in Bucharest was used. In OSM, the park is represented as a series of polygons, and one of the polygons was partitioned in parcels, each measuring 8000 square meters.

On the left side you can see the actual partition. On the right side, the same partition is displayed, except we also have some of the objects used throughout development, including bounding boxes, and extreme boundary points.

And to visualize how the algorithm works, here's an animation of how it partitions the polygon, parcel by parcel:

Challenges

  • While working on this algorithm implementation, at one point, it was required to find the extreme boundary points and assign them cardinal orientation labels. A first approach to determining the cardinal direction of the boundary points was to order them by ST_Azimuth relative to the centroid of the polygon (as the circle center passed to ST_Azimuth). This is a pitfall because the clock-wise order doesn't guarantee in any way membership to the N,E,S,W edges of the bounding box. It's perfectly possible that two such extreme points belong to the north edge, and one belongs to the west edge, and one to the south edge.
  • As it turns out, splitting the polygon multiple times, and using ST_Union to glue the remaining parts back together will create MULTIPOLYGON structures and which ST_ExteriorRing is not compatible with. So a convex hull was used in order to get one single polygon and compute the bounding box for it. This problem will surface after repeatedly cutting the input polygon.
  • Due to lack of knowledge about QGIS, some logic for an SVG-based 1 visualization was written in order to check the position of the cuts and the shape of the polygon after the cuts were made. Although QGIS has virtual layer feature, that was quite hard to use.
  • Early on, there was a need to visualize some objects from OSM in order to assess if they are fit for test data. It was easier to create a VIEW of those objects in PostgreSQL and then visualize that in QGIS (it was a bit hard to deal with all the menus in QGIS)
  • It would've been nice to have some automatically generated test polygons for each position of the road relative to the polygon. And some polygons that would be thiner near the top and wider near the bottom. However, the polygons extracted from OSM were used instead.
  • To get the extreme points of the polygon, the first approach used was to intersect exterior ring of the polygon with the exterior ring of the polygon's envelope (envelope meaning bounding box). This proved to be complicated and would definitely add complexity if one of the polygon's edges were axis-parallel as the result of the intersection would've been a line, and not a point. So the second approach was much simpler, just decomposing the polygon into its vertices (using ST_DumpPoints) and then picking the northmost, eastmost, westmost and southmost vertices. While this works, it only works for POLYGON and MULTIPOLYGON but it's not expected to work for CURVEPOLYGON.

After-thoughts

While implementing the logic for each corner-cut, it became more obvious that the logic for one corner would be enough, and the entire drawing could be rotated, the cut made, and then the resulting pieces rotated back. This could have decreased the size of the implementation.

No strict constraints are made on the shape of the parcels. Most of them are square, but some diverge from that shape quite a lot.

Some polygons may have U-shapes and if the closest-point on the boundary is one of the two tips of the U shape, the parcel that gets cut could turn out to be one made up of disconnected pieces.

Another gap that was not covered by the initial spec and is not currently handled is the row-excess which can form at the end of a row (in the example it will be at the end of the row). This row-excess is another factor that can generate thin parcels because once formed, the next split will use this excess and add a very thin portion beneath it to form the new parcel. This is not very helpful, but it's a known issue.

Someone pointed out that there should also be access paths/roads but that wasn't part of the problem statement. Another aspect that was not taken into consideration was surface inclination, it's assumed that subdivided area is flat.

Conclusion

In this post we've seen how PostGIS can be leveraged to build a land subdivision algorithm. The implementation for this algorithm is available on Github under MIT license.

Footnotes:

1

Satoshi Koda has a great blog post about this, and that was very useful reading while writing the SVG visualization for this project. His new blog is hosted here.