Na przykladzie danych z powiatu Krakowskiego wyszukamy i zaznaczymy obiekty w QGIS.
Najpierw potrzebujemy danych startowych, możemy je znaleźć w https://epodgik.pl/ a dokladniej w https://wms02.epodgik.pl/walidator/uslugi.php?wfs#
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
Pełna dokumentacja QGIS Python dostępna jest na stronie PyQGIS Developer Cookbook — QGIS Documentation documentation. Ja z tego korzystam z doskoku, tj wyszukuje to czego aktualnie potrzebuje.
Jest wiele sposobów na identyfikacje warstw w QGIS. Ja używam QgsProject (https://qgis.org/pyqgis/3.40/core/QgsProject.html) i QgsMapLayerStore.
# 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()}")
to zwróci nam coś w stylu
layerid = ms_budynki_7ae65efe_c42b_4471_8841_540e71af475b
layersrc = pagingEnabled='default' preferCoordinatesForWfsT11='false' restrictToRequestBBOX='1' srsname='EPSG:2178' typename='ms:budynki' url='https://wms.powiat.krakow.pl:1518/iip/ows' version='auto', layername = ms:budynki
layerid = ms_dzialki_ad325524_0958_4c2e_94e7_0947075a76b3
layersrc = pagingEnabled='default' preferCoordinatesForWfsT11='false' restrictToRequestBBOX='1' srsname='EPSG:2178' typename='ms:dzialki' url='https://wms.powiat.krakow.pl:1518/iip/ows' version='auto', layername = ms:dzialki
layerid = msbudynki_66ce0e5c_421b_447e_8b4d_be6840885e5e
layersrc = D:/projekt1/budynki.gpkg|layername=msbudynki, layername = budynki — msbudynki
layerid = msdzialki_7cacff14_ce63_4c32_9191_3788eccf239d
layersrc = D:/projekt1/dzialki.gpkg|layername=msdzialki, layername = dzialki — msdzialki
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)
Zwróci to nam 2 znalezione warstwy ,na których będziemy wyszukiwać, pełna specyfikacja : https://qgis.org/pyqgis/3.40/core/QgsMapLayer.html
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.
Dla przykładu przyjmijmy, że interesują nas działki mniejsze od 10 arów, zobaczymy ile takich znajdziemy. Iterujemy liste features i sprawdzamy czy powierzchnia geometrii danej feature jest mniejsza lub równa 1000m2 (patrz specyfikacja https://qgis.org/pyqgis/3.40/core/QgsFeature.html i https://qgis.org/pyqgis/3.40/core/QgsGeometry.html)
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ę 🙂
Dodaj komentarz