9

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
5

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
1

Thanks @ustroetz for your excellent answer, it really helped me a lot. I will share an adjustment and extra step I implemented to remove the unsightly loops that appear as a consequence of the process.

Note that this only works if your network is tree-like, which fortunately for my use case, mine always are.

I introduce geopandas and networkx for this process, as I find them a bit more wieldly than GDAL.

My implementation is rough and can certainly be handled more elegantly, but the theoretical process is:

  1. (adjustment) Write result as many lines instead of a single multiline
  2. Open the resultant file as a GeoDataFrame (geopandas)
  3. Convert the features into an undirected graph (networkx)
  4. Pick one of the terminal nodes and find all the edges involved in routing from it to all the other terminal nodes
  5. Create a subgraph containing only these edges
  6. Convert the subgraph into a GeoDataFrame
  7. Save the GeoDataFrame to another file

There are a couple of caveats, however:

  1. Unless you are more explicit than I, the first terminal node (routing source) is chosen arbitrarily. This can have implications at forks. You can see this effect in the second image, there is quite clearly a bias in the chosen paths because the source terminal node is located at the bottom of the network.
  2. Due to not compiling line features into a single multiline feature, a new line feature is generated for every "cell" from the source polygons. For large datasets, this could be an irritant.

Overall, though, going from raster to vector features is always going to be very messy business, so I'm satisfied with what I got from this.

import ogr, gdal, osr, os
import numpy as np
import itertools
from math import sqrt,ceil
import networkx as nx
from shapely.geometry import Point
import pandas as pd
import geopandas as gpd

def gdf_to_graph(gdf):

    # Find all unique start & end points and assign them an id
    gdf["start_node_coords"] = gdf["geometry"].apply(lambda x: x.coords[0])
    gdf["end_node_coords"] = gdf["geometry"].apply(lambda x: x.coords[-1])
    node_ids = {}
    i = 0
    for index, row in gdf.iterrows():
        node_1 = row["start_node_coords"]
        node_2 = row["end_node_coords"]
        if node_1 not in node_ids:
            node_ids[node_1] = i
            i += 1
        if node_2 not in node_ids:
            node_ids[node_2] = i
            i += 1

    # Assign the unique id to each
    gdf["source"] = gdf["start_node_coords"].apply(lambda x: node_ids[x])
    gdf["target"] = gdf["end_node_coords"].apply(lambda x: node_ids[x])

    gdf["length"] = gdf["geometry"].apply(lambda x: x.length)

    node_ids = {node_ids[k]: Point(k) for k in node_ids}
    node_df = pd.DataFrame.from_dict(node_ids, orient="index", columns=["geometry"])
    node_gdf = gpd.GeoDataFrame(node_df, geometry="geometry")
    # node_gdf.to_file("data/out/rivers_nodes.gpkg", driver="GPKG")

    # Drop these columns because they cannot be exported to GeoJSON
    gdf.drop("start_node_coords", axis=1, inplace=True)
    gdf.drop("end_node_coords", axis=1, inplace=True)

    # gdf.to_file("data/out/rivers_unbraided.gpkg", driver="GPKG")

    graph = nx.from_pandas_edgelist(gdf, edge_attr=["length", "geometry"])
    for n in node_gdf.index:
        graph.nodes[n]["geometry"] = node_gdf.iloc[n]

    return graph

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
    shpDriver = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(outSHPfn):
        shpDriver.DeleteDataSource(outSHPfn)
    outDataSource = shpDriver.CreateDataSource(outSHPfn)
    outLayer = outDataSource.CreateLayer(outSHPfn, geom_type=ogr.wkbLineString)
    featureDefn = outLayer.GetLayerDefn()
    outFeature = ogr.Feature(featureDefn)
    for i in itertools.combinations(pointDict.values(), 2):
        point1 = ogr.Geometry(ogr.wkbPoint)
        point1.AddPoint(i[0][0],i[0][1])
        point2 = ogr.Geometry(ogr.wkbPoint)
        point2.AddPoint(i[1][0],i[1][1])

        distance = point1.Distance(point2)

        if distance < maxDistance:
            line = ogr.Geometry(ogr.wkbLineString)
            line.AddPoint(i[0][0],i[0][1])
            line.AddPoint(i[1][0],i[1][1])

            line.AddGeometry(line)
            outFeature.SetGeometry(line)
            outLayer.CreateFeature(outFeature)


def main(rasterfn,outSHPfn,pixelValue):
    array = raster2array(rasterfn)
    array2shp(array,outSHPfn,rasterfn,pixelValue)


def path_to_edges(path, cycle=False):
    edges = []
    for i in range(len(path)-1):
        edges.append((path[i], path[i+1]))
    if cycle:
        edges.append((path[-1], path[0]))
    return edges


if __name__ == "__main__":
    rasterfn = 'data/pathfinder_output.tiff'
    outSHPfn = 'data/test.shp'
    pixelValue = 3
    main(rasterfn,outSHPfn,pixelValue)

    graph = gdf_to_graph(gpd.GeoDataFrame.from_file("data/test.shp"))
    terminal_nodes = [n for n in graph.nodes if nx.degree(graph, n) == 1]

    cycle_paths = [path_to_edges(path, cycle=True) for path in nx.cycle_basis(graph, terminal_nodes[0])]
    cycle_edges = set([item for sublist in cycle_paths for item in sublist])

    selected_edges = set()
    for terminal_node in terminal_nodes[1:]:
        edge_path = path_to_edges(nx.shortest_path(graph, terminal_nodes[0], terminal_node, weight="length"))
        selected_edges.update(set(edge_path))

    subgraph = graph.edge_subgraph(selected_edges)
    out_df = nx.to_pandas_edgelist(subgraph)
    out_gdf = gpd.GeoDataFrame(out_df)
    out_gdf.set_geometry("geometry")
    out_gdf.to_file("data/route_test.shp", driver="ESRI Shapefile")

Results:

enter image description here

enter image description here

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.