Rasterize with MERGE_ALG=ADD burns too many values

[Mike Taves]

Using the option 'MERGE_ALG=ADD' with RasterizeLayer may yield results with too many values burned. Furthermore, in combination with other options also yield unexpected results.

For the following examples, I'm using GDAL 1.11.0 from Python, with some helper calls. The burned values are 10, except when obtaining from Z coordinate.

import numpy as np
from osgeo import gdal, ogr

def ogr_data_source(geom_wkt_list):
    ds = ogr.GetDriverByName('Memory').CreateDataSource('')
    ly = ds.CreateLayer('')
    dfn = ly.GetLayerDefn()
    for wkt in geom_wkt_list:
        feat = ogr.Feature(dfn)
        feat.SetGeometryDirectly(ogr.Geometry(wkt=wkt))
        ly.CreateFeature(feat)
        feat = None
    return ds

def gdal_empty_raster_ds(shape, res=10.0):
    ny, nx = shape
    ds = gdal.GetDriverByName('MEM').Create('', nx, ny, 1, gdal.GDT_Int32)
    ds.SetGeoTransform([0.0, res, 0.0, res * ny, 0.0, -res])
    band = ds.GetRasterBand(1)
    band.Fill(0)
    return ds

def gdal_burn(rast_ds, feat_ds, options=[], value=10):
    err = gdal.RasterizeLayer(rast_ds, [1], feat_ds.GetLayer(), None, None, [value], options)
    assert err == 0, err
    return rast_ds.ReadAsArray()

def run_option_configurations(wkts):
    feat_ds = ogr_data_source(wkts)
    ar_replace = gdal_burn(gdal_empty_raster_ds((4, 3)), feat_ds, ['MERGE_ALG=REPLACE'])
    print('REPLACE:\n%s' % ar_replace)
    ar_add = gdal_burn(gdal_empty_raster_ds((4, 3)), feat_ds, ['MERGE_ALG=ADD'])
    print('ADD:\n%s' % ar_add)
    ar_add_touched = gdal_burn(gdal_empty_raster_ds((4, 3)), feat_ds, ['MERGE_ALG=ADD', 'ALL_TOUCHED=TRUE'])
    print('ADD/ALL_TOUCHED:\n%s' % ar_add_touched)
    ar_add_touched = gdal_burn(gdal_empty_raster_ds((4, 3)), feat_ds, ['ALL_TOUCHED=TRUE'])
    print('ALL_TOUCHED:\n%s' % ar_add_touched)
    ar_replace = gdal_burn(gdal_empty_raster_ds((4, 3)), feat_ds, ['BURN_VALUE_FROM=Z', 'MERGE_ALG=REPLACE'], 0)
    print('Z/REPLACE:\n%s' % ar_replace)
    ar_add = gdal_burn(gdal_empty_raster_ds((4, 3)), feat_ds, ['BURN_VALUE_FROM=Z', 'MERGE_ALG=ADD'], 0)
    print('Z/ADD:\n%s' % ar_add)
    ar_add_touched = gdal_burn(gdal_empty_raster_ds((4, 3)), feat_ds, ['BURN_VALUE_FROM=Z', 'MERGE_ALG=ADD', 'ALL_TOUCHED=TRUE'], 0)
    print('Z/ADD/ALL_TOUCHED:\n%s' % ar_add_touched)
    ar_add_touched = gdal_burn(gdal_empty_raster_ds((4, 3)), feat_ds, ['BURN_VALUE_FROM=Z', 'ALL_TOUCHED=TRUE'], 0)
    print('Z/ALL_TOUCHED:\n%s' % ar_add_touched)

---

Point types work normally as expected, with nothing strange happening:

>>> run_option_configurations(['POINT Z(14 29 10)', 'POINT Z(16.3 24.5 100)', 'POINT Z(13 13 100)', 'POINT Z(26.5 15.6 10000)'])
REPLACE:
[[ 0  0  0]
 [ 0 10  0]
 [ 0 10 10]
 [ 0  0  0]]
ADD:
[[ 0  0  0]
 [ 0 20  0]
 [ 0 10 10]
 [ 0  0  0]]
ADD/ALL_TOUCHED:
[[ 0  0  0]
 [ 0 20  0]
 [ 0 10 10]
 [ 0  0  0]]
ALL_TOUCHED:
[[ 0  0  0]
 [ 0 10  0]
 [ 0 10 10]
 [ 0  0  0]]
Z/REPLACE:
[[    0     0     0]
 [    0   100     0]
 [    0   100 10000]
 [    0     0     0]]
Z/ADD:
[[    0     0     0]
 [    0   110     0]
 [    0   100 10000]
 [    0     0     0]]
Z/ADD/ALL_TOUCHED:
[[    0     0     0]
 [    0   110     0]
 [    0   100 10000]
 [    0     0     0]]
Z/ALL_TOUCHED:
[[    0     0     0]
 [    0   100     0]
 [    0   100 10000]
 [    0     0     0]]

---

There are several problems with one simple linestring, where with only using ADD uses each segment of the linestring, rather than treating this as one object. In combination with ALL_TOUCHED is less understood what is going on. And with ALL_TOUCHED/BURN_VALUE_FROM=Z much less is understood.

>>> run_option_configurations(['LINESTRING Z(6 6 1, 7 17 10, 13 22 100, 14 28 1000, 28 37 10000)'])
REPLACE:
[[ 0  0 10]
 [ 0 10  0]
 [10  0  0]
 [10  0  0]]
ADD:
[[ 0  0 10]
 [ 0 30  0]
 [20  0  0]
 [10  0  0]]
ADD/ALL_TOUCHED:
[[ 0 10 10]
 [ 0 40  0]
 [20 20  0]
 [10  0  0]]
ALL_TOUCHED:
[[ 0 10 10]
 [ 0 10  0]
 [10 10  0]
 [10  0  0]]
Z/REPLACE:
[[    0     0 10000]
 [    0   100     0]
 [   10     0     0]
 [    1     0     0]]
Z/ADD:
[[    0     0 10000]
 [    0  2100     0]
 [   20     0     0]
 [    1     0     0]]
Z/ADD/ALL_TOUCHED:
[[   0 3000 4857]
 [   0 5364    0]
 [  22  119    0]
 [   4    0    0]]
Z/ALL_TOUCHED:
[[   0 3000 4857]
 [   0   64    0]
 [  12   64    0]
 [   4    0    0]]

---

And with one polygon, the ADD option works as expected, however the combination of ADD/ALL_TOUCHED does not make sense. It looks like the Z value from the last coordinate is used to represent the Z value for the whole polygon, which I agree is the best guess, since PolygonZ shapes are ambiguously defined in 3D space (i.e. interpolating between LinearRingZ coordinates is non-sense).

>>> run_option_configurations(['POLYGON Z((12 29 3, 10.6 23.3 5, 1 17 7, 13 2 11, 26.1 11.5 13, 28.6 21.4 17, 19 22 19, 17 29 23, 12 29 29))'])
REPLACE:
[[ 0  0  0]
 [ 0 10  0]
 [10 10 10]
 [ 0 10  0]]
ADD:
[[ 0  0  0]
 [ 0 10  0]
 [10 10 10]
 [ 0 10  0]]
ADD/ALL_TOUCHED:
[[ 0  0  0]
 [10 60 20]
 [40 10 30]
 [10 30 20]]
ALL_TOUCHED:
[[ 0  0  0]
 [10 10 10]
 [10 10 10]
 [10 10 10]]
Z/REPLACE:
[[ 0  0  0]
 [ 0 29  0]
 [29 29 29]
 [ 0 29  0]]
Z/ADD:
[[ 0  0  0]
 [ 0 29  0]
 [29 29 29]
 [ 0 29  0]]
Z/ADD/ALL_TOUCHED:
[[  0   0   0]
 [ 29 174  58]
 [116  29  87]
 [ 29  87  58]]
Z/ALL_TOUCHED:
[[ 0  0  0]
 [29 29 29]
 [29 29 29]
 [29 29 29]]

---

See related #5560 feature request.

[Mike Taves]

Polygonize unit tests for #5581

[Mike Taves]

I made some modifications to the above examples in the form of a ​unit test patch, where two tests fail, for the linestring and polygon examples.

[Even Rouault]

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to ​https://github.com/OSGeo/gdal/issues if it is not already reported there.