Problem
Informacje o charakterystyce odcinka (np. jakości nawierzchni) zapisano w tabeli z identyfikatorem drogi i odległością do punktów początku i końca odcinka, liczoną od początku drogi – jej 'zerowego’ kilometra. Należało podzielić drogę na zadane odcinki, aby w prosty następnie zwizualizować np. kategorię stanu nawierzchni na danym odcinku.
Dane wejściowe
- Warstwa liniowa z drogami.
- Tabela z wartościami położenia pikiet, licząc od początku drogi. Każdy wiersz zawierał wartość odległości pikiety od początku drogi, identyfikator tej drogi i inne cechy odcinka.
Przygotowanie danych
- Połączenie odcinków danej drogi w w jeden obiekt (ponieważ w tabeli dostępne są dane o położeniu pikiet, licząc od początku drogi). Z uwagi na charakter danych (dwujezdniowość itp.) analiza i połączenie zostały wykonano ręcznie – za pomocą narzędzia Merge Selected Features. Powstała warstwa drogi.shp z atrybutem id_od zawierającym identyfikator drogi.
- Sprawdzenie prawidłowości kierunków linii , poprzez wizualizacją za pomocą markera z kierunkiem (na obrazku poniżej). Ta weryfikacja była koniczna, aby odległości do początków odcinków były mierzone zgodnie z kierunkiem liczenia kilometrażu drogi.
- Przeformatowanie i uproszczenie tabeli xlsx tabela, w której każdy wiersz odpowiadał jednemu odcinkowi, z przynajmniej dwoma wymaganymi atrybutami:
koniec_m – kilometraż od początku drogi do końca danego odcinak;
Droga – identyfikator drogi (zgodny z identyfikatorem w warstwie liniowej).
Oczywiście poza tym, były tam atrybuty opisujące cechy odcinka, np. jakość powierzchni.
Skrypt
Ponieważ wykonanie „ręczne” całego zadania, tylko z wykorzystaniem gotowych funkcji i narzędzi QGIS jest niemożliwe, stworzono skrypt Python.
from qgis.core import (QgsFeature, QgsGeometry, QgsPointXY, QgsProject,
QgsVectorLayer, QgsField, QgsFeatureRequest, QgsVectorDataProvider)
from PyQt5.QtCore import QVariant
import processing
# Ładuje warstwy --------------------------------------------------
drogi_layer = QgsProject.instance().mapLayersByName('drogi')[0]
tabela38 = QgsProject.instance().mapLayersByName('tabela')[0]
# -----------------------------------------------------------------
uri = "Point?crs=epsg:2180"
punkty_layer = QgsVectorLayer(uri, "Endpoint odcinka", "memory")
fields = tabela38.fields()
punkty_layer.dataProvider().addAttributes(fields.toList())
punkty_layer.updateFields()
punkty_layer.startEditing()
for rekord in tabela38.getFeatures():
id_drogi = rekord['Droga']
odleglosc = rekord['koniec_m']
query = f'"id_od" = \'{id_drogi}\''
segment_iterator = drogi_layer.getFeatures(QgsFeatureRequest().setFilterExpression(query))
segment = next(segment_iterator, None)
if segment:
punkt_na_drodze = segment.geometry().interpolate(odleglosc)
if punkt_na_drodze and not punkt_na_drodze.isMultipart():
feat = QgsFeature()
feat.setGeometry(punkt_na_drodze)
feat.setAttributes(rekord.attributes())
punkty_layer.addFeature(feat)
punkty_layer.commitChanges()
QgsProject.instance().addMapLayer(punkty_layer)
tolerance = 0.001
snapped_output = processing.run("qgis:snapgeometries", {
'INPUT': drogi_layer,
'REFERENCE_LAYER': punkty_layer,
'TOLERANCE': tolerance,
'BEHAVIOR': 0,
'OUTPUT': 'memory:'
})
snapped_layer = snapped_output['OUTPUT']
def find_closest_point_attributes(line_end, point_layer):
min_distance = float('inf')
closest_point_attributes = None
line_end_geom = QgsGeometry.fromPointXY(line_end)
for point_feat in point_layer.getFeatures():
point_geom = point_feat.geometry()
distance = line_end_geom.distance(point_geom)
if distance < min_distance:
min_distance = distance
closest_point_attributes = point_feat.attributes()
return closest_point_attributes
def divide_line_at_points(line_layer, point_layer):
divided_features = []
buffer_distance = 0.0001
buffers = [feat.geometry().buffer(buffer_distance, 5) for feat in point_layer.getFeatures()]
for line_feat in line_layer.getFeatures():
line_geom = line_feat.geometry()
for buffer_geom in buffers:
line_geom = line_geom.difference(buffer_geom)
if line_geom.isMultipart():
segments = line_geom.asGeometryCollection()
else:
segments = [line_geom]
for seg in segments:
if not seg.isEmpty():
seg_feat = QgsFeature()
seg_feat.setGeometry(seg)
# Get attributes from the closest point
if seg.isMultipart():
line_end = seg.asMultiPolyline()[-1][-1]
else:
line_end = seg.asPolyline()[-1]
seg_feat.setAttributes(find_closest_point_attributes(line_end, point_layer))
divided_features.append(seg_feat)
return divided_features
divided_lines = divide_line_at_points(snapped_layer, punkty_layer)
divided_uri = "LineString?crs=epsg:2180"
divided_layer = QgsVectorLayer(divided_uri, "Divided Lines", "memory")
divided_layer.dataProvider().addAttributes(punkty_layer.fields().toList())
divided_layer.updateFields()
divided_layer.startEditing()
for feat in divided_lines:
divided_layer.addFeature(feat)
divided_layer.commitChanges()
divided_layer.setName("Odcinki")
QgsProject.instance().addMapLayer(divided_layer)