I want to polygonize .tiff raster file for specific values. For example, using one method I firstly reclassified the raster by gdal_calc.py --calc="A*logical_and(A>0,A<1)"
then used gdal.polygonize
function. But I know I can achieve this by providing a mask to the gdal.polygonize
function. How do I do that?
Update: I managed to do this sort of. I check the output GeoJSON with QGIS (raster to polygons -> Filter by value). The output has some issues: 1) it doesn't always compute, i.e. it gives GeoJSON of size 73 bytes, and 2) and it has more polygons than QGIS's result, which is odd.
from osgeo import gdal, ogr
import math
import os
sourceRaster = gdal.Open('ndvi.tiff')
band = sourceRaster.GetRasterBand(1) # Access the only band
outShapefile = "polygonized"
driver = ogr.GetDriverByName("GeoJSON")
if os.path.exists(outShapefile+".geojson"):
driver.DeleteDataSource(outShapefile+".geojson")
# We first create a new vector dataset and layer in it
outDatasource = driver.CreateDataSource(outShapefile+ ".geojson")
outLayer = outDatasource.CreateLayer("polygonized", srs=None)
# We add a new field ‘DN’ to the layer for storing the raster values for each of the polygons created
newField = ogr.FieldDefn('MYFLD', ogr.OFTInteger)
outLayer.CreateField(newField)
# Creating a copy datasource
driver2 = gdal.GetDriverByName( 'MEM' )
ds = driver2.CreateCopy('', sourceRaster)
band2 = ds.GetRasterBand(1)
needed_values = band2.ReadAsArray()
for row in needed_values:
for i, item in enumerate(row):
if (math.isnan(item) or not (0.0947260633111 < item < 0.09498295187950134)):
row[i] = -99
band2.WriteArray(needed_values)
band2.SetNoDataValue(-99)
for row in band2.ReadAsArray():
for item in row:
print(item)
# The fourth parameter is the index of the field to which the raster values shall be written (the index of the newly added ‘DN’ field in this case) (line 21)
gdal.Polygonize( band2, band2, outLayer, 0, [], callback=None )
band2.FlushCache()
ds = None
outDatasource.Destroy()
sourceRaster = None