[gdal-dev] having problem with polygonize

[devendra dahal]

That's it. I was not closing dst_ds, I  was closing only newArray.

Thank you Even. You are wonderful.


On Sat, Jun 28, 2014 at 3:14 PM, Even Rouault <even.rouault at mines-paris.org>
wrote:

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