7

I want to import the name of streets from OSM shapefile layer to another line layer based on neighborhood distance. The distances between two layers (line layers) are variable between 0 to 10 meter. The segments are not matched.

I already tried the Select by attribute tool in qgis under vector>data-management tools, but the result is not complete, and setting the parameter is tricky. I used overlap and set the distance differently but didn't t work really. for example, if I set the buffer distance higher I expect to get more matching and therefore more attributes but it is not so. Is there any solution?

Images: 1- The OSM shapefile with the name of streets 2- The Line layer which I want to import the names from OSM layer 3- The result

The OSM shapefile with the name of streets

The Line layer which I want to import the names from OSM layer

The result

  • Have you tried turning your line layer into a polygon layer using a 10m buffer and then spatial join? – user3386170 Mar 16 '18 at 13:54
1

You can benefit from "Hausdorff distance". There is a script available at github. But, it needs to be updated.

  • Create new script (Processing Toolbox > Scripts > Tools > Create new script),
  • Copy-paste this script in Script Editor,
  • Save the script to script folder (.qgis2/processing/scripts/) as any name

(Updated Script)

#Definition of inputs and outputs
#==================================
##[my scripts]=group
##origin_layer=vector
##target_layer=vector
##Road_name_field_in_OSM=field target_layer
##interval=number 1.0
##hausdorff_distance_weight=number 1.0
##length_difference_weight=number 1.0
##output=output vector

#Algorithm body
#==================================
from qgis.core import *
from PyQt4.QtCore import *
from processing.tools.vector import VectorWriter
import processing 
from scipy.spatial.distance import cdist
import numpy as np
from math import sqrt

def densify(polyline, interval):
    # densify the polyline using the given interval
    output = []
    for i in xrange(len(polyline) - 1):
        p1 = polyline[i]
        p2 = polyline[i + 1]
        output.append(p1)
        # calculate necessary number of points between p1 and p2
        pointsNumber = sqrt(p1.sqrDist(p2)) / interval
        if pointsNumber > 1:
            multiplier = 1.0 / float(pointsNumber)
        else:
            multiplier = 1
        for j in xrange(int(pointsNumber)):
            delta = multiplier * (j + 1)
            x = p1.x() + delta * (p2.x() - p1.x())
            y = p1.y() + delta * (p2.y() - p1.y())
            output.append(QgsPoint(x, y))
            if j + 1 == pointsNumber:
                break
    output.append(polyline[len(polyline) - 1])
    return output

def calculateHausdorffDistance(geom1,geom2):
    # calculate Hausdorff distance between two polylines
    distances=[]
    # calculate distances between origin and target feature
    D = cdist(geom1,geom2,'euclidean')
    H1 = np.max(np.min(D, axis=1))
    H2 = np.max(np.min(D, axis=0))
    distances.append( max(H1,H2) )
    # repeat the calculation in reverse order
    D = cdist(geom2,geom1,'euclidean')
    H1 = np.max(np.min(D, axis=1))
    H2 = np.max(np.min(D, axis=0))
    distances.append( max(H1,H2) )

    hausdorff = max(distances)
    return hausdorff


origin_layer = processing.getObject(origin_layer)
target_layer = processing.getObject(target_layer)
target_id_column_index = target_layer.fieldNameIndex(Road_name_field_in_OSM)

target_spatial_index = QgsSpatialIndex()
target_features = processing.features(target_layer)

origin_fields = origin_layer.pendingFields().toList()
origin_fields.append( QgsField("ROAD_NAME", QVariant.String ))
origin_fields.append( QgsField("HAUSDORFF", QVariant.Double ))
origin_fields.append( QgsField("LEN_DIFF", QVariant.Double ))
writer = VectorWriter(output, None, origin_fields, origin_layer.dataProvider().geometryType(), origin_layer.crs() )

outFeat = QgsFeature()

# populate the spatial index
for feat in target_features: 
    target_spatial_index.insertFeature(feat)

origin_features = processing.features(origin_layer)
for origin_feature in origin_features:
    center = origin_feature.geometry().centroid().asPoint()
    nearest_ids = target_spatial_index.nearestNeighbor(center,10)

    best_fit_id = None
    min_weight = None

    origin_geom = densify(origin_feature.geometry().asPolyline(), interval)

    for id in nearest_ids:
        target_feature = target_layer.getFeatures(QgsFeatureRequest().setFilterFid(id)).next()
        target_geom = densify(target_feature.geometry().asPolyline(), interval)
        hausdorff = calculateHausdorffDistance(origin_geom,target_geom)
        length_difference = abs(origin_feature.geometry().length() - target_feature.geometry().length())
        weight = hausdorff * hausdorff_distance_weight + length_difference * length_difference_weight

        if min_weight == None or weight < min_weight:
            min_weight = weight
            best_hausdorff_distance = hausdorff
            best_fit_id = target_feature.attributes()[target_id_column_index]
            best_length_difference = length_difference

    outFeat.setGeometry( origin_feature.geometry() )
    atMap = origin_feature.attributes()
    atMap.append(str(best_fit_id))
    atMap.append(float(best_hausdorff_distance))
    atMap.append(float(best_length_difference))
    outFeat.setAttributes( atMap )
    writer.addFeature( outFeat )

del writer
  • Run with parameters like below.

enter image description here

Output layer will include nearest road/street names.

enter image description here

  • And then, join output and lines by id (or by any matching field)

This is not a 100% solution for all data structure. You can get incorrect results. (Rarely I hope)

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.