<div dir="ltr">Thank you David and Even.<div><br></div><div>I really apreciate your suggestions. <br><div><br></div><div>Even's suggestion of including quotation did help running the script without problem. As you suggested I replaced the os.system with subprocess.call. But still the shapefile has just one polygon with zero value and it does not have the coordinate.</div>
<div><br></div><div>the <span style="font-family:arial,sans-serif;font-size:13px">"wetlands.tif" looks like this with two values 1 = wetlands (white) and 0= other lands (black)</span></div><div><br></div><div><img src="cid:ii_146e405a4c0cb694" alt="Inline image 1" width="339" height="281"><br>
</div><div><br></div><div>This is what "<span style="font-family:arial,sans-serif;font-size:13px">wetland_poly.shp" looks and without projection.</span></div><div><br></div><div class="gmail_extra"><img src="cid:ii_146e40898ac683fc" alt="Inline image 2" width="349" height="283"><br>
<br><div class="gmail_quote">On Sat, Jun 28, 2014 at 1:55 AM, Even Rouault <span dir="ltr"><<a href="mailto:even.rouault@mines-paris.org" target="_blank">even.rouault@mines-paris.org</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
Le samedi 28 juin 2014 05:26:00, devendra dahal a écrit :<br>
<div><div class="h5">> Dear All;<br>
><br>
> I have some annual landcovder raster datasets (preciously 20) with 7<br>
> different types.<br>
> Since the datasets covers large area and I just want to see trends of some<br>
> wetlands area changes over the year. I planned to subset the wetlands from<br>
> the dataset. Then using polygonize.py to convert the raster to polygon and<br>
> compute the area of each wetlands extent annually. The below is my concise<br>
> code.<br>
><br>
><br>
> file1 = "my_raster.tif"<br>
> OutFile = "wetlands.tif"<br>
> rastDriver = gdal.GetDriverByName('Gtiff') # reading gdal driver<br>
> rastDriver.Register()<br>
><br>
> src_ds = gdal.Open(file1)<br>
><br>
> srcband = src_ds.GetRasterBand(1)<br>
> cols = src_ds.RasterXSize<br>
> rows = src_ds.RasterYSize<br>
><br>
> Array1 = srcband.ReadAsArray(0, 0, cols, rows)<br>
> newArray = numpy.where(Array1 != 5, 0,Array1)<br>
><br>
><br>
> dst_ds = rastDriver.Create(OutFile, cols, rows, bands, GDT_Int16)<br>
> dst_ds.SetGeoTransform(geot)<br>
> dst_ds.SetProjection(proj)<br>
> dst_ds.GetRasterBand(1).WriteArray(newArray)<br>
> newArray = None<br>
><br>
><br>
> From the code above I created wetlands only raster and saved a file. Now<br>
> using the line below I tried to convert the wetland raster to polygon.<br>
><br>
> PolyName = "wetland_poly"<br>
> ConvertPoly= "gdal_polygonize.py %s -f %s %s %s Value" %(OutFile, "ESRI<br>
> Shapefile", PolyName + ".shp", PolyName)<br>
<br>
</div></div>Deven,<br>
<br>
The issue here is the lack of quoting around ESRI Shapefile once it is ingested<br>
into ConvertPoly.<br>
<br>
You likely need something like<br>
<div class=""><br>
ConvertPoly= "gdal_polygonize.py %s -f %s %s %s Value" %(OutFile, "\"ESRI<br>
Shapefile\"", PolyName + ".shp", PolyName)<br>
<br>
</div>Using the Python subprocess module might also be less error prone where you<br>
can pass an array of arguments instead of a string<br>
<br>
Hint: shapefile is the default vector format of all OGR utilities. So you don't<br>
actually need to specify -f "ESRI Shapefile" at all.<br>
<br>
Best regards,<br>
<br>
Even<br>
<div class=""><div class="h5"><br>
><br>
> os.system(ConvertPoly)<br>
><br>
><br>
> But it throws the error below.<br>
><br>
> dst_ds = drv.CreateDataSource( dst_filename )<br>
> AttributeError: 'NoneType' object has no attribute 'CreateDataSource'<br>
><br>
> When I tested qgis raster to vector conversion tool, I get the polygons as<br>
> I expected from the wetlands.tif but the polygon file lost the coordinate<br>
> system. if I use my_raster.tif in the code above it runs successfully but<br>
> runs for more than a hour so obviously I don't want to do this.<br>
><br>
> Any elegant suggestions for alternative idea or help to fix the script<br>
> would be highly appreciated.<br>
><br>
> Thank you very much<br>
><br>
> Deven<br>
<br>
</div></div><span class=""><font color="#888888">--<br>
Geospatial professional services<br>
<a href="http://even.rouault.free.fr/services.html" target="_blank">http://even.rouault.free.fr/services.html</a><br>
</font></span></blockquote></div><br></div></div></div>