8

I want to polygonize a raster file. The raster file looks like this.

Here is the raster file

raster

Here is my codes in Python to polygonize the raster file. The codes source from Polygonize a Raster Band

#raster to polygons

from osgeo import gdal, ogr

#  get raster datasource


src_ds = gdal.Open( "C:\\Users\\select3.tif" )

srcband = src_ds.GetRasterBand(1)

#  create output datasource

dst_layername = "output_select3"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )
#The following two steps are very important, otherwise will generate error like (Coordinate not defined)
dst_ds.Destroy()
src_ds=None

I got the result which looks like this.

And here is the output shapefile Polygons

I do not satisfy with this, because I only want to have the polygon areas (Check mark) and exclude the extent (Cross mark) of the polygons. Could anybody tell me how to fix this in python?

6

This can be achieved by applying a mask as the second argument in the Polygonize function, as stated in the GDAL documentation. The mask needs to be a separate raster layer, that has 0 where you don't want the algorithm to polygonize.

With your data, follow these steps to implement:

1) Run raster calculator on your original raster ("select3.tif") with this expression to set the value of the areas you want to 1, and all the rest to 0:

enter image description here

2) Edit your code with these lines:

from osgeo import gdal, ogr

#  get raster datasource
src_ds = gdal.Open( "C:\\Users\\select3.tif" )
srcband = src_ds.GetRasterBand(1)

mask_ds = gdal.Open( "C:\\Users\\mask.tif" ) #Path to the mask layer generated above
maskband = mask_ds.GetRasterBand(1) 

#  create output datasource
dst_layername = "output_select3"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, maskband, dst_layer, -1, [], callback=None ) #Mask's band as second argument in Polygonize function

dst_ds.Destroy()
src_ds=None

Result:

enter image description here

5

Normally you don't need a mask for that task, but RGB values of the polygon areas equals 0. So you need to change them into a value (for example 125) and to change others into 0. Since all bands have same value, so it's sufficient to change one band to polygonize the raster as you did.

from osgeo import gdal, ogr

src_ds = gdal.Open( "C:\\Users\\select3.tif" )
srcband = src_ds.GetRasterBand(1)

# Invert pixel value to use srcband as a mask.
numpy_band = srcband.ReadAsArray() ##
numpy_band[numpy_band<125] = 125   ##
numpy_band[numpy_band>125] = 0     ##
srcband.WriteArray(numpy_band)     ##

dst_layername = "C:\\Users\\mask.tif"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, srcband, dst_layer, -1, [], callback=None )  ##
dst_ds.Destroy()
src_ds=None

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.