Simplification of forest polygons

[wipfli]

Forests are nicely mapped in OpenStreetMap and I would like to show them on my map basically at all zoom levels. The problem with this is just that forest polygons have too much structure: forests have often boundaries which remind of fractals, forests have holes and they are cut into pieces by roads.

Here is an example of how forests look like at zoom z9 in Switzerland:

This was produced with something like the following planetiler code:

// in processFeatures()

    if (sourceFeature.canBePolygon() && (
        sourceFeature.hasTag("landuse", "forest") ||
        sourceFeature.hasTag("natural", "wood")
    )) {
      features.polygon("wood")
        .setMinZoom(0);
    }

// in posProcessFeatures()

    if ("wood".equals(layer)) {
      try {
        return FeatureMerge.mergeOverlappingPolygons(items, 4);
      }
      catch (GeometryException e) {
        return null;
      }
    }

The problem with the above map is that from my point of view, the forest polygons have too much detail. I would like to simplify them but I don't know of a good way of doing this. I played a bit around with mergeNearbyPolygons() but could not get satisfying results.

Basically what I would like to achieve is that at low zoom levels I get an idea of what areas are mostly forest or mostly not forest, and as I zoom in, the forest polygons should get more and more details.

For example, I would like to simplify the above forest polygons to something like the hand-drawn blobs:

Is that possible at all?

[wipfli]

To me it feels almost like what would help here would be to blur the forests...

[phanecak-maptiler]

Maybe looking at this may help:

https://github.com/openmaptiles/planetiler-openmaptiles/pull/50/files#diff-f73b4a64267e736ad7b3d0bed5e524502c0d9b83b258bc4f8f62f6b6bf5837a4

It is adjustment of the forest merging in planetiler-openmaptiles for z8-z9. Merging is there, "blurring" not. But as first step, merging might be good.

[farfromrefug]

@wipfli what do you think about this https://data.maptiler.com/downloads/dataset/landcover ?
I like what they have done. Wish they would share the way they did it

[farfromrefug]

Indeed but we will never know how bings and google do it while with maptiler we might as they are part of openmaptiler

[msbarry]

I'm guessing that global landcover datasets like that (and what google/apple maps have at low zooms) are generated by vectorizing raster datasets. There are a few existing global raster landcover data datasets (here's a good survey: https://gisgeography.com/free-global-land-cover-land-use-data/) but they might also be inferring their own low zoom vector layers from raw global satellite imagery since you're basically just trying to approximate what satellite imagery would look like.

[acalcutt]

you guys had me thinking about this map I made with ESA Climate Change Initiative NetCDF. Each of the colors in this aligns to a different land type (though the key said they should have been color)
https://tiles.wifidb.net/data/landcover/#4/0.00000/-0.00275
https://www.techidiots.net/notes/mapbox-gl-tileserver-gl/esa-climate-change-initiative-landcover-to-raster-mbtiles

I was wondering if it would be possible to 'vectorize' something like that...maybe by color. or possibly right from the NETCDF?

[farfromrefug]

@msbarry i think i will try to ask maptiler guys to see if they are willing to share how they generated it

[wipfli]

https://worldcover2020.esa.int/downloader is a nice data set as well, it looks like this:

To me it feels like it should be possible to first smooth out the image with a zoom-level dependent gaussian filter, then apply a step function which keeps only values above a certain threashold, and then vectorize the resulting image...

[msbarry]

What kind of tweaks did you try with mergeNearbyPolygons? That's basically grouping polygons that are within minDist of one another, then expanding those polygons by buffer pixels, shrinking them back down by that same buffer, and filtering-out any polygons less than minArea.

Here's what I get if I change postProcess to:

  @Override
  public List<VectorTile.Feature> postProcess(int zoom, List<VectorTile.Feature> items) throws GeometryException {
    for (var item : items) {
      item.attrs().remove(TEMP_NUM_POINTS_ATTR);
    }
    return FeatureMerge.mergeNearbyPolygons(items,
      128, // minArea
      128, // minHoleArea
      1, // minDist
      1 // buffer/unbuffer
    );
  }

That merges any polygons with the same features within 1px of eachother and the large values for minArea/minHoleArea filter filter-out any small polygons/holes that are left over.

You could also try tweaking other parameters within FeatureMerge, like switching to BufferParameters.JOIN_ROUND then using a larger buffer/unbuffer value to smooth the polygons even more, although the corner-circles start to become very obvious...

JTS provides a lot of other vector geometry utilities that we could use as building blocks for operations like this (anything you can do in postgis since it uses geos which is ported from JTS). Let me know if you have any other ideas to try based on jts/postgis geom utils?

Another route that might make sense here would be to rasterize the geometries, do some raster operations, then go back to vector geometries?

[wipfli]

Thanks for the input, @msbarry, this is really interesting. What I am worried about at the moment is that features appear, disappear, and appear again as you zoom into the map.

For example, here: https://wipfli.github.io/planetiler-forest-simplification/#map=10.02/46.9921/7.9663

At z9, we have this arm:

But at z10, it is gone again:

[wipfli]

I wonder if it would be better to increase the size of forest polygons before merging by some characteristic length scale in meters rather than using a fixed pixel value. If I had to guess I would say that the holes in a forest in Switzerland are on a scale of 100 m to 1000 m. Would be interesting to look at spatial frequencies of land cover raster data...

[msbarry]

postProcess takes a zoom level, so you could compute inputs to mergeNearbyPolygons based on that to keep a consistent buffer/min dist/min size ? GeoUtils has metersToPixelAtEquator to help with something like that.