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 <[email protected]> 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 >
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
