🌿GRASS is everything❤️‍🔥 (or using GRASS GIS to build elevation models in cities)

In the early days of OpenDroneMap, Dakota Benjamin and I discussed the possibility of integrating GRASS into OpenDroneMap to unlock more approaches for processing and post-processing data. We didn’t pull the trigger at that time, in part due to concerns about expanding the universe of dependencies, but later Piero Toffanin added GRASS to WebODM, and we haven’t looked back. Other libraries and tools at the heart of the ecosystem are OpenSfM, OpenMVS, MVSTexturing, PDAL, the venerable GDAL, and a range of other critical tools that make the OpenDroneMap ecosystem possible.

But let’s be honest, GRASS is everything. What other GIS product has William Shatner providing a fantastic introduction?

GRASS is an important part of the OpenDroneMap ecosystem, but not just for the functionality it provides, but also the quick testing of algorithms that it affords.

The Problem

Photogrammetrically derived point clouds don’t behave the same way as LiDAR datasets. Usually, when we want to filter these for e.g. finding ground points, we nevertheless use algorithms developed for lidar, such as Zhang, Bartels, or Pingel’s Progressive or Simple Morphological, Cloth Simulation, or Skewness Balancing Filters. These approaches work, but are underoptimized. I still need to test the Cloth Simulation Filter which holds promise for photogrammetric point clouds (ht Adam Steer), but in general, these approaches are very effective for LiDAR and less effective in photogrammetry. There are also some interesting, emerging, and sophisticated algorithms using deep learning floating around these days, including Gevaert, et al’s 2018 paper using morphological filters to train a deep learning model, a topic we should explore more deeply another time.

One of the limitations of algorithm’s developed for LiDAR is that they are often predicated on an assumption of a 3D dataset. Photogrammetric point clouds are rarely fully 3D in the same way that LiDAR is: by which I mean, we don’t capture data in the small gaps that LiDAR does, so we usually don’t have the layers of data we get from LiDAR. This limitation prompted the realization: maybe our approach to photogrammetric elevation model classification needs to be 2.5D, instead of 3D. (Side note, I’d love to calculate the actual dimensionality to photogrammetric point clouds, but that’s a whole other diversion).

Another advantage to using 2.5D info for classification, we can use our typical approaches for constructing a digital surface model as the starting place for further classification. This approach leverages filtering applied by OpenMVS based on the visibility calculations, some additional noise filtering in OpenDroneMap, but just uses the toolchain and the DSM products as-is.

Digital surface model over Zanzibar City, derived from Zanzibar Mapping Initiative data processed in OpenDroneMap 2.1.x.
Digital surface model over Zanzibar City, derived from Zanzibar Mapping Initiative data processed in OpenDroneMap 2.1.x.

A solution?

If you have followed my blog for a long time, you’ll know something of my fascination with landscape classification. I think it’s a category of analysis that is underutilized: there’s so much opportunity with the current generations of elevation models to measure, model, and understand the landscape using modern data and landscape classification: from predictions of habitat to improved soil mapping, to a host of problems we haven’t even realized we can solve with these tools and their derivatives. Combined with recent advances in machine learning approaches for automated classification of landscape, we have some real opportunities to advance our understanding and improve the management of our natural landscapes. But, to focus on the topic of the day, with the emergence of one such tool, GRASS GIS’ fantastic r.geomorphon tool, I became curious: can we use this for a different kind of “landform” analysis? It seems to be attempting to answer similar questions with respect to classification, specifically (and simplistically) what is the summit, ridge, shoulder, spur, slope, hollow, and footslope, which might be similar to buildings and tall vegetation. And what is a valley or depression, which may be similar to shape to ground. Geomorphon-based landform analysis, like the LiDAR morphological filters, is trying to answer similar landscape questions with respect to shape. However, deeper than that and intriguingly, it is an approach developed for a generic understanding of 2.5D cell elevation differential relationships. Thus we can potentially employ as much or as little of the ~500 geomorphons that are classifiable to develop an algorithm for ground classification of photogrammetrical derived point clouds.

Geomorphon calculation example using the EU DEM 25m (with search=11) from GRASS manual https://grass.osgeo.org/grass78/manuals/addons/r.geomorphon.html
Geomorphon calculation example using the EU DEM 25m (with search=11) from GRASS manual https://grass.osgeo.org/grass78/manuals/addons/r.geomorphon.html

A first pass

Geomorphons give us a lot of space to play in, from the calculation of them via parameterization, through the use of the raw, unclassified values in our own classification. It may be possible for us to create a stack of geomorphon raw calculations with different parameterizations and use the stack of the differently parameterized calculations in our own classification. For our initial analysis, we will use the classified forms as per the example above, and see how well the existing classes apply to our use case.

Of all the parameters available with geomorphons, I chose to only test the search distance, finding an optimum search distance around 100m for the test dataset. Below we compare the raw DSM data with the geomorphon classification.

If we filter this down to just valleys and depressions, we get a pretty promising result:

From here, we can create a mask of just the valleys and depressions and use that to mask our elevation data as follows:

The promise

This approach is barely optimized: I twisted one knob on the parameters and used default classification. Exploring the parameter space and coming up with a custom classification are on the docket for future exploration. I am fortunate enough with this dataset that I am in receipt of building footprint data that we can use for validation, but also leverage a combination of the human digitized portion with the ML classified version.

Looking at it this way, we can see that some valleys in the roofs are getting classified as valleys and propagating through to our elevation dataset. Further work in running geomorphons at multiple scales or diving into the unclassified raw data (or a combination) may help reduce this noise. We could also pretty safely set a size threshold for our areas of interest in order to remove these small inclusions. A great many of these are also areas where the algorithm is finding courtyards and similar that the building digitizers (understandably) ignored.

And now the cheat…

Since we have both the ML derived data and building footprint data, we can take advantage of the combination of learnings to both to create a good hydrology dataset. We could just do a simple interpolation of our filtered elevation model, and then use that for hydrological analysis.

The above approach is not uncommon. But aside from some edge cases, water doesn’t really flow through buildings, so more realistic modeling goes a bit deeper. First, since we haven’t optimized our geomorphon approach to do the best filtering we could, let’s leverage the building footprints to adugment our ML model:

Next we need to interpolate our remaining data. We could, as above, just use a basic inverse distance weighted interpolation to cover the entire area. This is a legitimate approach, but we might get better results by further constraining with an additional assumption: the direction of those alley ways should inform the weight and nature of interpolation. In other words: rather than use distance as our weight, use cost, and use the buildings as hard breaks.

Introducing v.surf.icw

GRASS has an add-on with a cost interpolator: v.surf.icw, in the following example being used for concentration interpolation in a bay:

File:Inlets 03 SurfSal icw big.png

How does this work with elevation data? The results are promising whether in our example in Zanzibar City:

or even the tall urban canyons of Stone Town:

Work remaining

We have already identified that the above approach would be more generalized by improvements to the ML approach, specifically by using multiple scales in combinations with a custom classification scheme. Additionally, for our use case where we have digitized building footprints, we might improve the fidelity of the approach by recombining the filtered data with the interpolated data so that we have the accuracy of the measured data in every place available, and interpolated data elsewhere. An advanced approach to recombine might weight the available measured and interpolated data based on a set of metrics of believability.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.