<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<meta name="Generator" content="Microsoft Word 14 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";
mso-fareast-language:EN-US;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
p.MsoAcetate, li.MsoAcetate, div.MsoAcetate
{mso-style-priority:99;
mso-style-link:"Balloon Text Char";
margin:0cm;
margin-bottom:.0001pt;
font-size:8.0pt;
font-family:"Tahoma","sans-serif";
mso-fareast-language:EN-US;}
span.EmailStyle17
{mso-style-type:personal-compose;
font-family:"Calibri","sans-serif";
color:windowtext;}
span.BalloonTextChar
{mso-style-name:"Balloon Text Char";
mso-style-priority:99;
mso-style-link:"Balloon Text";
font-family:"Tahoma","sans-serif";}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri","sans-serif";
mso-fareast-language:EN-US;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:3.0cm 2.0cm 3.0cm 2.0cm;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="DA" link="blue" vlink="purple">
<div class="WordSection1">
<p class="MsoNormal">Hi devs<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><span lang="EN-US">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).
<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US">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.<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US">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.<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US">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.<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US">I have tried this using GDAL 2.2.0 x64 running on Windows 7.<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US">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.<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US">Regards, Casper<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> # -*- coding: utf-8 -*-<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> import os<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> from osgeo import gdal, ogr<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> src_polygons = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'polygon.shp')<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> data_source = ogr.Open(src_polygons, 0)<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> layer = data_source.GetLayer(0)<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> projection = layer.GetSpatialRef()<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> current_feature = layer.GetNextFeature()<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> polygon = current_feature.GetGeometryRef().Clone()<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> layer = None<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> data_source = None<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> print polygon.ExportToWkt()<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> geotransform = [499000, 0.4, 0, 6095000, 0, -0.4]
<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> data_source = ogr.GetDriverByName('MEMORY').CreateDataSource('')<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> layer = data_source.CreateLayer('', projection, geom_type=ogr.wkbPolygon)<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> feature = ogr.Feature(layer.GetLayerDefn())
<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> feature.SetGeometry(polygon.Clone())<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> layer.CreateFeature(feature)<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> mask_ds = gdal.GetDriverByName('Mem').Create('', 5000, 5000, 1, gdal.GDT_Byte)<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> mask_ds.SetGeoTransform(geotransform)<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> mask_ds.SetProjection(projection.ExportToWkt())<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> gdal.RasterizeLayer(mask_ds, [1], layer, burn_values=[1], options=["ALL_TOUCHED"])<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""><o:p> </o:p></span></p>
<p class="MsoNormal"><span lang="EN-US" style="font-size:9.0pt;font-family:"Courier New""> mask = mask_ds.ReadAsArray().astype(int)<o:p></o:p></span></p>
<p class="MsoNormal"><span lang="EN-US"><o:p> </o:p></span></p>
</div>
</body>
</html>