Opened 8 years ago

Closed 3 years ago

#5581 closed defect (wontfix)

Rasterize with MERGE_ALG=ADD burns too many values

Reported by: Mike Taves Owned by: warmerdam
Priority: normal Milestone: closed_because_of_github_migration
Component: Algorithms Version: svn-trunk
Severity: normal Keywords:
Cc:

Description

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:

http://imgur.com/Lgnb1Pr.png

>>> 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.

http://imgur.com/40FYgdv.png

>>> 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).

http://imgur.com/MHA3JpX.png

>>> 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.

Attachments (1)

polygonize_unit_tests.patch (7.9 KB ) - added by Mike Taves 8 years ago.
Polygonize unit tests for #5581

Download all attachments as: .zip

Change History (3)

by Mike Taves, 8 years ago

Attachment: polygonize_unit_tests.patch added

Polygonize unit tests for #5581

comment:1 by Mike Taves, 8 years ago

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.

comment:2 by Even Rouault, 3 years ago

Milestone: closed_because_of_github_migration
Resolution: wontfix
Status: newclosed

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.

Note: See TracTickets for help on using tickets.