Skomplikowana (chyba) analiza przestrzenna działki - budynki

Centrum pomocy
1

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!
2

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.")
3

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
4

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
5

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.”)

6

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