4

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

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.

|improve this answer|||||
  • Hi @RafDouglas, Thanks for your time. Based on GDAL documentation, there are two options for polygonize "FPolygonize" and "Polygonize". In first time for vectorize, I used "Polygonize" class and result was more different from my expected. I edited my question to show result of "Polygonize" class. – nickan Oct 31 '19 at 15:42
  • 1
    Another (not smart but feasible) approach is to create a unique "cell ID" field for the pixels using the lat and long coordinates then run "Fpolygonize" on the resulting raster – Kartograaf Oct 31 '19 at 17:33
  • 1
    Right, but that would loose the original band information for the pixel. One could however transfer values in a subsequent step. @Cory G.: your answer is a better option. – RafDouglas Oct 31 '19 at 18:25
1

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.

|improve this answer|||||
  • Hi @Cory G, Unfortunately your method is not simple method and it needs another processes to produce first image in question. Also, the points from this code are in 4 corners of raster pixels and not in center of each pixel. – nickan Oct 31 '19 at 19:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.