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)
polygon.dbf
Description: polygon.dbf
polygon.prj
Description: polygon.prj
polygon.shp
Description: polygon.shp
polygon.shx
Description: polygon.shx
test.py
Description: test.py
_______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
