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
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:
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.
This will only work on small rasters due to the size constraints of ESRI's .shp file format.