[gdal-dev] Problem with polygon and gdal.RasterizeLayer causing endless cpu usage

[Casper Børgesen]

Hi devs

I have stumpled upon an issue where I need some help to investigate further. I have a simple polygon which I would like to rasterise into a small quadratic raster (see the code below and also attached ready to run).

I think the problem is the geometry. Using FME and QGIS I have tried to validate the polygon file, but it came back valid. The attached file is created using "Save As..." in QGIS, but the issue also occurs on the original dataset. Running the script prints the WKT and its pretty simple with only 4 points.

Running the script below results in python getting stuck at the line with gdal.RasterizeLayer(...) with high CPU usage. No exceptions, no python crash and looking in task manager I see no change in memory consumption.

In the original file I had more polygons and the three neighbouring polygons (north, north east and east of the polygon) all have the same problems - but not the rest of the polygons in the dataset.

I have tried this using GDAL 2.2.0 x64 running on Windows 7.

Can anyone give me a hint to whats going on in this case? If this is a bug somehow, I will gladly report it. But I think it would help if I was able to narrow down the problem a bit more.

Regards, Casper


    # -*- coding: utf-8 -*-
    import os

    from osgeo import gdal, ogr

    src_polygons = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'polygon.shp')

    data_source = ogr.Open(src_polygons, 0)
    layer = data_source.GetLayer(0)
    projection = layer.GetSpatialRef()
    current_feature = layer.GetNextFeature()
    polygon = current_feature.GetGeometryRef().Clone()
    layer = None
    data_source = None

    print polygon.ExportToWkt()

    geotransform = [499000, 0.4, 0, 6095000, 0, -0.4]

    data_source = ogr.GetDriverByName('MEMORY').CreateDataSource('')
    layer = data_source.CreateLayer('', projection, geom_type=ogr.wkbPolygon)
    feature = ogr.Feature(layer.GetLayerDefn())
    feature.SetGeometry(polygon.Clone())
    layer.CreateFeature(feature)

    mask_ds = gdal.GetDriverByName('Mem').Create('', 5000, 5000, 1, gdal.GDT_Byte)
    mask_ds.SetGeoTransform(geotransform)
    mask_ds.SetProjection(projection.ExportToWkt())

    gdal.RasterizeLayer(mask_ds, [1], layer, burn_values=[1], options=["ALL_TOUCHED"])

    mask = mask_ds.ReadAsArray().astype(int)

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170602/328da37b/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: polygon.dbf
Type: application/octet-stream
Size: 78 bytes
Desc: polygon.dbf
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170602/328da37b/attachment-0005.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: polygon.prj
Type: application/octet-stream
Size: 388 bytes
Desc: polygon.prj
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170602/328da37b/attachment-0006.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: polygon.shp
Type: application/octet-stream
Size: 236 bytes
Desc: polygon.shp
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170602/328da37b/attachment-0007.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: polygon.shx
Type: application/octet-stream
Size: 108 bytes
Desc: polygon.shx
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170602/328da37b/attachment-0008.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test.py
Type: application/octet-stream
Size: 1044 bytes
Desc: test.py
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20170602/328da37b/attachment-0009.obj>