Dziś zajmiemy się demonstracją automatyzacji działań z poprzedniego wpisu pt.QGIS i tworzenie mapy do slużebności przesyłu.
Ściągnijmy najpierw paczkę danych przygotowaną w ramach poprzedniego wpisu. W załączeniu projekt z danymi (spakowany 7zip)
Otwieramy projekt:
Mamy tu 3 warstwy, najpierw w konsoli Pythona wyszukajmy te wartstwy i przepiszmy do zmiennych. Skoro mamy tylko 3 warstwy i nie musimy sie upewniać że nie istnieje więcej niż 1 warstwa o tej samej nazwie użyjemy prostego :
https://qgis.org/pyqgis/3.44/core/QgsProject.html#qgis.core.QgsProject.mapLayersByName
woda_lyr = QgsProject.instance().mapLayersByName('woda')[0]
kanal_lyr = QgsProject.instance().mapLayersByName('kanal')[0]
dzialki_lyr = QgsProject.instance().mapLayersByName('dzialki_temp')[0]
print(woda_lyr)
print(kanal_lyr)
print(dzialki_lyr)
output:
<QgsVectorLayer: 'woda' (ogr)>
<QgsVectorLayer: 'kanal' (ogr)>
<QgsVectorLayer: 'dzialki_temp' (ogr)>
Dlaczego [0] bo mapLayersByName zwraca listę warstw spełniających kryteria.
Stwórzmy także 2 warstwy tymczasowe wyjściowe wynikowy_poligon_lyr (poligonową) i wynikowy_linia_lyr (liniową). Po specyfikację odsyłam do https://qgis.org/pyqgis/3.44/core/QgsVectorLayer.html#module-QgsVectorLayer. My stworzymy 2 warstwy jedną (poligonową) z samym polem id, drugą (liniową) z polem dn (średnica). Na końcu pobieramy instancje projektu i dodajemy stworzone warstwy do listy warstw.
uri_poligon_lyr = "polygon?crs=epsg:2178&field=id:integer"
wynikowy_poligon_lyr = QgsVectorLayer(uri_poligon_lyr, "wynikowy_poligon", "memory")
uri_linia_lyr = "linestring?crs=epsg:2178&field=id:integer&field=dn:string"
wynikowy_linia_lyr = QgsVectorLayer(uri_linia_lyr, "wynikowy_linia", "memory")
QgsProject.instance().addMapLayer(wynikowy_poligon_lyr)
QgsProject.instance().addMapLayer(wynikowy_linia_lyr)
Zaznaczymy teraz ręcznie obiekty na warstwach, wszystkie odcinki kanałowe, wodociągowe i wybraną działkę.
Nasz skrypt będzie pobierał zaznaczone obiekty z danej warstwy iterował je i przetwarzał. Dodajmy jeszcze zmienną aggr_buf, która w przyszłości będzie przechowywała zagregowany bufor.
aggr_buf = None
Użyjemy do tego prostej funkcji selectedFeatures()tworzącej listę zaznaczonych obiektów. : https://qgis.org/pyqgis/3.44/core/QgsVectorLayer.html#qgis.core.QgsVectorLayer.selectedFeatures
for woda_lyr_feats in woda_lyr.selectedFeatures():
print(woda_lyr_feats['id'])
for kanal_lyr_feats in kanal_lyr.selectedFeatures():
print(kanal_lyr_feats['id'])
output:
1
2
3
4
5
1
2
3
4
5
Potem utworzymy bufor wokół każdego obiektu, plus zagregujemy wszystko w jednej geometrii. Ale najpierw na podstawie średnicy obliczmy szerokość bufora. Przydałaby się jakaś kontrola błędów. Dodajmy w pętli „for” zamiast linii print(…
try:
dn = int(woda_lyr_feats['DN'])
except:
dn = 0
zakres_m = 1+dn/2000
To wyeliminuje wartości nieprzetwarzalne na liczbę całkowitą takie jak NULL, tekstowe. Ogólnie wszystko czego nie znamy uznajemy ze ma średnice 0.Teraz utworzymy bufor wokół obiektu. Odsyłam do specyfikacji : https://qgis.org/pyqgis/3.44/core/QgsGeometry.html#qgis.core.QgsGeometry.buffer, i dodamy to do geometrii aggr_buf funkcją combine https://qgis.org/pyqgis/3.44/core/QgsGeometry.html#qgis.core.QgsGeometry.combine.
bufor_temp = woda_lyr_feats.geometry().buffer(zakres_m,5)
if aggr_buf is None:
aggr_buf = bufor_temp
else:
aggr_buf = aggr_buf.combine(bufor_temp)
tak samo robimy dla kanal_lyr
try:
dn = int(kanal_lyr_feats['DN'])
except:
dn = 0
zakres_m = 1+dn/2000
bufor_temp = kanal_lyr_feats.geometry().buffer(zakres_m,5)
if aggr_buf is None:
aggr_buf = bufor_temp
else:
aggr_buf = aggr_buf.combine(bufor_temp)
Teraz bierzemy wszystkie zaznaczone geometrie z dzialki_temp, iterując wycinamy geometrię działki, zagregowaną wcześniej geometrią. Patrz: https://qgis.org/pyqgis/3.44/core/QgsGeometry.html#qgis.core.QgsGeometry.intersection, i jeśli wyjdzie nam geometria niepusta to …
if aggr_buf is not None:
for dzialki_lyr_feats in dzialki_lyr.selectedFeatures():
aggr_buf_dzialka = aggr_buf.intersection(dzialki_lyr_feats.geometry())
if not aggr_buf_dzialka.isEmpty():
… dodajemy to do naszej tymczasowej warstwy poligonowej. Uwaga, nie dzielę poligonów, nie sprawdzam czy nie ma szczątkowych poligonów, a czasem nawet punkty mogą powstać. To wymaga osobnego potraktowania, ale nie chce gmatwać kodu. Features dodajemy przez dataProvider-a. Wyjaśnienie poniżej.
feature = QgsFeature(wynikowy_poligon_lyr.fields())
feature.setGeometry(aggr_buf_dzialka)
feature["id"] = wynikowy_poligon_lyr.featureCount() + 1
wynikowy_poligon_lyr.dataProvider().addFeatures([feature])
No i po uruchomieniu kodu mamy:
Pozostaje nam wyciąć wszystkie domiary. Narzędzia i działanie takie samo jak w przypadku buforów, wycinamy linie przy pomocy geometrii zaznaczonych działek, poniżej kod:
for kanal_lyr_feats in kanal_lyr.selectedFeatures():
domiar_wyciety = kanal_lyr_feats.geometry().intersection(dzialki_lyr_feats.geometry())
if not domiar_wyciety.isEmpty():
feature = QgsFeature(wynikowy_linia_lyr.fields())
feature.setGeometry(domiar_wyciety)
feature["id"] = wynikowy_linia_lyr.featureCount() + 1
try:
dn = int(kanal_lyr_feats['dn'])
except:
dn = 0
feature["dn"] = dn
wynikowy_linia_lyr.dataProvider().addFeatures([feature])
for woda_lyr_feats in woda_lyr.selectedFeatures():
domiar_wyciety = woda_lyr_feats.geometry().intersection(dzialki_lyr_feats.geometry())
if not domiar_wyciety.isEmpty():
feature = QgsFeature(wynikowy_linia_lyr.fields())
feature.setGeometry(domiar_wyciety)
feature["id"] = wynikowy_linia_lyr.featureCount() + 1
try:
dn = int(woda_lyr_feats['dn'])
except:
dn = 0
feature["dn"] = dn
wynikowy_linia_lyr.dataProvider().addFeatures([feature])
Wynikowo mamy:
Obiekty mają przepisane dn – średnice (wynik informacji na screenie powyżej). A poniżej cały kod…
woda_lyr = QgsProject.instance().mapLayersByName('woda')[0]
kanal_lyr = QgsProject.instance().mapLayersByName('kanal')[0]
dzialki_lyr = QgsProject.instance().mapLayersByName('dzialki_temp')[0]
uri_poligon_lyr = "polygon?crs=epsg:2178&field=id:integer"
wynikowy_poligon_lyr = QgsVectorLayer(uri_poligon_lyr, "wynikowy_poligon", "memory")
uri_linia_lyr = "linestring?crs=epsg:2178&field=id:integer&field=dn:string"
wynikowy_linia_lyr = QgsVectorLayer(uri_linia_lyr, "wynikowy_linia", "memory")
QgsProject.instance().addMapLayer(wynikowy_poligon_lyr)
QgsProject.instance().addMapLayer(wynikowy_linia_lyr)
aggr_buf = None
for woda_lyr_feats in woda_lyr.selectedFeatures():
try:
dn = int(woda_lyr_feats['dn'])
except:
dn = 0
zakres_m = 1+dn/2000
bufor_temp = woda_lyr_feats.geometry().buffer(zakres_m,5)
if aggr_buf is None:
aggr_buf = bufor_temp
else:
aggr_buf = aggr_buf.combine(bufor_temp)
for kanal_lyr_feats in kanal_lyr.selectedFeatures():
try:
dn = int(kanal_lyr_feats['dn'])
except:
dn = 0
zakres_m = 1+dn/2000
bufor_temp = kanal_lyr_feats.geometry().buffer(zakres_m,5)
if aggr_buf is None:
aggr_buf = bufor_temp
else:
aggr_buf = aggr_buf.combine(bufor_temp)
aggr_buf_dzialka = None
if aggr_buf is not None:
for dzialki_lyr_feats in dzialki_lyr.selectedFeatures():
aggr_buf_dzialka = aggr_buf.intersection(dzialki_lyr_feats.geometry())
if not aggr_buf_dzialka.isEmpty():
feature = QgsFeature(wynikowy_poligon_lyr.fields())
feature.setGeometry(aggr_buf_dzialka)
feature["id"] = wynikowy_poligon_lyr.featureCount() + 1
wynikowy_poligon_lyr.dataProvider().addFeatures([feature])
for dzialki_lyr_feats in dzialki_lyr.selectedFeatures():
for kanal_lyr_feats in kanal_lyr.selectedFeatures():
domiar_wyciety = kanal_lyr_feats.geometry().intersection(dzialki_lyr_feats.geometry())
if not domiar_wyciety.isEmpty():
feature = QgsFeature(wynikowy_linia_lyr.fields())
feature.setGeometry(domiar_wyciety)
feature["id"] = wynikowy_linia_lyr.featureCount() + 1
try:
dn = int(kanal_lyr_feats['dn'])
except:
dn = 0
feature["dn"] = dn
wynikowy_linia_lyr.dataProvider().addFeatures([feature])
for woda_lyr_feats in woda_lyr.selectedFeatures():
domiar_wyciety = woda_lyr_feats.geometry().intersection(dzialki_lyr_feats.geometry())
if not domiar_wyciety.isEmpty():
feature = QgsFeature(wynikowy_linia_lyr.fields())
feature.setGeometry(domiar_wyciety)
feature["id"] = wynikowy_linia_lyr.featureCount() + 1
try:
dn = int(woda_lyr_feats['dn'])
except:
dn = 0
feature["dn"] = dn
wynikowy_linia_lyr.dataProvider().addFeatures([feature])
Na dodatkowe wyjaśnienie zasługuje …dataProvider(), QgsVectorLayer dziedziczy tą metodę z QgsMapLayer (https://qgis.org/pyqgis/3.44/core/QgsMapLayer.html#qgis.core.QgsMapLayer.dataProvider). Zwraca ona klasę Providera danej usługi (https://qgis.org/pyqgis/3.44/core/QgsDataProvider.html#qgis.core.QgsDataProvider), i już z tej klasy wywołujemy addFeatures(list). Da się to zrobić też bezpośrednio z poziomu QgsVectorLayer, ale trzeba by dodać startEditing() a potem zakończyć commit lub rollbaclk (patrz trzeci akapit specyfikacji https://qgis.org/pyqgis/3.44/core/QgsVectorLayer.html ). Podsumowując Ja wybrałem takie rozwiązanie, ale są też inne…
Czego nie zrobiliśmy to kompletna kontrola błędów. W trakcie generowania i wycinania może powstawać wiele rodzajów geometrii, niekoniecznie potrzebnych. Do automatyki przydałoby się wyrzucenie szczątkowych poligonów, geometrii punktowych i innych błędnych geometrii (artefaktów niedoskonałości obliczeń wektorowych).
Jak zawsze zapraszam do dyskusji i kontaktu…
