March 01, 2023

Easily center a marker on a complex line in GeoDjango

GeoDjango comes with a native way of retrieving a simple line's geographic center, but will put us off in the presence of complex MultiLineStrings. Let's help ourselves to more advanced GeoDjango ORM tools for a flexible approximation of a marker placement.

Those naughty multiple LineStrings

While getting the center of a LineString in GeoDjango is as convenient as invoking the centroid property on a Geometry instance with pretty reliable results, a MultiLineString will have us hit this wall where the point falls outside of the Geometry boundaries. The following figures show a couple of ways a MultiLineString can be represented:

Fig. 1 - Single MultiLineString with centroid

Fig. 2 - Stacked LineStrings with centroid

The reason behind this problem is that the centroid is calculated from the convex hull of the geometry. While there may be cases in which we welcome this behaviour, most of the time it turns out to be an unwanted side-effect.

Embracing simplicity and sticking to PostGIS

While we could resort to a pythonic solution such as sci-py's K-Neighbour implementation, let's try and see what we can accomplish with GeoDjango only.

Our gut feeling can take us to get the geometry center through a coeffiecient-based solution. After some research in the PostGIS realm, we stumble upon the ST_LineInterpolatePoint function. As our research progresses, we find out its input has to be forced to a LineString type of geometry, therefore making ST_LineMerge() the perfect companion to the aforementioned function.

GeoDjango provides a nifty API to incorporate PostGIS geographic functions in its own ORM:

from django.contrib.gis.db import models

 
class STLineInterpolatePoint(models.functions.GeoFunc):
    """A function returning a point interpolated along a line at a fractional location."""

    function = "ST_LineInterpolatePoint" 
    output_field = models.PointField()


class STLineMerge(models.functions.GeoFunc):
    """A function returning a LineString or MultiLineString formed by joining together the line elements of a MultiLineString."""

    function = "ST_LineMerge"
    output_field = models.LineStringField()

We can then use the two functions on the queryset where the geometry field lives: 

GeoQuerySet.objects.annotate(
    geometry_center=STLineInterpolatePoint(STLineMerge("geometry"), 0.5)
)

This implementation does give interesting results with the line on figure 1:

Fig. 3 - Single MultiLineString with interpolated center

But what about figure 2?

django.db.utils.InternalError: line_interpolate_point: 1st arg isn't a line

Uh-oh... we're getting an error from the PostGIS backend! This is because figure 2 - remember? - is actually multiple stacked lines.

Let's not throw in the towel just yet and backtrack a little bit: we can still calculate the point on the line that's as close to the centroid as possible!

Luckily, PostGIS provides us with the Closest Point function, which, unsurprisingly, returns the 2-dimensional point on a Geometry that is closest to another Geometry, based on the shortest line between these two points. Let's implement it:

from django.contrib.gis.db import models

 
class STClosestPoint(models.functions.GeoFunc):
    """A function returning the 2-dimensional point on a geometry that is closest to another geometry."""

    function = "ST_ClosestPoint" 
    output_field = models.PointField()

The function can then be used in combination with the Centroid one that GeoDjango provides natively, as follows:

from django.contrib.gis.db.models.functions import Centroid


GeoQuerySet.objects.annotate(
    geometry_center=STClosestPoint("geometry", Centroid("geometry"))
)

This approximation, while being a little less accurate than the interpolation method, works with our multiple stacked lines too:

Fig. 4 - Single MultiLineString with point on geometry closest to centroid

Fig. 5 - Stacked LineStrings with point on geometry closest to centroid

As demonstrated, an operation such as centering a marker does not require anything more complex than one spatial query even for a harsh scenario. If you're curious as to how to refine this feature even further, you can look into implementations of the nearest neighbour algorithm. On the other hand, they may easily turn out cumbersome and slow, especially in comparison with simpler solutions as those outlined above.

Aftermath

This research led me to propose the ClosestPoint function for inclusion in the main Django project, with support to PostGIS e Spatialite. The contribution was ultimately accepted as part of the contrib.gis package.