[gdal-dev] having problem with polygonize

[devendra dahal]

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)

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140627/a3afa9d4/attachment.html>