Densify and smooth out categorical raster with diagonal neighbors

Question

I'm classifying vegetation types in QGIS (3.10). I have a pretty coarse raster that I want to smooth out. The standard approach seems to be to polygonize and then apply various smoothing algorithms (e.g. Chaiken or sliding_averaging in Grass' v.generalize). This works fine for rectangular raster polygons, but doesn't smooth out thin diagonal features well. For instance, if the raster is

then the output has 3 islands of A's and 4 of C's. Technically, they may be one polygon pinched where the A's (resp C's) meet at the diagonal, but they appear as islands, not truly connected. Smoothing keeps sharp angles at the diagonal touchpoints rather than "building a bridge".

I'm thinking of a solution which involves building these diagonal bridges in a densified raster (for computational efficiency) prior to polygonizing. The output would be something like

It's a 4x4 nearest-neighbor densification of the raster, with then 2:2 and 2:1:1 diagonal touchpoints identified (marked with outline) and smoothed out by changing the cells in bold. Note the categories need to be prioritized, since in the example above at the centre-right, A "beats out" C in terms of building a bridge. This densified raster can then be polygonized with enough breadcrumbs provided for vector smoothing algorithms (if desired) to further adjust the curves while keeping desired connectedness.

Before I muck around in python or build some fine-tuned combo of QGIS/Grass processing algorithms to implement this, am I reinventing the wheel? Or are there equivalent simpler solutions?

I'd be OK doing the diagonal smoothing and bridge-building after polygonization rather than in raster form, but solutions with e.g. + and - buffering the staircase-shaped polygons seem unnecessarily complicated, especially in then filling gaps and resolving overlaps.

Answer 1

I would approach this problem using r.neighbors. This module assigns each pixel a value depending on values in some window around the pixel. By using method=maximum you can achieve what you refer to as smoothing. You would first have to assign values to each vegetation type such that the higher values "win" over neighboring vegetation types. In you example you mention the "A" beats out "C". So if "A" has a hgiher value than "C" then r.neighbors would assign the value of A to all pixels where A appears anywhere in the window surrounding those pixels.

Although, this would also alter the whole border between A and C regions, expanding the A region by window size at the expense of the C region.

Answer 2

Building on @Micha's answer, I was able to get this to work within the QGIS GIU with a careful stacking of Grass r.neighbors processing runs. In the below, 1418r is my base classification raster, but (as per Micha's suggestion) reclassified so higher-numbered classes are the ones that should "win".

1. Upsample each cell to 4x4 by Exporting as with 4x the resolution (so 7.5m since it's a satellite-derived 30m to start). I called this 1418r_u.

2. Build a majority layer 1418ru_mo5s running r.neighbors method=mode, square-5 kernel, from 1418r_u. This rounds off corner pixels of the 4x4s.

3. Also build bridge layer 1418ru_ma3c running r.neighbors method=max, circular-3 kernel. This ends up being precisely what finds the dominant class in the corner pixel where there is a 2:2 diagonal.

4. Build a decision mask layer in 2 steps.

   • 1st step, build 1418ru_i using r.neighbors method=interspersion, 3-square kernel (from 1418r_u)
   • 2nd step, build 1418ru_ia using r.neighbors method=average, 3-circular kernel, 1418ru_i

   This averaging of the interspersion is a kluge to pick out and differentiate the corner subpixels in exactly the right diagonal meeting points. Without it, a couple of pixels with different replacement requirements come up with the same mask value.

5. Build a stitched layer 1418rus according to the decision mask, e.g. with the QGIS Raster calculator as follows (could doubtless use the GDAL one, or use masks in Grass instead)

    ("1418ru_ia@1"=38.8)*"1418ru_mo5s@1"
    +("1418ru_ia@1"=48.8)*"1418ru_mo5s@1"
    +("1418ru_ia@1"=49)*"1418r_u@1"
    +("1418ru_ia@1"=46.2)*"1418ru_ma3c@1"
    +("1418ru_ia@1"<38)*"1418r_u@1"
    +("1418ru_ia@1">50)*"1418r_u@1"

This magical incantation works since 38.8 is the inside corner of a 3:1 meeting place in the original raster, 48.8 (48.75 but gets rounded) is where a bridge needs to be build in a 2:1:1 diagonal meeting, 49 is at a 2:1:1 straight-line meeting (could use mo5s there if you wanted an incursion of the 2-up class), and 46.2 is where there is a 2:2 diagonal meeting exactly where the bridge needs to be built.

Then can use usual techniques to smooth and/or polygonize, the above has made the topology be the desired one.

If I were doing this again, I might either do all the analysis in Grass (to use mask layers in r.neighbors rather than the raster calculator piecing together), or might bite the bullet and program the subsitution algorithm directly in R using raster or Python and rasterio. But it can be done natively with processing in QGIS...