GDAL polygonize líneas

Question

Tengo un archivo raster que contiene las carreteras (imagen de la izquierda). Yo polygonize el archivo raster con GDAL (ver la secuencia de comandos siguiente). En la final me gustaría tener líneas vectoriales. Sin embargo GDAL sólo me da algo como esto de nuevo (imagen de la derecha). Obviamente esto es correcto, ya que polygonize crea polígonos. Es allí una manera de "linize" una trama?

Aquí está mi código:

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) 

Answer 1

Al final me escribió la siguiente secuencia de comandos que resolver mi problema. El script convierte la trama de píxeles con un valor especificado de líneas vectoriales. Por ejemplo, el azul píxeles (valor = 0) se convierten en líneas vectoriales. Hay definitivamente habitación para mejorar la secuencia de comandos, como se puede ver en la imagen resultado. La secuencia de comandos se puede encontrar y editar aquí.

Imagen De Trama

Imagen Raster y Vector de carreteras

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)