<div dir="ltr">That's it. I was not closing dst_ds, I was closing only newArray.<div><br></div><div>Thank you Even. You are wonderful.</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Sat, Jun 28, 2014 at 3:14 PM, 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:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Le samedi 28 juin 2014 22:05:54, devendra dahal a écrit :<br>
<div class="">> Thank you David and Even.<br>
><br>
> I really apreciate your suggestions.<br>
><br>
> Even's suggestion of including quotation did help running the script<br>
> without problem. As you suggested I replaced the os.system with<br>
> subprocess.call. But still the shapefile has just one polygon with zero<br>
> value and it does not have the coordinate.<br>
><br>
> the "wetlands.tif" looks like this with two values 1 = wetlands (white) and<br>
> 0= other lands (black)<br>
><br>
> [image: Inline image 1]<br>
><br>
> This is what "wetland_poly.shp" looks and without projection.<br>
><br>
> [image: Inline image 2]<br>
<br>
</div>Hum, you didn't show your code, but a possibility is that your raster file<br>
isn't properly closed before gdal_polygonize.py is run.<br>
<br>
Make sure that dst_ds = None is executed before gdal_polygonize.<br>
Check that wetlands.tif has got the projection (with gdalinfo for example),<br>
and try running outsize of your python script gdal_polygonize to determine<br>
where issues are.<br>
<div class="HOEnZb"><div class="h5"><br>
><br>
> On Sat, Jun 28, 2014 at 1:55 AM, Even Rouault<br>
> <<a href="mailto:even.rouault@mines-paris.org">even.rouault@mines-paris.org</a>><br>
><br>
> wrote:<br>
> > Le samedi 28 juin 2014 05:26:00, devendra dahal a écrit :<br>
> > > 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<br>
> ><br>
> > some<br>
> ><br>
> > > wetlands area changes over the year. I planned to subset the wetlands<br>
> ><br>
> > from<br>
> ><br>
> > > the dataset. Then using polygonize.py to convert the raster to polygon<br>
> ><br>
> > and<br>
> ><br>
> > > compute the area of each wetlands extent annually. The below is my<br>
> ><br>
> > concise<br>
> ><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.<br>
> > > Now using the line below I tried to convert the wetland raster to<br>
> > > 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>
> > Deven,<br>
> ><br>
> > The issue here is the lack of quoting around ESRI Shapefile once it is<br>
> > ingested<br>
> > into ConvertPoly.<br>
> ><br>
> > You likely need something like<br>
> ><br>
> > ConvertPoly= "gdal_polygonize.py %s -f %s %s %s Value" %(OutFile, "\"ESRI<br>
> > Shapefile\"", PolyName + ".shp", PolyName)<br>
> ><br>
> > Using the Python subprocess module might also be less error prone where<br>
> > you 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<br>
> > don't<br>
> > actually need to specify -f "ESRI Shapefile" at all.<br>
> ><br>
> > Best regards,<br>
> ><br>
> > Even<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<br>
> ><br>
> > as<br>
> ><br>
> > > I expected from the wetlands.tif but the polygon file lost the<br>
> > > coordinate system. if I use my_raster.tif in the code above it runs<br>
> > > successfully but runs for more than a hour so obviously I don't want<br>
> > > 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>
> > --<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>
<br>
--<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>
</div></div></blockquote></div><br></div>