[gdal-dev] having problem with polygonize

devendra dahal kalu671 at gmail.com
Sat Jun 28 13:05:54 PDT 2014


Thank you David and Even.

I really apreciate your suggestions.

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.

the "wetlands.tif" looks like this with two values 1 = wetlands (white) and
0= other lands (black)

[image: Inline image 1]

This is what "wetland_poly.shp" looks and without projection.

[image: Inline image 2]

On Sat, Jun 28, 2014 at 1:55 AM, Even Rouault <even.rouault at mines-paris.org>
wrote:

> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140628/55cad394/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 13105 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140628/55cad394/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 2811 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140628/55cad394/attachment-0003.png>


More information about the gdal-dev mailing list