3

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

1 Answer 1

0
needed_values = band2.ReadAsArray()

needed_values[np.isnan(needed_values)] = -99

needed_values[((needed_values[:, ] <= start_number) | (needed_values[:, ] >= stop_number))] = -99

band2.WriteArray(needed_values)
band2.SetNoDataValue(-99)
1
  • It can help you in your case, try it) Here numpy array is used
    – A.Rauan
    Commented Sep 3, 2019 at 11:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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