Skomplikowana (chyba) analiza przestrzenna działki - budynki

Cześć wszystkim,
Proszę was o pomoc w temacie wykonania analizy przestrzennej z uzupełnieniem atrybutów.
Zęby nie było, że idę na łatwiznę - walczę z użyciem ChatGpt i DEEPSeek już od tygodnia. Staram się wycisnąć z ejaja skrypt do wykonania w konsoli Pythona . Poległem. Jesteście ostatnią deską ratunku jaka przyszła mi do głowy. Opiszę stan:
Dane są tu: Dane Ewidencji Gruntów i Budynków - Biuletyn Informacji Publicznej Urząd Miejski w Koszalinie

Co chcę osiągnąć : chcę na warstwie budynki dodać pole id_dzialek które w wyniku analizy ma być uzupełnione o id_dzialki z warstwy działek. Chodzi o zidentyfikowanie działek na których stoi budynek co pozwoli na wykrycie przekroczenia granicy (jeśli identyfikatorów będzie 4.5. więcej niż jeden).

Warunki:

  1. każdy budynek musi mieć uzupełnione pole id_dzialek,
  2. jeśli budynek stoi (przecina obrysem) kilka działek to chcę mieć te identyfikatory działek w polu id_dzialek ale:
    2.1 jeśli budynek dotyka jakiejś granicy (przylega do niej, styka się nie ważne czy z jednej, dwu czy trzech stron) - to nie interesują mnie te inne działki niż działka główna (bo faktycznie budynek nie narusza granicy - klasyczna dobrze wybudowana szeregówka lub garaże)
    2.2 z uwagi na błędy ludzkie a przy analizowaniu danych z dużym powiększeniem widać to jak na dłoni (brak snapowania) należy przyjąć, że przekroczenie granicy do 5 cm (odległość krawędzi budynku od linii granicznej a nie obiektu od obiektu) - to należy założyć błąd ludzki (z resztą mieści się to w błędzie pomiaru geodezyjnego) i też interesuje mnie tylko działka główna
  3. separatory identyfikatorów w polu id_dzialek - przecinek
  4. przy każdym identyfikatorze w polu id_dzialek, innym niż identyfikator działki głównej - powierzchnia zajęta przez budynek na obszarze tej działki co do której doszło do naruszenia granicy,
    Póki co nie jestem w stanie tego wymyśleć . HELP!

Spróbuj tego:

from qgis.PyQt.QtCore import QVariant
from qgis.core import QgsSpatialIndex

# Pobierz warstwy działki i budynki z projektu QGIS do zmiennych
dzialki = QgsProject.instance().mapLayersByName('EGB_ZP_DZIALKIPolygon')[0]
budynki = QgsProject.instance().mapLayersByName('EGB_ZP_BUDYNKIPolygon')[0]

# Dodaj pole 'id_dzialek' typu tekstowego do warstwy budynki, jeśli nie istnieje
def ensure_field(layer, field_name):
    fields = layer.fields()
    if fields.indexOf(field_name) == -1:
        from qgis.PyQt.QtCore import QVariant
        layer.startEditing()
        layer.addAttribute(QgsField(field_name, QVariant.String))
        layer.commitChanges()

ensure_field(budynki, 'id_dzialek')

# Utwórz spatial index dla działek
spatial_index = QgsSpatialIndex()
dzialka_features = {}
for dzialka in dzialki.getFeatures():
    spatial_index.insertFeature(dzialka)
    dzialka_features[dzialka.id()] = dzialka

# Przypisz działki do budynków
budynki.startEditing()
for i, budynek in enumerate(budynki.getFeatures()):
    if i > 100:
        break
    # Czytam geometrię budynku
    budynek_geom = budynek.geometry()
    # Buforuję geometrię budynku, aby uwzględnić ewentualne błędy we wrysowaniu w EGiB (do 5cm)
    budynek_geom_buf = budynek.geometry().buffer(-0.05, 5)
    # Twórzę tablicę na odnalezione działki
    id_dzialek = []
    # Wyznaczam bbox budynku i znajduję działki, które mogą zawierać budynek
    bbox = budynek_geom.boundingBox()
    dzialki_ids_temp = spatial_index.intersects(bbox)
    
    # Obliczam powierzchnie pokrycia dla wszystkich działek i sortuję
    dzialki_z_powierzchnia = []
    for dzialka_id in dzialki_ids_temp:
        dzialka = dzialka_features[dzialka_id]
        # Sprawdzam, czy budynek z buforem -0.05m jest całkowicie wewnątrz działki
        if dzialka.geometry().contains(budynek_geom_buf):
            # Działka zawiera cały budynek - największa możliwa powierzchnia
            powierzchnia = budynek_geom.area()
            dzialki_z_powierzchnia.append((powierzchnia, dzialka, True))  # True = zawiera cały budynek
        else:
            # Obliczam powierzchnię części budynku w działce
            czesc_w_dzialce = budynek_geom.intersection(dzialka.geometry())
            powierzchnia = czesc_w_dzialce.area()
            # Tylko działki z dodatnim pokryciem
            if not czesc_w_dzialce.isEmpty() and powierzchnia > 0:
                dzialki_z_powierzchnia.append((powierzchnia, dzialka, False))
    
    # Sortuję działki od największej powierzchni pokrycia
    dzialki_z_powierzchnia.sort(key=lambda x: x[0], reverse=True)
    
    # Przetwarzam posortowane działki
    for powierzchnia, dzialka, zawiera_caly in dzialki_z_powierzchnia:
        if zawiera_caly:
            # Jeśli działka zawiera cały budynek, dodaję jej ID i kończę
            id_dzialek.append(str(dzialka['IDENTYFIKA']))
            break
        else:
            # Dodaję ID działki z powierzchnią pokrycia
            powierzchnia_rounded = round(powierzchnia, 2)
            id_dzialek.append(f"{dzialka['IDENTYFIKA']}({powierzchnia_rounded})")
    # Ustawiam wartość atrybutu 'id_dzialek' w warstwie budynki
    budynki.changeAttributeValue(budynek.id(), budynki.fields().indexOf('id_dzialek'), ','.join(id_dzialek))
# Zatwierdzam zmiany w warstwie budynki
budynki.commitChanges()
print("Przypisanie działek do budynków zakończone.")

Bardzo koledze dziękuję!
Szacunek za wiedzę i chęć pomocy. Bardzo mi to pomogło. W życiu bym tego nie ogarnął.
Pozdrawiam z Koszalina

1 polubienie

Tam jest takie jedno założenie, które teoretycznie mogłoby być groźne. Chodzi o bufor -0.05m wykonany na geometrii budynku. Na podstawie życiowego doświadczenia przyjąłem zasadę, że żaden budynek nie ma mniej niż 11cm szerokości, więc nie spowoduje to “anihilacji” żadnej z geometrii budynków. Ale co do zasady najpierw należałoby to sprawdzić, albo pójść inną drogą :wink:

1 polubienie

Nie ma takich budynków, Doświadczenie życiowe dobrze podpowiada. Zauważyłem tam test na iteracji (100) obiektów. Zmodyfikowałem. Dodałem też meldunek po każdych 500 obiektach by obserwować postęp. Teraz wygląda to tak:

Jeszcze raz bardzo dziękuję. Na prawdę szacunek. :slight_smile:

from qgis.PyQt.QtCore import QVariant
from qgis.core import QgsSpatialIndex

Pobierz warstwy działki i budynki z projektu QGIS do zmiennych

dzialki = QgsProject.instance().mapLayersByName(‘EGB_ZP_DZIALKIPolygon’)[0]
budynki = QgsProject.instance().mapLayersByName(‘EGB_ZP_BUDYNKIPolygon’)[0]

Dodaj pole ‘id_dzialek’ typu tekstowego do warstwy budynki, jeśli nie istnieje

def ensure_field(layer, field_name):
fields = layer.fields()
if fields.indexOf(field_name) == -1:
from qgis.PyQt.QtCore import QVariant
layer.startEditing()
layer.addAttribute(QgsField(field_name, QVariant.String))
layer.commitChanges()

ensure_field(budynki, ‘id_dzialek’)

Utwórz spatial index dla działek

spatial_index = QgsSpatialIndex()
dzialka_features = {}
for dzialka in dzialki.getFeatures():
spatial_index.insertFeature(dzialka)
dzialka_features[dzialka.id()] = dzialka

Przypisz działki do budynków

budynki.startEditing()
for i, budynek in enumerate(budynki.getFeatures()):
if i % 500 == 0:
print(f"Przetworzono {i} budynków…")

# Czytam geometrię budynku
budynek_geom = budynek.geometry()
# Buforuję geometrię budynku, aby uwzględnić ewentualne błędy we wrysowaniu w EGiB (do 5cm)
budynek_geom_buf = budynek.geometry().buffer(-0.05, 5)
# Twórzę tablicę na odnalezione działki
id_dzialek = []
# Wyznaczam bbox budynku i znajduję działki, które mogą zawierać budynek
bbox = budynek_geom.boundingBox()
dzialki_ids_temp = spatial_index.intersects(bbox)

# Obliczam powierzchnie pokrycia dla wszystkich działek i sortuję
dzialki_z_powierzchnia = []
for dzialka_id in dzialki_ids_temp:
    dzialka = dzialka_features[dzialka_id]
    # Sprawdzam, czy budynek z buforem -0.05m jest całkowicie wewnątrz działki
    if dzialka.geometry().contains(budynek_geom_buf):
        # Działka zawiera cały budynek - największa możliwa powierzchnia
        powierzchnia = budynek_geom.area()
        dzialki_z_powierzchnia.append((powierzchnia, dzialka, True))  # True = zawiera cały budynek
    else:
        # Obliczam powierzchnię części budynku w działce
        czesc_w_dzialce = budynek_geom.intersection(dzialka.geometry())
        powierzchnia = czesc_w_dzialce.area()
        # Tylko działki z dodatnim pokryciem
        if not czesc_w_dzialce.isEmpty() and powierzchnia > 0:
            dzialki_z_powierzchnia.append((powierzchnia, dzialka, False))

# Sortuję działki od największej powierzchni pokrycia
dzialki_z_powierzchnia.sort(key=lambda x: x[0], reverse=True)

# Przetwarzam posortowane działki
for powierzchnia, dzialka, zawiera_caly in dzialki_z_powierzchnia:
    if zawiera_caly:
        # Jeśli działka zawiera cały budynek, dodaję jej ID i kończę
        id_dzialek.append(str(dzialka['IDENTYFIKA']))
        break
    else:
        # Dodaję ID działki z powierzchnią pokrycia
        powierzchnia_rounded = round(powierzchnia, 2)
        id_dzialek.append(f"{dzialka['IDENTYFIKA']}({powierzchnia_rounded})")
# Ustawiam wartość atrybutu 'id_dzialek' w warstwie budynki
budynki.changeAttributeValue(budynek.id(), budynki.fields().indexOf('id_dzialek'), ','.join(id_dzialek))

Zatwierdzam zmiany w warstwie budynki

budynki.commitChanges()
print(“Przypisanie działek do budynków zakończone.”)

Prawda, zapomniałem usunąć ten test :wink: Nie ma za co, całkiem przyjemna łamigłówka to była :slight_smile: