Zaloguj się

Tag: SŁUŻEBNOŚĆ PRZESYŁU

  • QGIS i automatyzacja tworzenia mapy do służebności przesyłu

    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&amp;field=id:integer"
wynikowy_poligon_lyr = QgsVectorLayer(uri_poligon_lyr, "wynikowy_poligon",  "memory")

uri_linia_lyr = "linestring?crs=epsg:2178&amp;field=id:integer&amp;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…

  • QGIS i tworzenie mapy do slużebności przesyłu

    Dzisiaj trochę inny model podejścia do problemu. Jak już pisałem w poprzednich wpisach Polska nie publikuje w formie wektorowej obiektów, którymi chciałem się dzisiaj zająć (czyli uzbrojenia terenu). Stworzyłem więc swój własny obszar prac. Działki zostały wycięte w losowym miejscu mapy udostępnianym w ramach https://epodgik.pl/, natomiast sieć wodociągową i kanalizacyjną „wymyśliłem” sam na potrzeby przedstawienia rozwiązania problemu. Chciałbym zaznaczyć, że dane do „prawdziwych” służebności przesyłu powinny wywodzić się z danych urzędowych, a niniejszy przykład przedstawia demonstracyjne rozwiązanie problemu.

    W załączeniu plik z danymi do dzisiejszej zabawy : projekt z danymi (spakowany 7zip) (fix: dodałem brakującą warstwę do pliku)

    Mamy tu 3 warstwy z czego dzialki_temp zawiera pola ID i NR_DZIALKI, a warstwy woda i kanal zawieraja ID i DN (srednice sieci).

    Do stworzenia załącznika do służebności przesyłu potrzebujemy

    • pas i jego powierzchnie
    • długości i średnice sieci objętych

    Zajmijmy się pasem, dla przykładowej działki :

    Zaznaczamy wszystkie obiekty na warstwie woda i kanał które dotykają działki. Z mojego doświadczenia lepiej zaznaczyć więcej sieci żeby nie wynikły problemy, że nagle pas na granicy działki zmniejsza się(powstaje łuk), chociaż wiemy że sieć ta biegnie dalej i pas powinien być równo oddalony po osi aż do granicy.

    Dla przykładu zaznaczymy przez lokalizację wszystko w warstwie kanal co przecina dzialki_temp (tylko zaznaczone) – 5 features selected, tak samo robimy dla warstwy woda, także 5 features selected, czyli mamy po 5 obiektów wodnych i kanalizacyjnych.

    Stwórzmy teraz pasy wokół zaznaczonych sieci. Załóżmy że pas ma mieć metr + połowę średnicy sieci od osi przewodu (jak w niektórych urzędach się przyjmuje). Użyjemy do tego narzędzia Bufor

    Wyjaśnijmy parametry po kolei.

    • najpierw warstwa wejściowa, tutaj kanał,
    • potem powinniśmy zaznaczyć „Tylko zaznaczone obiekty”, tutaj mamy dokładnie wszystkie obiekty zaznaczone wiec nie robi to różnicy.
    • odległość, tu możemy podać jedną stałą odległość bufora, my użyjemy bardziej zaawansowanego podejscia i klikniemy w pole Edytuj…

    i wpisujemy :

    Już wyjaśniam dlaczego tak robimy. Możemy wymusić aby dla każdego obiektu dla którego będzie powstawał bufor, była to zmienna a nie stała. Tutaj jest to:

    dla przykładu DN = 100:

    • Potrzebujemy średnice wyrażoną w metrach zatem =100mm / 1000 = 0.1m
    • połowa średnicy czyli (/ 2) = 0.05m
    • dodajmy metr od krawędzi przewodu czyli 0.05m + 1 m
    • wynikowo dla DN=100 bufor będzie wynosił 1.05m

    Klikamy Ok

    Wróćmy do Bufora

    • Pole Odcinki określa stopień zaokrąglenia bufora na końcach(zostawmy 5).
    • Styl zakończenia „zaokrąglony” – chcemy mieć równo odsunięty bufor na końcach
    • Styl połączenia „zaokrąglony” – chcemy mieć równo odsunięty bufor na łączeniach
    • Limit fazy (uciosu) – nie mamy połączeń ostrych więc nas nie interesuje
    • Zaznaczamy Agreguj wyniki, spowoduje to powstanie jednego połączonego obiektu.

    Wybieramy Uruchom

    Otrzymaliśmy 1 tymczasowy bufor dla kanału.

    Tak samo skonfigurujemy to dla wody i otrzymujemy:

    Łączymy teraz warstwy, albo narzędziem albo kopiując obiekty z jednej warstwy do drugiej. Ja użyję złączenia:

    Klikamy … i zaznaczamy warstwy Bufor i Bufor, Klikamy Uruchom, Powstała nam warstwa Złączone którą przenosimy na samą górę, odznaczamy też niepotrzebne warstwy Bufor.

    używamy teraz Narzędzia Agreguj

    Da nam to warstwę zagregowaną

    Pozostaje tylko przyciąć do granic działki, narzędzie Przytnij:

    Wybieramy kolejno:

    • Warstwa wejsciowa : Zgrupowane obiekty
    • Warstwa nakładki : dzialki_temp
    • Zaznaczamy „Tylko zaznaczone obiekty” – bo chcemy tylko dla zaznaczonej działki wycięcie
    • I klikamy Uruchom…

    Po odznaczeniu niepotrzebnych warstw mamy warstwę Przycięte – nasz pas dla działki.

    Teraz spróbujmy przyciąć sieci:

    Przechodzimy na warstwę woda, wybieramy znowu Przytnij i tym razem:

    • Warstwa wejsciowa : woda
    • Zaznaczamy : Tylko zaznaczone obiekty
    • Warstwa nakładki: dzialki_temp
    • Zaznaczamy „Tylko zaznaczone obiekty” – bo chcemy tylko dla zaznaczonej działki wycięcie
    • I znowu Uruchom…

    Tak samo robimy dla warstwy kanal

    Odznaczamy warstwę woda i kanał, Odznaczamy też wszystkie obiekty zaznaczone. Mamy w projekcie 3 warstwy Przycięte, warstwę pasa i 2 warstwy domiarów odrębnie dla kanału i wody

    Pozostaje nam skonfigurować wyglad i etykiety i wynikowo otrzymamy:

    Czego się dzisiaj nauczyliśmy, co zrobiliśmy:

    • zaznaczyliśmy po lokalizacji warstwy liniowe
    • stworzyliśmy i zagregowaliśmy bufory poszczególnych warstw liniowych
    • docięliśmy bufory do wymaganego zakresu
    • wynikowo powstał zagregowany pas służebności ze zmienną długością pasa dla 2 odrębnych warstw.

    W przyszłych wpisach pokażę jak zautomatyzować to za pomocą Pythona.

    Jak zawsze zapraszam do dyskusji…