[GIS] How to create a line layer which connects a point layer to a line layer

geoprocessingpyqgispythonqgisvector

In the map shown below, I have manually created the dummy lines (in green).
Using python, and without using PostGIS or sextante, I would like to automate the creation of these dummy links between a point and a line layer, based on the following criteria:

  • Depending on the polygon size, I should be able to add an X amount of dummy links between the centroid (point) and its nearest projection(s) onto the roadway layer.
  • The script so that the largest polygons get the most amount of dummy links and the smallest get the least amount.
  • The script should connect the newly created dummy link to nearest highway based on a classifications scheme, i.e. connect to "busier" highways first
  • If a polygon does not have any roads/links in it, the script should allow the centroid of that polygon to connect to roads outside of its boundaries within Y units of distance
  • The dummy links must be connected to the roadway layer in such a way that this roadway layer could be potentially used for routing including the dummy links.
  • Therefore, the dummy link layer must have the same attributes as the road layer
  • Create a new layer with all the new dummy links

Question:
Is this possible using Qgis and Python alone? any pointers as to where to start or can you help in providing a pseudocode to assist me to solve this problem?

enter image description here

Best Answer

I'd add my answer as a comment, but since I can't comment (no reputation enough), I'll add as an answer. It seems that you are trying to add centroid connectors to feed into a traffic assignment algorithm. Right? I should say that doing it automatically is not considered best practice in this industry, but it you are going to do it, I'd like to suggest you to make it proportional to the amount of network comprised within the polygon (roadway stock) and the displacement of this network as well.

I'll also mix some of the transportation modelling jargon with the GIS jargon...

First of all, are you considering splitting links or just connecting the new centroid connectors (dummy links) to the existing extremities of the links? Since I don't know what is the case, I'll assume that you want to connect them to the existing extremities. And it DOES make a huge difference in the algorithm if anybody is wondering.

Now is where I think there are two major questions you need to think about: 1st - How to weigh the criteria closest links versus busiest links 2nd - How to spread the connectors through the TAZ so you don't have them all clustered together?

I'll not answer thse questions, so when it comes to the part of the algorithm to actually add centroids, I'll add a single one.

The algorithm I am suggesting is:

  1. Compute the number of connectors to be added in each polygon (for each TAZ)
  2. Associate this number of connectors with the centroids file
  3. Extract all nodes from the road network (link extremities) to a new layer
  4. Tag TAZ (polygon) ID to each one of these nodes

  5. Loop through the centroids adding the connectors to a new layer

    Now let's put some codes to these steps (I hope the variables names are self explanatory):

Steps 1 & 2:

Let's assume there is TAZ field in both the centroids layer and the TAZ (polygon) layer. Something like this:

   connec={}

   features = taz_layer.getFeatures()
   idx=taz_layer.fieldNameIndex("TAZ")
   for f in features:
        taz_id=feat.attributes()[idx]
        geom = f.geometry()
        area= geom.area()
        connectors=area*number_of_connectors_per_unit_of_area
        connec[taz_id]=connectors

   features = centroids_layer.getFeatures()
   idx=centroids_layer.fieldNameIndex("TAZ")
   for f in features:
        taz_id=feat.attributes()[idx]
        a=centroids_layer.dataProvider().changeAttributeValues({f.id(): {idx:connec[taz_id]}})

   centroids_layer.updateFields()

Step 3:

    #Add node layer as suggested in the PyQGIS cookbook
    writer = QgsVectorFileWriter("node_layer.shp", "CP1250", fields, \
         QGis.WKBPoint, None, "ESRI Shapefile")

    idx=taz_layer.fieldNameIndex("TAZ")
    for feat in line_layer.getFeatures():
            link=list(feat.geometry().asPolyline())
            node_a=(link[0][0],link[0][1])
            node_b=(link[-1][0],link[-1][1])
            new_feat = QgsFeature()
            new_feat.setGeometry(QgsGeometry.fromPoint(QgsPoint(node_a)))
            new_feat.setAttributes(["how busy is the link", "name of the street", "etc", "etc])"#Here you need to capture and add all the characteristics of the links you want to consider when adding connectors.
            writer.addFeature(new_feat)

Step 4:

#This code is very slow, but will suffice #You should look into using spatial index, but this is what I could write by heart...

    idx=taz_layer.fieldNameIndex("TAZ")
    idx2=node_layer.fieldNameIndex("TAZ")
    for feat in node_layer.getFeatures(): #The extremities
        for feature in taz_layer.getFeatures():
             enclosed = feature.geometry().contains(feat.geom())
             if enclosed:
                 taz_id=feature.attributes()[idx]
                 a=node_layer.dataProvider().changeAttributeValues({feat.id(): {idx2:taz_id}})

    centroids_layer.updateFields()

Step 5:

     #Add link layer as suggested in the PyQGIS cookbook
    writer = QgsVectorFileWriter("node_layer.shp", "CP1250", fields, \
         QGis.WKBLinestring, None, "ESRI Shapefile")


    index = QgsSpatialIndex()
    for feature in node_layer.values():
       index.insertFeature(feature)

    for feature in centroids_layer.values():
        geom=pt.feature.geometry()
        pt = QgsPoint(geom.asPoint())
        nearest = index_from.nearestNeighbor(pt, 100)  #100 is a guess, you might need to do more...

        #Nearest contains the IDS of the nodes, so it is wise to build dictionaries or numpy arrays with all the relevant criteria (TAZ, if it is busy or not, etc) before you even start looping through centroids
        for i in nearest:  
           if criteria1(i) and criteria2(i) and ... and criterian(i):
              P=[pt, node_i.geometry().asPoint()]
              line_geom=QgsGeometry.fromPolyline(P)
              new_feat = QgsFeature()
              new_feat.setGeometry(line_geom)
              new_feat.setAttributes['Whatever','attributes','you','want']
              writer.addFeatures([new_feat])

Well... It ended up being an answer longer than I had anticipated, but it should get you started... And just as a warning, I did NOT test this code, so it is likely to have typos and other mistakes...

Related Question