[gdal-dev] having problem with polygonize

[Even Rouault]

Le samedi 28 juin 2014 05:26:00, devendra dahal a écrit :
> Dear All;
> 
> I have some annual landcovder raster datasets (preciously 20) with 7
> different types.
> Since the datasets covers large area and I just want to see trends of some
> wetlands area changes over the year. I planned to subset the wetlands from
> the dataset. Then using polygonize.py to convert the raster to polygon and
> compute the area of each wetlands extent annually. The below is my concise
> code.
> 
> 
> file1 = "my_raster.tif"
> OutFile = "wetlands.tif"
> rastDriver = gdal.GetDriverByName('Gtiff') # reading gdal driver
> rastDriver.Register()
> 
> src_ds = gdal.Open(file1)
> 
> srcband = src_ds.GetRasterBand(1)
> cols = src_ds.RasterXSize
> rows = src_ds.RasterYSize
> 
> Array1 = srcband.ReadAsArray(0, 0, cols, rows)
> newArray = numpy.where(Array1 != 5, 0,Array1)
> 
> 
> dst_ds = rastDriver.Create(OutFile, cols, rows, bands, GDT_Int16)
> dst_ds.SetGeoTransform(geot)
> dst_ds.SetProjection(proj)
> dst_ds.GetRasterBand(1).WriteArray(newArray)
> newArray = None
> 
> 
> From the code above I created wetlands only raster and saved a file. Now
> using the line below I tried to convert the wetland raster to polygon.
> 
> PolyName = "wetland_poly"
> ConvertPoly= "gdal_polygonize.py %s -f %s %s %s Value" %(OutFile, "ESRI
> Shapefile", PolyName + ".shp", PolyName)

Deven,

The issue here is the lack of quoting around ESRI Shapefile once it is ingested 
into ConvertPoly.

You likely need something like

ConvertPoly= "gdal_polygonize.py %s -f %s %s %s Value" %(OutFile, "\"ESRI 
Shapefile\"", PolyName + ".shp", PolyName)

Using the Python subprocess module might also be less error prone where you 
can pass an array of arguments instead of a string

Hint: shapefile is the default vector format of all OGR utilities. So you don't 
actually need to specify -f "ESRI Shapefile" at all.

Best regards,

Even

> 
> os.system(ConvertPoly)
> 
> 
> But it throws the error below.
> 
> dst_ds = drv.CreateDataSource( dst_filename )
> AttributeError: 'NoneType' object has no attribute 'CreateDataSource'
> 
> When I tested qgis raster to vector conversion tool, I get the polygons as
> I expected from the wetlands.tif but the polygon file lost the coordinate
> system. if I use my_raster.tif in the code above it runs successfully but
> runs for more than a hour so obviously I don't want to do this.
> 
> Any elegant suggestions for alternative idea or help to fix the script
> would be highly appreciated.
> 
> Thank you very much
> 
> Deven

-- 
Geospatial professional services
http://even.rouault.free.fr/services.html