[GIS] Converting raster pixels to polygons with GDAL python

gdalgdal-polygonizepythonqgisraster-conversion

I want to convert raster pixels to vector polygons as below image that is output of "Processing Toolbox–> Vector Creation–> Raster pixels to polygons" in QGIS 3.8:
enter image description here

I tried to do this vectorize with below code:

from osgeo import gdal, ogr, osr
import sys
import os


#gdal.UseExceptions()
os.chdir("H:/aa")
fileName = "H:/aa/image.tif"
src_ds = gdal.Open(fileName)
srs = osr.SpatialReference()
srs.ImportFromWkt(src_ds.GetProjection())
if src_ds is None:
    print('Unable to open %s' % src_fileName)
    sys.exit(1)
srcband = src_ds.GetRasterBand(1)
print (srcband)
dst_layername = "H:/aa/aapoly"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(dst_layername + ".shp")
dst_layer = dst_ds.CreateLayer(dst_layername , srs = srs, geom_type=ogr.wkbMultiPolygon)
newField = ogr.FieldDefn('DN', ogr.OFTReal)
dst_layer.CreateField(newField)
gdal.FPolygonize(srcband, None, dst_layer, 0, [], callback=None)
dst_ds.SyncToDisk()
dst_ds=None

But output of this code is:
enter image description here

In second image, pixels with same value that are neighborhood were merged in vector output, but I don't want to merge them. First image is desirable but second image as output of the code is not desirable.
It is important to do it with GDAL code and not another tools such as in Convert Raster to vector – create polygons based on each pixel.

With another class, "Polygonize", its result is more different from what I expected (Based on @RafDouglas suggestion in answer section):

from osgeo import gdal, ogr, osr
import sys
import os

#gdal.UseExceptions()
os.chdir("H:/aa")
fileName = "H:/aa/image.tif"
src_ds = gdal.Open(fileName)
srs = osr.SpatialReference()
srs.ImportFromWkt(src_ds.GetProjection())
if src_ds is None:
    print('Unable to open %s' % src_fileName)
    sys.exit(1)
srcband = src_ds.GetRasterBand(1)
print (srcband)
dst_layername = "H:/aa/aapoly"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(dst_layername + ".shp")
dst_layer = dst_ds.CreateLayer(dst_layername , srs = srs, geom_type=ogr.wkbMultiPolygon)
newField = ogr.FieldDefn('DN', ogr.OFTReal)
dst_layer.CreateField(newField)
gdal.Polygonize(srcband, None, dst_layer, 0, [], callback=None)
dst_ds.SyncToDisk()
dst_ds=None

Unfavorable Result:

enter image description here

Best Answer

This code will make a point vector layer with points in the center of all of your pixels and having an 'ID' field that contains a unique value for each. If you take the output layer and rasterize it using the bounding area and scale from the input (i.e. with gdal.Translate) you should have the raster you want.

import ogr, gdal, osr, os
import numpy as np
import itertools

def pixelOffset2coord(raster, xOffset,yOffset):
    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):

    # max distance between points
    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    pixelWidth = geotransform[1]

    # wkbPoint
    shpDriver = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(outSHPfn):
        shpDriver.DeleteDataSource(outSHPfn)
    outDataSource = shpDriver.CreateDataSource(outSHPfn)
    outLayer = outDataSource.CreateLayer(outSHPfn, geom_type=ogr.wkbPoint )
    featureDefn = outLayer.GetLayerDefn()
    outLayer.CreateField(ogr.FieldDefn("ID", ogr.OFTInteger))

    # array2dict
    point = ogr.Geometry(ogr.wkbPoint)
    row_count = array.shape[0]
    for ridx, row in enumerate(array):
        if ridx % 100 == 0:
            print("{0} of {1} rows processed" .format(ridx, row_count))
        for cidx, value in enumerate(row):
            Xcoord, Ycoord = pixelOffset2coord(raster,cidx,ridx)
            point.AddPoint(Xcoord, Ycoord)
            outFeature = ogr.Feature(featureDefn)
            outFeature.SetGeometry(point)
            outLayer.CreateFeature(outFeature)
            outFeature.SetField("ID", outFeature.GetFID())
            outLayer.SetFeature(outFeature)
            outFeature.Destroy()
    outDataSource.Destroy()

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

if __name__ == "__main__":
    rasterfn = r'C:/image.tif'
    outSHPfn = r'C:/points.shp'
    main(rasterfn,outSHPfn)

This will only work on small rasters due to the size constraints of ESRI's .shp file format.

Related Question