Converting raster pixels to polygons with GDAL python

Question

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:

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:

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:

Answer 1

Fpolygonize is not the right tool to do that.

From the documentation:

This utility creates vector polygons for all connected regions of pixels in the raster sharing a common pixel value.

The idea of any polygonization tool is to minimize the number of polygons used to describe the original raster, and to minimize the overall length.

This is usually accomplished by aggregating first (i.e. merging, which you instead want to avoid) the pixel with similar value in the band to be analysed, then by running the actual polygonization algorithm.

Using GDAL it is definitely possible to obtain what you need, but not with Fpolygonize.

A (not very smart, maybe, but feasible) approach would be to iterate over each pixel and create a corresponding quadrangle in the output shapefile.

Answer 2

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.