isoband generation

[dsogari]

This ticket is intended to be a knowledge-sharing/improvement-seeking attempt to implement the isoband generation functionality in GDAL. My initial approach is simple. It accepts the following as input:

• the raster band handle (allowed to be of any integer or real data type, but not complex)
• the contour layer handle (assumed to be the output of GDALContourGenerate)
• a boundary layer handle (assumed to have one polygon feature, or it can be NULL)
• the output layer handle (for the isoband polygons, assumed to be initially empty)
• the kinds of information to gather for each isoband polygon (mean, max, min, etc.)
• a flag to inidicate that all pixels touching an isoband polygon will be considered when computing the statistics for that isoband
• a list of options (for instance, alternative names for the output fields can be specified)
• the progress function pointer
• the arguments for the progress function

The boundary can have any number of holes. If no boundary is provided, a default rectangular one will be created based on the raster extent.

Basically, the algorithm works like this:

1. Union contour lines
2. Split the boundary with the unioned conotur lines
3. Rasterize the resulting isobands in a memory dataset, to serve as key map*
4. Use this map to perform zonal statistic computation*
5. Write the computed info to appropriate fields*

*steps 3-5 are only needed if the caller actually requested computation of some statistic

The pros are:

1. it reuses the very-well tested GDALContourGenerate function
2. it provides additional functionality that can be easily tested
3. its performance depends mostly on the number/topology of contour lines and on the cost of the rasterization process, which is pretty fast

The downsides are:

1. the polygon splitting algorithm was adapted from QGIS, which is GPL-licensed
2. the rasterization has a memory cost equivalent to the size of the raster, since it is performed in-memory
3. this is not the optimal solution. A more theoretically correct way of implementing it would be to use the Marching-Squares algorithm for isobands

Some sample files for testing are attached. This code should be used in conjunction with the latest version of alg/contour.cpp (see #6519). Feel free to comment with ideas/fixes/improvements. Also, I tried to follow the code naming conventions used in GDAL as much as possible, but I may have missed something.

Anyway, I hope we can arrive at a good solution. Thanks

[Even Rouault]

I'm not familiar with the isoband functionality. A few remarks :

• regarding the part adapted from QGIS, no way to reimplement that without deriving from QGIS code ? We don't want GPL code in GDAL.
• Some algorithms offer the possibility to use a on-disk GTiff temp file, e.g. : alg/rasterfill.cpp

[dsogari]

An isoband is a filled contour band, delimited by the same limits used in the contour generation. Visually, the only difference is that the regions divided by the contour lines become actual polygons. These can be rendered in a visualization tool with arbitrary resolution, since they are vector data.

Regarding the part of polygon splitting, after adaptation to this purpose that code looks pretty different than the original, but still the idea remains the same. That's the reason why it is fast. It decomposes the boundary polygon into line strings, union that with all contour line strings, then polygonize it to obtain a list of candidate polygons. From these, only the ones that are actually part of the boundary (not in a hole) gets added to the output layer. I can't think of another way to do it.

As for the on-disk GTiff temp file, I'll look into that. Thanks

[vmo]

Sample code for closing contours in the contour algorithm

[vmo]

What I did in the contour.patch (no for merge yet) is to add an option to duplicate contours and give them the proper orientation and class, such that the contour merging part creates all the needed rings to create the polygons.

Except for the x2 in line storage, my guess it that it should perform reasonably well.

I'd like to try an avoid geos-polygonize because we have a lot of information in the current algorithm that we should be able to make use of.

[Even Rouault]

Implemented per ​https://github.com/OSGeo/gdal/pull/834