🚀

Automatyzacja inicjalizacji projektu

dawid.huczynski.pl/lectures/project-startup

🔥 Cel

• Świadomość
• Dobre praktyki
• Przykłady zastosowania

🔥 Program

1. Wprowadzenie do GIS
2. Problem
3. Dobre praktyki
4. Zbieranie danych
5. Zarządzanie danymi
6. Pomocne biblioteki

🌍 1. Wprowadzenie do GIS

Projekcje map

Co nas zwodzi?

• engaging-data.com: projekcje map
• thetruesize.com

Układy wysokościowe

• PL-KRON86-NH – najpóźniej do grudnia 2023
• PL-EVRF2007-NH

Geodezyjne układy odniesienia w Polsce

będące matematyczną i fizyczną realizacją
europejskiego ziemskiego systemu odniesienia ETRS89,
oznaczone symbolami:

• PL-ETRF2000
• PL-ETRF89

Układy współrzędnych płaskich prostokątnych

• PL-LAEA
• PL-LCC
• PL-UTM
• PL-2000 – 15ºE, 18ºE, 21ºE i 24ºE, oznaczane numerami 5, 6, 7 i 8
• PL-1992

Strefy

• EPSG:2180 – cała Polska, PL-1992
• EPSG:2176 – PL-2000/15
• EPSG:2177 – PL-2000/18
• EPSG:2178 – PL-2000/21
• EPSG:2179 – PL-2000/24

Rodzaje danych

• rastrowe
• wektorowe:
  • punkty
  • linie
  • poligony
• bazy (tabele)

Usługi

• XYZ
• WMS / WMTS
• WFS

🙈 🙉 🙊 2. Zdefiniujmy problem

🛫 Problemy przy starcie

• Świadomość informacji
• Powtarzalność zbierania danych
• Powtarzalność przetwarzania danych
• Brak spójności
• Czas
• Przerost informacji
• Aktualizacja

💡 3. Dobre praktyki

Jednostki

CAD czy BIM: zawsze ustawmy jednostki w pliku.

Dostajemy plik: sprawdzamy i ewentualnie uzupełniamy jednostki.

Wydajemy plik: sprawdzamy i uzupełniamy jednostki.

TERYT

226101_1.0303.769/16

województwo · powiat · gmina _ rodzaj . obręb . nr działki

Sześcian koordynacyjny

• lokalizacja – punkt bazowy (a nawet 8)
• skala – widać różnice w każdym kierunku
• jednostki – 1 m × 1 m × 1 m
• każdy budynek w jednym miejscu (ściany pionowe)
• każda branża w budynku (od góry)
• przyciąganie, wyrównywanie, skalowanie

Optymalizacja pliku startowego

• podział na mniejsze części
• kontrola zmian
• minimum niezbędnych informacji
• aktualizacje

Co mogę zrobić?

Jeżeli mogę, staram się zmniejszyć liczbę rodzin i elementów.

Jeżeli poziom skomplikowania rośnie, lepiej rozłożyć to na części i szybko aktualizować proste rodziny.

Proste jest lepsze niż złożone.
Złożone jest lepsze niż skomplikowane.

Zen of Python

Śledzenie zmian

• git
• śledzenie zmian użytkowników
• IFC

Śledzenie zmian

Śledzenie zmian

IFC a zwykły tekst

⛏️ 4. Zbieranie danych

Działka

Działka

https://uldk.gugik.gov.pl/?request=GetParcelByXY&xy={x},{y},{epsg}
https://uldk.gugik.gov.pl/?request=GetParcelById&id={teryt}
https://uldk.gugik.gov.pl/?request=GetAggregateArea&id={teryt},{teryt},{teryt}

Działka

def get_plot(project_id, dzialki=None, x=None, y=None, epsg=None, target_epsg=None):
    # https://uldk.gugik.gov.pl/opis.html#
    if epsg and x and y:  # przez punkt
        req = requests.get(
            f"https://uldk.gugik.gov.pl/?request=GetParcelByXY&xy={x},{y},{epsg}")
    elif isinstance(dzialki, str):  # przez pojedynczą działkę, np. '141201_1.0001.1867/2'
        req = requests.get(
            f"https://uldk.gugik.gov.pl/?request=GetParcelById&id={dzialki}")
    elif isinstance(dzialki, list):  # przez listę działek
        req = requests.get(
            f'https://uldk.gugik.gov.pl/?request=GetAggregateArea&id={",".join(dzialki)}')
    else:
        raise ValueError("Źle wprowadzone dane.")

    wkb = req.text.split("\n")[1]
    shape = shapely.from_wkb(wkb)
    if target_epsg:
        epsg = epsg if epsg else 2180
        project = pyproj.Transformer.from_crs(epsg, target_epsg, always_xy=True).transform
        shape = transform(project, shape)
    buffer = shape.buffer(200)

Chmura punktów

Zapytanie: jakie chmury są w danej okolicy

URL = "https://mapy.geoportal.gov.pl/wss/service/PZGIK/DanePomiaroweNMT/WMS/"
x = point.x  # koordynaty
y = point.y  # koordynaty
url = URL + "SkorowidzeUkladEVRF2007?"
LAYERS = [
    # 'Układ wysokościowy PL-EVRF2007-NH',
    "SkorowidzeLIDAR2023",
    "SkorowidzeLIDAR2022",
    "SkorowidzeLIDAR2021",
    "SkorowidzeLIDAR2020iStarsze",
]

Chmura punktów

Zapytanie: jakie chmury są w danej okolicy

PARAMS = {
    "SERVICE": "WMS",
    "request": "GetFeatureInfo",
    "version": "1.3.0",
    "layers": ",".join(LAYERS),
    "styles": "",
    "crs": f"EPSG:{epsg}",
    "bbox": "%f,%f,%f,%f"
    % (y - distance, x - distance, y + distance, x + distance),
    "width": "101",
    "height": "101",
    "format": "image/png",
    "transparent": "true",
    "query_layers": ",".join(LAYERS),
    "i": "50",
    "j": "50",
    "INFO_FORMAT": "text/html",
}
r = requests.get(url=url, params=PARAMS)

Teren / toposurface

Budynki 3D

url = "http://mapy.geoportal.gov.pl/wss/service/PZGIK/PRG/WFS/AdministrativeBoundaries?language=pol&"
distance = 50
PARAMS = {
    "SERVICE": "WFS",
    "request": "GetFeature",
    "version": "2.0.0",
    "typenames": "ms:A02_Granice_powiatow",
    "crs": f"EPSG:{epsg}",
    "bbox": "{:f},{:f},{:f},{:f}".format(y - distance, x - distance, y + distance, x + distance),
    "FORMAT": "application/gml+xml",
    "WIDTH": 101,
    "HEIGHT": 101,
    "i": 50,
    "j": 50,
}
r = requests.get(url=url, params=PARAMS)

Budynki 3D

url = f"https://opendata.geoportal.gov.pl/InneDane/Budynki3D/LOD1/2022/{id_powiatu[:2]}/{id_powiatu}.zip"

WMS / WFS

• Ortofotomapa
• Sieci
• MPZP – wektor / obraz
• Drogi – wektor
• BDOT – wektor

WMS / WFS

Powiaty a geoportal

integracja.gugik.gov.pl/eziudp

WMS / WFS

WMS / WFS

WMS / WFS

🧹 5. Zarządzanie danymi

shapely

shape = get_plot(project_id, dzialki=dzialki)
buffer = shape.buffer(200)
center = shape.centroid

shapely

shapely

shapely

• polygon.intersection(other_polygon)
• polygon.difference(other_polygon)
• polygon.union(other_polygon)
• polygon.buffer(distance)

geopandas

zarządzanie danymi

geopandas

laspy / PDAL

• przetwarzanie chmur punktów i plików las/laz
• filtrowanie
• geolokalizacja
• transpozycja
• reprojekcja – przekształcenie
• klasyfikacja

laspy / PDAL

def join_and_split_point_clouds(las_fp, shape, buffer, proj_id):
    points, header = merge_clouds(las_fp)
    las = laspy.LasData(header)
    las.points = points
    clouds = {}
    buffer_mask = get_site_mask(points, buffer)
    buffer_points = las.points[buffer_mask]
    clouds["buffer-teren"] = buffer_points[buffer_points.classification == 2]
    print("Filtered buffer-teren")
    site_mask = get_site_mask(points, shape)
    dzialka_points = las.points[site_mask]
    buffer_without_site_mask = np.array([buffer_mask, ~site_mask]).min(axis=0)
    clouds["buffer"] = las.points[buffer_without_site_mask]
    print("Filtered buffer")
    cl = dzialka_points.classification
    clouds["dzialka-teren"] = dzialka_points[cl == 2]
    print("Filtered dzialka teren")
    zielen_mask = np.array([cl == 3, cl == 4, cl == 5]).max(axis=0)
    clouds["dzialka-zielen"] = dzialka_points[zielen_mask]
    clouds["dzialka-budynki"] = dzialka_points[cl == 6]
    print("Filtered dzialka budynki dzialka zielen")
    fpaths = []
    for x in clouds:
        las.points = clouds[x]
        fp = os.path.join("output", f"{proj_id}-{x}.laz")
        las.write(fp)
        fpaths.append(fp)
    return fpaths, clouds, las

laspy / PDAL

laspy / PDAL

laspy / PDAL

laspy / PDAL

laspy / PDAL

laspy / PDAL

📚 6. Pomocne biblioteki

• shapely
• pandas
• geopandas
• pyproj
• laspy

shapely

• 2D
• bufor / offset
• wycinanie / przenikanie
• operacje na poligonach
• import / eksport różnych formatów 2D

pandas

• Excel na sterydach
• DataFrame
• dowolne dane
• zaawansowane i automatyczne filtrowanie
• import / eksport do Excela, CSV itd.

🫂 Dzięki za wytrwałość

Zapraszam do kontaktu

dawid.huczynski.pl · [email protected]