7

I have a raster file containing roads (left image). I polygonize the raster file with GDAL (see the script below). In the end I like to have vector lines. However GDAL only gives me something like this back (right image). Obviously this is correct, since polygonize creates polygons. Is there a way to "linize" a raster?

enter image description here

Here is my code:

import gdal,ogr,os

# open raster file
raster = gdal.Open('test.tif')
band = raster.GetRasterBand(1)

#create new shp file
newSHPfn = 'test.shp'
shpDriver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(newSHPfn):
    shpDriver.DeleteDataSource(newSHPfn)
outDataSource = shpDriver.CreateDataSource(newSHPfn)
outLayer = outDataSource.CreateLayer(newSHPfn, geom_type=ogr.wkbLineString )

# polygonize
gdal.Polygonize(band, None, outLayer, 1) 
  • Does 'gdal.ForceToMultiLineString' provide anything useful from your polygon? – rickhg12hs Nov 24 '13 at 16:06
  • @rickhg12hs I can't find to ForceToMultiLineString for GDAL. I only find it for OGR. Can you link to the documentation or explain a little more what you mean to do with it? – ustroetz Nov 25 '13 at 1:33
  • You're right, it does seem to be part of OGR. OGRGeometryH OGR_G_ForceToMultiLineString ( OGRGeometryH hGeom ). Is there a python interface for this? Will it "linize" your polygons? – rickhg12hs Nov 25 '13 at 21:10
  • 1
    Did you ever get anywhere with this? – oskarlin Oct 13 '14 at 11:18
  • 1
    @OskarKarlin Yes, I did. I posted my solution below. – ustroetz Oct 13 '14 at 11:35
4

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)
  • This works, but its very slow at for i in itertools.combinations(pointDict.values(), 2): , do you have any tips to make it faster? – sagarr Jan 28 '18 at 17:59
  • I would also be interested in a faster way to convert a raster to a vector line layer. Wouldn't it be an option to speed it up by creating single line geometries and iteratively writing them to a shapefile, instead of creating one large mulitline layer that is afterwards written once to a shapefile? – Sophie Crommelinck Apr 16 '18 at 10:02
  • Why don't you use the centroid of each polygon for the lines? I think it would give a much accurate result – ImanolUr Jul 4 '18 at 14:29

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.