[gdal-dev] having problem with polygonize

[Even Rouault]

Le samedi 28 juin 2014 22:05:54, devendra dahal a écrit :
> 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]

Hum, you didn't show your code, but a possibility is that your raster file 
isn't properly closed before gdal_polygonize.py is run.

Make sure that dst_ds = None is executed before gdal_polygonize.
Check that wetlands.tif has got the projection (with gdalinfo for example), 
and try running outsize of your python script gdal_polygonize to determine 
where issues are.

> 
> 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

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