GDAL – How to Polygonize NDVI Raster

gdalpolygonizepythonvectorization

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

Best Answer

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

Related Solutions

[GIS] GDAL polygonize lines

In the end I wrote the following script that solved my problem. The script converts raster pixels with a specified value to vector lines. For example the blue pixels (value = 0) are converted to vector lines. There is definitly room to improve the script, as you can see in the result image. The script can be found and edited here.

Raster Image enter image description here

Raster Image and Vector roads enter image description here

import ogr, gdal, osr, os
import numpy as np
import itertools
from math import sqrt,ceil

def pixelOffset2coord(rasterfn,xOffset,yOffset):
    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    coordX = originX+pixelWidth*xOffset
    coordY = originY+pixelHeight*yOffset
    return coordX, coordY

def raster2array(rasterfn):
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    array = band.ReadAsArray()
    return array

def array2shp(array,outSHPfn,rasterfn,pixelValue):

    # max distance between points
    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    pixelWidth = geotransform[1]
    maxDistance = ceil(sqrt(2*pixelWidth*pixelWidth))
    print maxDistance

    # array2dict
    count = 0
    roadList = np.where(array == pixelValue)
    multipoint = ogr.Geometry(ogr.wkbMultiLineString)
    pointDict = {}
    for indexY in roadList[0]:
        indexX = roadList[1][count]
        Xcoord, Ycoord = pixelOffset2coord(rasterfn,indexX,indexY)
        pointDict[count] = (Xcoord, Ycoord)
        count += 1

    # dict2wkbMultiLineString
    multiline = ogr.Geometry(ogr.wkbMultiLineString)
    for i in itertools.combinations(pointDict.values(), 2):
        point1 = ogr.Geometry(ogr.wkbPoint)
        point1.AddPoint(i[0][0],i[0][4])
        point2 = ogr.Geometry(ogr.wkbPoint)
        point2.AddPoint(i[1][0],i[1][5])

        distance = point1.Distance(point2)

        if distance < maxDistance:
            line = ogr.Geometry(ogr.wkbLineString)
            line.AddPoint(i[0][0],i[0][6])
            line.AddPoint(i[1][0],i[1][7])
            multiline.AddGeometry(line)

    # wkbMultiLineString2shp
    shpDriver = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(outSHPfn):
        shpDriver.DeleteDataSource(outSHPfn)
    outDataSource = shpDriver.CreateDataSource(outSHPfn)
    outLayer = outDataSource.CreateLayer(outSHPfn, geom_type=ogr.wkbMultiLineString )
    featureDefn = outLayer.GetLayerDefn()
    outFeature = ogr.Feature(featureDefn)
    outFeature.SetGeometry(multiline)
    outLayer.CreateFeature(outFeature)


def main(rasterfn,outSHPfn,pixelValue):
    array = raster2array(rasterfn)
    array2shp(array,outSHPfn,rasterfn,pixelValue)

if __name__ == "__main__":
    rasterfn = 'test.tif'
    outSHPfn = 'test.shp'
    pixelValue = 0
    main(rasterfn,outSHPfn,pixelValue)

Python – Using GDAL Polygonize

You are confusing the use of the gdal_polygonize command line utility with the python function gdal.Polygonize().

As you mentioned, you've managed to use the command line utility successfully; however, the Python function works differently and expects different arguments than those specified in the utility. The first argument should be a GDAL Band object, not a string, so this is why you get your error.

To get the Band object you need to open the input file using gdal.Open() and use the GetRasterBand() method to get your intended band. Additionally, you need to create an output layer in which the resulting polygons will be created.

The Python GDAL/OGR Cookbook has a good example on how to use this function. The required parameters are explained in a bit more detail here.

Alternative

Following on from your comment, if you would prefer to keep using the command line utility, one solution is to call it from within a Python script using the subprocess module i.e.

import subprocess

script = 'PATH_TO_GDAL_POLYGONIZE'
for in_file in glob.glob("dir/*.asc"):
    out_file = in_file[:-4] + ".json"
    subprocess.call(["python",script,in_file,'-f','GeoJSON',out_file])

This loops through the files and updates the input/output paths.

This way you get the result that you are use to getting from the utility, but the ability to loop through your files. If you do plan on going this way, the subprocess documention will be useful.

Best Answer

Related Solutions

[GIS] GDAL polygonize lines

Python – Using GDAL Polygonize

Related Question