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 <[email protected]> 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 > > <[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 > > -- > Geospatial professional services > http://even.rouault.free.fr/services.html >
_______________________________________________ gdal-dev mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/gdal-dev
