Dziś wyjaśnię jak uruchomić QGIS z poziomu Pythona bez włączania GUI (Graphical User Interface) QGIS. Brzmi to dziwnie, ale jest przydatne gdy robimy coś co ma działać w tle, być uruchamiane regularnie jako task, wykonywać złożone analizy nie wymagające gui QGIS.
Zacznijmy od konfiguracji Pythona w Visual Studio, w końcu na czymś pisać musimy.
Najpierw tworzymy klasyczną aplikacją Pythona, nazwałem to qgisnogui.
Utworzyło nam pusty projekt, musimy teraz dodać naszą instalację Pythona, tj Python Environment z zainstalowanego QGIS-a. Klikamy prawym klawiszem myszy na Python Environments, Add Environment…
Potem wybieramy Existing Environment i wybieramy <Custom>
Rozwinie się dialog, musimy tylko wskazać lokalizację PYTHONA (Prefix path) w folderze QGIS, reszta powinna się sama wypełnić…
Klikamy Add. Nasze nowo utworzona konfiguracja środowiska się automatycznie zmieni w projekcie.
Teraz tylko musimy jeszcze dodać do Search Paths folder z bibliotekami QGIS pythona. Prawy klawisz na Search Paths. Add Folder to Search Path i u mnie to był : C:\qgis3421\apps\qgis\python.
No i teraz już można uruchomić QGIS-a z poziomu pythona (i tutaj cenna uwaga, nie do wszystkiego qgisapp jest potrzebna, część bibliotek pracuje zupełnie niezależnie, np jak chcemy tylko „obrobić” geometrię to wystarczy sam Python i api qgis, ale jak już chcemy skorzystać z providera, załadować warstwę, wydrukować pdf, a może nawet napisać własną nakładkę gui bazującą na api QGIS to już musimy zainicjować aplikację qgis)
Dla przykładu prosty kod :
# zacznijmy od czegoś prostego czyli zdefiniowania geometrii
from qgis.core import QgsGeometry
geom1 = QgsGeometry.fromWkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))")
print(geom1.area())
po kliknięciu uruchom oblicza nam powierzchnię poligonu.
Teraz uruchommy QGIS i załadujmy warstwę:
# uruchamiamy qgisapp bez gui, ale z możliwością korzystania z jego funkcji
import sys
from qgis.core import QgsApplication
# Inicjalizacja aplikacji QGIS bez GUI
qgs = QgsApplication([], False)
qgs.initQgis()
# Teraz możemy korzystać z funkcji QGIS, np. wczytać warstwę
from qgis.core import QgsVectorLayer
layer = QgsVectorLayer(r"D:\dane\dzialki_krakow_2180.shp", "dzialki", "ogr")
if not layer.isValid():
print("Layer failed to load!")
print(f"Ilość obiektów : {layer.featureCount()}")
i = 0
for feature in layer.getFeatures():
i+=1
print(f"Obiekt {feature.id()}, ID_DZIALKI = {feature['ID_DZIALKI']}")
if i > 10:
break
No i wynik :
Tutaj jeszcze cenna uwaga dla początkujących, jak używamy polskich znaków w Pythonie to pamiętajmy o tym aby przed uruchomieniem kodu zrobić Save As…, potem Save with Encoding… i wybrać UTF-8 (65001) 🙂
Na tym zakończymy dzisiejszy wywód, pozdrawiam, przepraszając równocześnie za długą nieobecność. Praca praca i praca, ale będzie o czym pisać w przyszłości 🙂
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 :
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.
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.
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(…
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.
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])
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).
Kopiujemy adres WFS z Usługi pobierania i dodajemy w warstwach QGIS / WFS-OGC
Dodajemy 2 przykładowe warstwy
Eksportujemy je np. do gpkg.
Wynikowo mamy wyeksportowane 2 warstwy, będziemy wyszukiwać w konsoli pythona. Reszte odznaczymy bo serwisy online są nieprzewidywalne, tj teraz są, za chwilę może ich nie być. Zapisałem też projekt bo QGIS też czasem potrafi się zwiesić bez powodu.
Teraz przechodzimy do konsoli pythona w QGIS,
Menu -> Wtyczki -> Konsola Pythona -> w otwartej konsoli klikamy ikonkę Pokaż Edytor.
Konsola Pythona składa się z kilku obszarów:
ŻÓŁTY to Output Pythona czyli tu będziemy wypisywać wyniki jak i błędy naszego kodu
RÓŻOWY to input pojedynczych komend Pythona, tu możemy wpisać komendę która ma być uruchomiona po naciśnięciu Enter
ZIELONY (edytor otwieramy poprzez naciśniecie przycisku zaznaczonego na fioletowo) w tym edytorze możemy napisać już złożony skrypt, uruchamiamy go poprzez naciśniecie przycisku PLAY zaznaczonego na niebiesko
# to zwraca nam zainicjowaną klasę aktualnie
# otwartego projektu https://qgis.org/pyqgis/3.40/core/QgsProject.html
qgsprj = QgsProject.instance()
# to zwraca nam klase
# https://qgis.org/pyqgis/3.40/core/QgsMapLayerStore.html
qgsls = qgsprj.layerStore()
layers = qgsls.mapLayers() # lista niepowtarzalnych id warstw
for layer in layers: # iterujemy listę niepowtarzalnych id warstw
# potrzebujemy wyszukać specyficzną warstwę
# wiele warstw może mieć podobne nazwy, źródła etc
# musimy byc pewni ze wyszukalismy specyficzną warstwę
# przyglądnijmy się danym w klasie QgsMapLayer i wróćmy do QgsProject
# https://qgis.org/pyqgis/3.40/core/QgsMapLayer.html
print(f"layerid = {layer}")
print(f"layersrc = {qgsprj.mapLayer(layer).source()}, layername = {qgsprj.mapLayer(layer).name()}")
nam chodzi tylko o nasze gpkg czyli layersrc in 'dzialki.gpkg’ i layerbname = 'dzialki – msdzialki’ oraz layersrc in 'budynki.gpkg’ and layername = 'budynki – msbudynki’, zmodyfikujmy nasz skrypt :
#zmienna lyr budynki
budynki = None
#zmienna lyr dzialki
dzialki = None
# to zwraca nam zainicjowaną klasę aktualnie
# otwartego projektu https://qgis.org/pyqgis/3.40/core/QgsProject.html
qgsprj = QgsProject.instance()
# to zwraca nam klase
# https://qgis.org/pyqgis/3.40/core/QgsMapLayerStore.html
qgsls = qgsprj.layerStore()
layers = qgsls.mapLayers() # listę warstw (obie)
for layer in layers: # iterujemy listę niepowtarzalych id warstwy
# potrzebujemy wyszukać specyficzną warstwę
# wiele warstw może mieć podobne nazwy
# musimy byc pewni ze wyszukalismy specyficzną warstwę
# przyglądnijmy się danym w klasie QgsMapLayer i wróćmyu do QgsProject
# https://qgis.org/pyqgis/3.40/core/QgsMapLayer.html
print(f"layerid = {layer}")
print(f"layersrc = {qgsprj.mapLayer(layer).source()}, layername = {qgsprj.mapLayer(layer).name()}")
# jeżeli spelniaja nasze kryteria, przepisz zmienna budynki, dzialki
if 'budynki.gpkg' in qgsprj.mapLayer(layer).source() and qgsprj.mapLayer(layer).name() == 'budynki — msbudynki':
budynki = qgsprj.mapLayer(layer)
if 'dzialki.gpkg' in qgsprj.mapLayer(layer).source() and qgsprj.mapLayer(layer).name() == 'dzialki — msdzialki':
dzialki = qgsprj.mapLayer(layer)
print(dzialki)
print(budynki)
najpierw wyszukamy po danych tj. uzyjemy selectByExpression, przygladnijmy się danym które oferuje nam warstwa dzialki
Wyszukajmy wszystkie dzialki które są w Giebułtowie i działka zaczyna się od 434, czyli expression musi brzmieć : NAZWA_OBREBU = 'Giebułtów’ AND NUMER_DZIALKI LIKE '434%’, bazując na qgisowej wersji zaznaczania będzię to ok 51 obiektów, ale sprawdźmy jak to będzie działać w Pythonie.
dzialki.selectByExpression('NAZWA_OBREBU = \'Giebułtów\' AND NUMER_DZIALKI LIKE \'434%\'')
Jest to najprostszy sposób zaznaczania. W QGIS oczywiscie mozemy pobrać features bez zaznaczania i mieć je w iteratorze, liście etc. do tego celu musimy użyć funkcji getFeatures i stworzyć QgsFeatureRequest np:
features = dzialki.getFeatures(QgsFeatureRequest(QgsExpression('NAZWA_OBREBU = \'Giebułtów\' AND NUMER_DZIALKI LIKE \'434%\'')))
flist = list(features)
print(len(flist))
Da nam to jak już mówiłem około 51 features selected. Powód dla którego przetwarzam to na listę jest prosty, iteratory można użyć tylko raz, a listę wiele razy. Możemy teraz zrobic liste feature id i wyselekcjonować recznie wybrane id.
fid = [] # pusta lista samych id
flist2 = [] # lista samych wybranych features
for f in flist:
if f.geometry().area() <= 1000:
fid.append(f.id())
flist2.append(f)
dzialki.selectByIds(fid)
Wynikowo mamy wyselekcjonowane po Id wszystkie obiekty mniejsze lub równe 1000m2
No a po co była nam 2 warstwa Budynki ? Zaznaczmy teraz wszyskie budynki znajdujące się na wyselekcjonowanych działkach. Aby tego dokonać i zrobić to wydajnie musimy wykonać kilka kroków.
Po pierwsze żeby nie porównywać wszystkich obiektów na warstwie Budynki potrzebujemy pobrać tylko obiekty w zakresie features. Ja preferuje uzycie funkcji boundingBox ktora tworzy prostokąt z którego zostaną pobrane dane od providera (dostawcy warstwy tutaj ogr). Najczesciej providerzy mają zoptymalizowaną funkcję pobierania features z prostokątnych obszarów. Potem z wybranych juz features sprawdzamy czy pojedyncze features dotykają naszego obiektu geometrycznego. Moglibyśmy iterować wszystkie obiekty na warstwie budynki ale będzie to mało wydajne (czytaj bardzo wolne). Finalny kod.
#zmienna lyr budynki
budynki = None
#zmienna lyr dzialki
dzialki = None
# to zwraca nam zainicjowaną klasę aktualnie
# otwartego projektu https://qgis.org/pyqgis/3.40/core/QgsProject.html
qgsprj = QgsProject.instance()
# to zwraca nam klase
# https://qgis.org/pyqgis/3.40/core/QgsMapLayerStore.html
qgsls = qgsprj.layerStore()
layers = qgsls.mapLayers() # listę warstw (obie)
for layer in layers: # iterujemy listę niepowtarzalych id warstwy
# potrzebujemy wyszukać specyficzną warstwę
# wiele warstw może mieć podobne nazwy
# musimy byc pewni ze wyszukalismy specyficzną warstwę
# przyglądnijmy się danym w klasie QgsMapLayer i wróćmyu do QgsProject
# https://qgis.org/pyqgis/3.40/core/QgsMapLayer.html
print(f"layerid = {layer}")
print(f"layersrc = {qgsprj.mapLayer(layer).source()}, layername = {qgsprj.mapLayer(layer).name()}")
# jeżeli spelniaja nasze kryteria, przepisz zmienna budynki, dzialki
if 'budynki.gpkg' in qgsprj.mapLayer(layer).source() and qgsprj.mapLayer(layer).name() == 'budynki — msbudynki':
budynki = qgsprj.mapLayer(layer)
if 'dzialki.gpkg' in qgsprj.mapLayer(layer).source() and qgsprj.mapLayer(layer).name() == 'dzialki — msdzialki':
dzialki = qgsprj.mapLayer(layer)
print(dzialki)
#dzialki.selectByExpression('NAZWA_OBREBU = \'Giebułtów\' AND NUMER_DZIALKI LIKE \'434%\'')
features = dzialki.getFeatures(QgsFeatureRequest(QgsExpression('NAZWA_OBREBU = \'Giebułtów\' AND NUMER_DZIALKI LIKE \'434%\'')))
flist = list(features)
print(len(flist))
fid = []
flist2 = []
# latwiej nam bedzie selekcjonowac 1 geometria niz 51
# wiec polaczymy to w jedna geom
geom = None
for f in flist:
if f.geometry().area() <= 1000:
fid.append(f.id())
flist2.append(f)
# tu sie dzieje magia łączenia
if geom is None:
geom = f.geometry()
else:
geom = geom.combine(f.geometry())
dzialki.selectByIds(fid)
# teraz w geom mamy polaczona geometrie wszystkich dzialek
# tworzymy boundingBox czyli prostokąt który zawiera całą geometrie geom
geombb = geom.boundingBox()
# no i pobieramu features z budynkow po boundingBox
feat_budynki = budynki.getFeatures(QgsFeatureRequest().setFilterRect(geombb))
fbudynki_id = []
for f in feat_budynki:
# teraz musimy sprawdzic czy rzeczywiscie geometrie sie przecinaja gdziekolwiek
if f.geometry().intersects(geom):
fbudynki_id.append(f.id())
# no i finalnie zaznaczamy budynki na mapie
budynki.selectByIds(fbudynki_id)
i finalny wynik w QGIS
32 features selected on layer budynki.
Podsumowując :
dodalismy warstwy wfs
wyeksportowalismy je
wyszukalismy potrzebne warstwy z poziomu Pythona
wyszukalismy dane po expression z warstwy dzialki
wyfiltrowalismy tylko wyniki spelniajace kryteria wielkosci obiektu
i na koniec zaznaczylismy tymi obiektami, obiekty z warstwy budynki
Zapraszam do kontaktu i komentowania, jestem też otwarty na konstruktywną krytykę 🙂
