0

I have raster with NDVI values. I would like to polygonize the raster to get vector where each polygon is a pixel. I have tried gdal.Polygonize as described in the documentation and that runs without any error, but the result is just the number 0. I'm not sure if there is way I can display /read the result, or even save it. This is how I did it:

src_ds = gdal.Open( "imgs/response.tiff" )
if src_ds is None:
    print ('Unable to open')
    sys.exit(1)
    

srcband = src_ds.GetRasterBand(1)
    

    
    
dst_layername = "POLYGONIZED_STUFF"
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 )


>>>
0

I believe maybe there is problem with the GetRasterBand, but couldn't display it either.

My end goal is to be able to access the polygonized results

2 Answers 2

2

Firstly, I would question the wisdom of vectorizing a raster containing NDVI values. The vector dataset is likely to be huge and unwieldy.

Having said that, a fast solution would be to convert your raster into XYZ format (point coordinates) using GDAL translate:

gdal_translate -of XYZ NDVI.tif NDVI.xyz

If you really need polygons, you can do the following:

from osgeo import ogr, gdal
import os

# open image:
im = gdal.Open('filename.tif')
srcband = im.GetRasterBand(1)

if srcband.GetNoDataValue() is None:
    mask = None
else:
    mask = srcband

# create output vector:
outfile = 'POLYGONIZED_STUFF.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')

if os.path.exists(outfile):
    driver.DeleteDataSource(outfile)

vector = driver.CreateDataSource(outfile)
layer = vector.CreateLayer(outfile.replace('.shp', ''), im.GetSpatialRef(), geom_type=ogr.wkbPolygon)

# create field to write NDVI values:
field = ogr.FieldDefn('NDVI', ogr.OFTReal)
layer.CreateField(field)
del field

# polygonize:
gdal.Polygonize(srcband, mask, layer, 0, options=[], callback=gdal.TermProgress_nocb)

# close files:
del im, srcband, vector, layer
1

Slightly altering the accepted answer created polygons for negative pixel values for me:

mask = srcband.GetMaskBand()

I don't have enough reputation to just comment...

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.