2

I have one huge problem and that is to get the value from one layer into another. It is a network with over a million links and I need to update a table from a layer of + 100 000 links. I am using QGIS and tried intersect, overlap on smaller data but I don't get the result I want. The two layers have identical geometry but one layer for example can have one long link compare to layer two that can have several links (road intersections) so overlap and intersection does not work. I did try the "Select by location" and "Include input features completely within the selection features" It did work as I wanted but when trying on the bigger data it just sit there for many hours (or overnight) and I don't know if it is working or just freezed in the process. All I want is to match the data and update the values from one layer to the other. As there is no identical ID I can not use the "Join" function but have to use geometry match instead. Is there another way around this problem?

  • What are you updating? Just attributes, just geometry, or both? Which is the target layer? The one with the longer geometry or the one with the intersected geometry? If the target is the longer geometry, what do you want to occur if two features that join to it have different attributes? Anyway, I would probably Spatial Join and keep only features that touch to limit the set of possible feature matches, then use a dictionary to store the join features and an update cursor to process the target features by Spatial Join IDs, compare the geometry, and make the updates where there is a match. – Richard Fairhurst Oct 24 '15 at 18:21
  • I am updating attributes only, target is the one with the intersected geometry so there will not be two features and both layers are shape files. – QGIS-user Oct 24 '15 at 19:55
  • I'd compute shorter dataset midpoints first. Spatial join points to long features, nearest. This might take long with 1 mil features. Points will get attributes of both sets. Count number of points (summarise) per long feature ID, remove ones with frequency=1. Remaining points = bridge between 2 sets – FelixIP Oct 24 '15 at 20:18
  • Any idea how to compute the midpoints for polylines? No python experience here for me and did not find any plugin neither... – QGIS-user Oct 24 '15 at 23:11
  • If you want to transfer attributes from one layer to another layer, when line segments are identical, my answer to a network related question may be helpful. Instead of adding attribut values you would do an assignment. Afterwards you have to dissolve the segments based on an unique attribute of the target layer. To customize the script some more details about your data structures are necessary. – Detlev Oct 26 '15 at 8:54
2

Given 6 layers the Python code executed in Python console (or executed from the editor) attributes will be transfered to features from layer 1. Layer 1 contains the features which provide the master geometry. They must be equal to, contained by, or disjoint from geometries of features in the other layers.

Layer 1 (differentiation of roads by number of lanes; attribute: Lanes; identified by feature id); in this example the master geometry

Layer 1: differentiation of roads by number of lanes

Layer 2 ( (differentiation by street name; attribute: street_name)

differentiation by street name

And so on.

Now the Python code:

from PyQt4.QtCore import QVariant
import datetime
print 'Starting at: ', datetime.datetime.now()

# get layers; one line for each layer to evaluate, with layer1 containing the master geometry
layer1 = QgsMapLayerRegistry.instance().mapLayersByName('thtn-vbg')[0]
layer2 = QgsMapLayerRegistry.instance().mapLayersByName('thtn-vbg_funk')[0]
layer3 = QgsMapLayerRegistry.instance().mapLayersByName('thtn-vbg_speed')[0]
layer4 = QgsMapLayerRegistry.instance().mapLayersByName('thtn-vbg_street_name')[0]
layer5 = QgsMapLayerRegistry.instance().mapLayersByName('thtn-vbg_one_way')[0]
layer6 = QgsMapLayerRegistry.instance().mapLayersByName('thtn-vbg_urban')[0]

# specify for each layer the attribute to transfer to the master geometry
attributes_dict = {layer2: 'KLASS', layer3: 'HTHAST', layer4: 'NAMN', layer5: 'RIKTNING', layer6: 'RLID'}

# split lines into segments and save all segments into a dict along with provided attributes
# master geometry; must be equal to, contained by, or disjoint from other geometries
segments = {}
for line in layer1.getFeatures():
    c = 1
    for i in range(len(line.geometry().asPolyline())-1):
        segments[((line.geometry().asPolyline()[i][0], line.geometry().asPolyline()[i][1]),
                        (line.geometry().asPolyline()[i+1][0], line.geometry().asPolyline()[i+1][1]))] = {'fid': line.id(), 'seq': c}
        c += 1

# other geometries with attributes
for the_layer, the_attrib in attributes_dict.items():
    for line in the_layer.getFeatures():
        for i in range(len(line.geometry().asPolyline())-1):
            key = ((line.geometry().asPolyline()[i][0], line.geometry().asPolyline()[i][1]),
                       (line.geometry().asPolyline()[i+1][0], line.geometry().asPolyline()[i+1][1]))
            if key in segments:
                segments[key][the_attrib] = line.attribute(the_attrib)

print 'Segments built: ', datetime.datetime.now()

# prepare segments for dissolving
composed = {}
for segment in segments:
    fid = segments[segment].get('fid', None)
    if fid:
        if fid not in composed:
            composed[fid] = {}
            _attributes = [segments[segment]['fid']]
            _attributes.extend([segments[segment].get(the_attrib, None) for the_layer, the_attrib in attributes_dict.items()])
            composed[fid]['attributes'] = _attributes
            composed[fid]['vertices'] = [(segments[segment]['seq'], (segment[0][0], segment[0][1]), (segment[1][0], segment[1][1]))]
        else:
            composed[fid]['vertices'].append((segments[segment]['seq'], (segment[0][0], segment[0][1]), (segment[1][0], segment[1][1])))

print 'Segments ordered: ', datetime.datetime.now()

# define QGIS layer with attributes
layer = QgsVectorLayer('LineString?crs=EPSG:3006', 'Result', 'memory')
prov = layer.dataProvider()
attributes = [QgsField('fid', QVariant.Int)]
attributes.extend([QgsField(the_attrib, QVariant.String) for the_layer, the_attrib in attributes_dict.items()])
attributes.append(QgsField('length', QVariant.Double))
prov.addAttributes(attributes)
QgsMapLayerRegistry.instance().addMapLayer(layer)

# create features
feats = []
layer.startEditing()

for compose in composed:
    composed[compose]['vertices'].sort()
    feat = QgsFeature()
    vertices = [QgsPoint(composed[compose]['vertices'][0][1][0], composed[compose]['vertices'][0][1][1]), 
    QgsPoint(composed[compose]['vertices'][0][2][0], composed[compose]['vertices'][0][2][1])]
    for idx in range(1, len(composed[compose]['vertices'])):
        vertices.append(QgsPoint(composed[compose]['vertices'][idx][2][0], composed[compose]['vertices'][idx][2][1]))
    feat.setGeometry(QgsGeometry.fromPolyline(vertices))
    composed[compose]['attributes'].append(feat.geometry().length())
    feat.setAttributes(composed[compose]['attributes'] )
    feats.append(feat)

print 'Lines built: ', datetime.datetime.now()

# finally add all created features and save edits
prov.addFeatures(feats)
layer.updateExtents()
layer.commitChanges()

print 'Finished: ', datetime.datetime.now()

The result is written to a memory layer, which you may save to disk. And this is the final result:

final result

Roads symbolized by number of lanes and labeled with streetname.

The script was tested with 16.000 master features, which receive attributes from 30.000 features in 5 layers. Running time was 1 minute on 3-years old PC with QGIS 2.8.2

  • Great! I will get back when I have tried. I will start with a smaller set first before getting into the big thing. Thanks – QGIS-user Oct 28 '15 at 20:42
  • Seems to have problem with python, I did install networkx module.. Python Console Use iface to access QGIS API interface or Type help(iface) for more info rom PyQt4.QtCore import QVariant File "<input>", line 1 rom PyQt4.QtCore import QVariant ^ SyntaxError: invalid syntax import networkx as nx Traceback (most recent call last): File "<input>", line 1, in <module> File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 478, in _import mod = _builtin_import(name, globals, locals, fromlist, level) ImportError: No module named networkx – QGIS-user Oct 30 '15 at 16:27
  • ImportError: No module named networkx - did your networkx installation went to the Python instance installed with QGIS? If so you should have a directory like this: `d:\OSGeo4W64\apps\Python27\Lib\site-packages\networkx`. I am not familiar with MacOS, don't know what path it has to be in your case. I suppose if the installation went somewhere else just copy the directory to the correct place. Not fine, but it works on another workstation of mine, where the same mistake happens. – Detlev Oct 30 '15 at 16:52
  • I did tried on PC, I have to check on MacOS later. It worked better with no error but the Result layer is empty it has "fid", "lanes" and "Type" (one I picked, attribute that has numbers only). It seems to run fine on this smaller set… but it gives empty tables – QGIS-user Oct 30 '15 at 17:06
  • Is it possible for you to share this small dataset, along with your code version, so I could check the process? – Detlev Oct 30 '15 at 17:30

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.