3

I am currently working on the problem of generating roads from Openstreetmap files. I use these Osm files to create a shapefile from all different road types with different roadway numbers, widths and so on. Later on, I would also like to use digital orthophotos (DOP) to eliminate inconsistencies and errors in the Osm files.

In the appendix you can see an example image of a DOP (with corrected shadows) and a generated mask, in addition to the mask there is the generated shapefile (with the different information contained).

So much for the introduction, now for the main question: In order to generate road markings on the generated roads it is first necessary to split shapefiles according to individual lanes. How is that possible in Python?

I want the generated polygons to be split up in such a way that they represent each lane and then map road markers to these polygons. So in my code a shapefile is generated, which buffers the reference lines of OSM depending on factors like number of tracks, track width etc.. (Further, the corresponding DOP is clipped in my code, so that only the relevant roads remain on the DOP. However, the result of this step will only be needed in the future and should become less relevant for this question.

Main Application

from qgis.core import *
from qgis.utils import *

from ui.main_window import MainWindow
from utils.files import getPathToTempFile
from utils.qgis_utils import loadVectorLayer, loadRasterLayer
from utils.osm import OsmBuildings, OsmStreets, OsmWater
from utils.vector import vector2raster, clipping, CopyExtentandEPSG


def main():
    app = QgsApplication([], True)
    app.setPrefixPath("C:/OSGeo4W64/apps/qgis", True)
    app.initQgis()
    window = MainWindow()
    window.show()
    window.raise_()
    bb = (635432, 5616219, 635939, 5616575) # WestSouthEastNorth

    # # streets
    streetsPath = ("files/cir/streets.shp")
    osmStreets = OsmStreets(extent=bb, extentEPSG=25832)
    osmStreets.getBufferedStreets(savePath=streetsPath)
    loadVectorLayer(streetsPath, "streets")
    streetMaskPath = ("files/cir/streets_mask.tiff")
    vector2raster(bb, streetsPath, streetMaskPath)
    loadRasterLayer(streetMaskPath, "streetsMask")

    # clipping dop with streetMask
    dop = ("files/colour/streets_dop.tif")
    clippedPathTemp = ("files/colour/streets_clipped_temp.tiff")
    clipping(streetMaskPath, dop, clippedPathTemp)

    # copy extent and EPSG from streetMaskPath
    clippedPath = ("files/colour/streets_clipped.tiff")
    CopyExtentandEPSG(streetMaskPath, clippedPathTemp, clippedPath)

    loadRasterLayer(clippedPath, "clippedSteets")

    # set CRS to the one of image and zoom to extent
    crs = QgsCoordinateReferenceSystem()
    crs.createFromId(25832)
    window.canvas.setDestinationCrs(crs)

    exitcode = app.exec_()
    QgsApplication.exitQgis()
    sys.exit(exitcode)

if __name__ == "__main__":
    main()

Needed functions (utils.osm)

class OsmStreets(Osm):
    def __init__(self, extent, extentEPSG=25832):
        Osm.__init__(self, extent, extentEPSG)

    def getBufferedStreets(self, savePath):
        """
        Loads streets form OSM and converts them to buffers (polygons).
        The size of the buffers will be calculated based on the OSM attributes (street width, sidewalk etc)

        :type savePath: str
        """
        filePathWithoutExtension, fileExtension = os.path.splitext(savePath)

        tempLinesPath = filePathWithoutExtension + "_temp_lines.geojson"
        self.loadFromOSM(tempLinesPath)

        self.createStreetBuffers(
            inputPath=tempLinesPath,
            outputPath=savePath)

        # cleanup
        # deleteFile(tempLinesPath)

    @staticmethod
    def createStreetBuffers(inputPath, outputPath):
        """
        Creates buffers of given streets (vector layer with lines).

        :param inputPath: path to OSM file (.geojson) with correct epsg
        :type inputPath: str
        :param outputPath: shape file that will be created
        :type outputPath: str
        """
        inputDataSource = ogr.Open(inputPath)  # type: ogr.DataSource
        inputLayer = inputDataSource.GetLayerByIndex(0)  # type: ogr.Layer

        # create output shape file
        outputDriver = ogr.GetDriverByName("ESRI Shapefile")
        deleteFile(outputPath)
        outputDataSource = outputDriver.CreateDataSource(outputPath)  # type: ogr.DataSource
        outputLayer = outputDataSource.CreateLayer(
            "streetBuffers", srs=inputLayer.GetSpatialRef(),
            geom_type=ogr.wkbPolygon
        )  # type: ogr.Layer

        # create fields of output layer

        # outputFieldName = 'turn:lanes'
        # fieldDefn = ogr.FieldDefn(outputFieldName, ogr.OFTString)
        # outputLayer.CreateField(fieldDefn)

        outputFieldName = 'highway'
        fieldDefn = ogr.FieldDefn(outputFieldName, ogr.OFTString)
        outputLayer.CreateField(fieldDefn)

        outputFieldName = 'lanes'
        fieldDefn = ogr.FieldDefn(outputFieldName, ogr.OFTString)
        outputLayer.CreateField(fieldDefn)

        outputFieldName = 'width'
        fieldDefn = ogr.FieldDefn(outputFieldName, ogr.OFTReal)
        fieldDefn.SetPrecision(2)
        fieldDefn.SetWidth(6)
        outputLayer.CreateField(fieldDefn)

        outputFieldName = 'sidewalk'
        fieldDefn = ogr.FieldDefn(outputFieldName, ogr.OFTString)
        outputLayer.CreateField(fieldDefn)

        outputFieldName = 'tunnel'
        fieldDefn = ogr.FieldDefn(outputFieldName, ogr.OFTString)
        outputLayer.CreateField(fieldDefn)

        outputFieldName = 'surface'
        fieldDefn = ogr.FieldDefn(outputFieldName, ogr.OFTString)
        outputLayer.CreateField(fieldDefn)

        outputFieldName = 'junction'
        fieldDefn = ogr.FieldDefn(outputFieldName, ogr.OFTString)
        outputLayer.CreateField(fieldDefn)


        featureDefn = outputLayer.GetLayerDefn()  # type: ogr.FeatureDefn
        for feature in inputLayer:  # type: ogr.Feature
            geom = feature.geometry()  # type: ogr.Geometry
            osmParameters = OsmStreetParameters(feature)
            streetWidth = osmParameters.getStreetWidth()
            bufferGeom = geom.Buffer(distance=streetWidth / 2.0)
            outputFeature = ogr.Feature(featureDefn)

            outputFeature.SetField("highway", osmParameters.highway)

            highway_select = outputFeature.GetFieldAsString("highway")
            # selection of "highways"
            if highway_select == "motorway" or \
               highway_select == "motorway_link" or \
               highway_select == "trunk" or \
               highway_select == "trunk_link" or \
               highway_select == "primary" or \
               highway_select == "primary_link" or \
               highway_select == "secondary" or \
               highway_select == "secondary_link" or \
               highway_select == "tertiary" or \
               highway_select == "tertiary_link" or \
               highway_select == "residential":

                outputFeature.SetField("width", streetWidth)
                outputFeature.SetField("sidewalk", osmParameters.sidewalk)
                outputFeature.SetField("lanes", osmParameters.lanes)
                outputFeature.SetField("tunnel", osmParameters.tunnel)
                outputFeature.SetField("surface", osmParameters.surface)
                outputFeature.SetField("junction", osmParameters.junction)

            # outputFeature.SetField("turn:lanes", osmParameters.turnlanes)
            else:
                continue

            outputFeature.SetGeometry(bufferGeom)
            outputLayer.CreateFeature(outputFeature)

        # invalid = outputFeature.GetFieldAsString("turn:lanes")
        # print invalid

        # save and cleanup
        outputDataSource.SyncToDisk()
        outputDataSource = None
        inputDataSource = None

    def loadFromOSM(self, savePath):
        """
        Loads streets ("highway") from OpenStreetMap and converts them to the epsg of the extent.

        :type savePath: str
        """
        minX, minY, maxX, maxY = self._transformExtent(self.extent, inputEPSG=self.extentEPSG, outputEPSG=4326)
        boundingBox = str((minY, minX, maxY, maxX))  # Note: y must be first here
        query = "way[highway]%s;" % boundingBox
        self.loadGeoJsonFromOverpassAndConvertToExtentEPSG(savePath, query)


class OsmStreetParameters:
    def __init__(self, feature):
        """
        Reads out attributes of given feature (feature will no be saved because reference are only temporally)

        :type feature: ogr.Feature
        """
        self.highway = feature.GetFieldAsString("highway")
        self.lanes = feature.GetFieldAsString("lanes")
        self.sidewalk = feature.GetFieldAsString("sidewalk")
        self.tunnel = feature.GetFieldAsString("tunnel")
        self.surface = feature.GetFieldAsString("surface")
        self.junction = feature.GetFieldAsString("junction")
        # self.turnlanes = feature.GetFieldAsString("turn" + ":" + "lanes")
        # print self.turnlanes

    def getStreetWidth(self):
        """
        :rtype: float
        """
        return self._getLaneWidth() * self._getLaneCount() + self._getSidewalkCount() * self._getSidewalkWidth()

    def _getLaneWidth(self):
        """
        Interpreted for germany from:
        http://wiki.openstreetmap.org/wiki/DE:Key:highway
        https://de.wikipedia.org/wiki/Richtlinien_f%C3%BCr_die_Anlage_von_Stra%C3%9Fen_%E2%80%93_Querschnitt
        """
        if self.highway == "motorway":
            return 3.5
        elif self.highway == "motorway_link":
            return 2.75
        elif self.highway == "trunk":
            return 3.5
        elif self.highway == "trunk_link":
            return 2.75
        elif self.highway == "primary":
            return 3.5
        elif self.highway == "primary_link":
            return 3.5
        elif self.highway == "secondary":
            return 3.0
        elif self.highway == "secondary_link":
            return 3.0
        elif self.highway == "tertiary":
            return 2.75
        elif self.highway == "tertiary_link":
            return 2.75
        elif self.highway == "unclassified":
            return 2.75
        elif self.highway == "residential":
            return 2.75
        elif self.highway == "service":
            return 2.75
        elif self.highway == "living_street":
            return 2.75
        else:
            return 2.0

    def _getLaneCount(self):
        """
        :rtype: int
        """
        if self.lanes is None or self.lanes == "":
            if self._getLaneWidth() == 2:
                return 1
            else:
                return 2
        else:
            try:
                return float(self.lanes)
            except ValueError:
                return 2

    def _getSidewalkCount(self):
        """
        http://wiki.openstreetmap.org/wiki/DE:Key:sidewalk
        :rtype: int
        """
        if self.sidewalk is None or self.sidewalk == "" or self.sidewalk == "none":
            return 0
        elif self.sidewalk == "both":
            return 2
        elif self.sidewalk == "left" or self.sidewalk == "right":  # todo: make difference between left and right sidewalks
            return 1
        else:
            return 0

    @staticmethod
    def _getSidewalkWidth():
        """
        Values based on https://de.wikipedia.org/wiki/Gehweg
        :rtype: float
        """
        return 2.5

enter image description here

  • I suspect thats actually quite a challenge: I can't recall any built-in solution in either pyshp or pyqgis, but I might be wrong here. there is the st_approximatemedialaxis in postgis that could help you make that work. ...you might be better off buffering the road lanes in the OSM data by its widths and try to apply your DOP enhancement on those somehow, if necessary – ThingumaBob Jan 8 '18 at 9:56
  • Well, the DOP were more like a description of where to go. The main focus will be on the generation of individual lanes, in order to map/generate lane markings on these lanes in further steps. With my algorithm I generate roads with different widths based on the number of tracks etc. I'm more interested in sharing the generated shapefiles, based on the number of lanes. – freddykrueger Jan 8 '18 at 10:05
  • ...well, then I got this wrong, I thought you'd like the generated polygons to be splitted to represent each lane. but, most of the time, actual splitting of shapefiles (if that is indeed what you try to accomplish) can (and probably should) be avoided by proper feature attributation, or be seen as a call for better data storage... – ThingumaBob Jan 8 '18 at 11:09
  • I think we were just passing each other by. I want the polygons created to be split to represent each lane. How is that possible? – freddykrueger Jan 8 '18 at 11:19
  • 3
    while I personally think it´s a bit harsh to put your question on hold (I get that you are looking for more general directions to where to start) you should definetely update your question with the info from your comments, correct the teminology (e.g. 'shapefiles' -> 'polygons'), share any code you have and/or be more specific and focussed on how to get what you want step by step (with some effort on your side in between and maybe new questions on isolated problems) without the need for others to write the whole script for you ,) – ThingumaBob Jan 8 '18 at 11:54
2

Do I get you right, that you want to clip the OSM road linetrings by your generated road polygons ('shapefile' in your illustration)? I assume, you are working from QGIS, so if QGIS processing module from PyQGIS is an option for you, you could use something like

import processing
clipped_roads = processing.runalg('qgis:clip', osm_road_linetrstrings, your_shapefile, None)
clipped_roads_layer = QgsVectorLayer(clipped_roads['OUTPUT'], 'Clipped OSM Roads', 'ogr')
QgsMapLayerRegistry.instance().addMapLayer(clipped_roads_layer, False)
  • The focus is on the generation of individual tracks, in order to map or generate the road markings on these tracks in further steps. With my algorithm I generate roads with different widths, based on the number of tracks, etc. My question refers to dividing the generated shapefiles based on the number of tracks, so that I can then map road markings to the individual tracks. – freddykrueger Jan 8 '18 at 10:41

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.