MODELOVÁNÍ MORFOLOGIE KOMÁŘÍCH LÍHNIŠŤ POMOCÍ GEOSTATICKÝCH NÁSTROJŮ

import arcpy
import os
import arcpy.sa  # Spatial Analyst
import math
# externí knihovny pro .las parsing
import laspy
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import numpy as np

# Zajistit 3D Analyst pro práci s .las
arcpy.CheckOutExtension("3D")
arcpy.env.overwriteOutput = True
arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

# --- DEFINICE PARAMETRŮ SCRIPT TOOLU ---
arcpy.Parameter(
    displayName="Vstupní feature class nebo LAS dataset",
    name="input_features",
    datatype="DEFeatureClass",
    parameterType="Required",
    direction="Input"
)
arcpy.Parameter(
    displayName="Pole hodnot pro rastr",
    name="value_field",
    datatype="Field",
    parameterType="Required",
    direction="Input"
)
arcpy.Parameter(
    displayName="Velikost buňky rastru",
    name="cell_size",
    datatype="GPString",
    parameterType="Required",
    direction="Input"
)
arcpy.Parameter(
    displayName="Vrstva řek (volitelné)",
    name="rivers_layer",
    datatype="DEFeatureClass",
    parameterType="Optional",
    direction="Input"
)
arcpy.Parameter(
    displayName="Vrstva budov (volitelné)",
    name="buildings_layer",
    datatype="DEFeatureClass",
    parameterType="Optional",
    direction="Input"
)
arcpy.Parameter(
    displayName="Výstupní vrstva tůní filtrovaná podle konvexity",
    name="filtered_output",
    datatype="DEFeatureClass",
    parameterType="Derived",
    direction="Output"
)
arcpy.Parameter(
    displayName="Minimální rozměr plochy polygonu",
    name="minimal_plocha_tun",
    datatype="GPString",
    parameterType="Required",
    direction="Input"
)
arcpy.Parameter(
    displayName="Maximální rozměr plochy polygonu",
    name="maximal_plocha_tun",
    datatype="GPString",
    parameterType="Required",
    direction="Input"
)
arcpy.Parameter(
    displayName="Hodnota PEAK algoritmu (v metrech)",
    name="peak_al",
    datatype="GPString",
    parameterType="Required",
    direction="Input"
)

# PŘIDANÝ PARAMETR PRO MAZÁNÍ MEZIKROKOVÝCH VRSTEV
arcpy.Parameter(
    displayName="Smazat mezikrokové vrstvy?",
    name="delete_intermediate_data",
    datatype="GPBoolean",
    parameterType="Optional",
    direction="Input"
)
# --- NAČTENÍ HODNOT PARAMETRŮ ---
input_features        = arcpy.GetParameterAsText(0)
value_field           = arcpy.GetParameterAsText(1)
cell_size             = arcpy.GetParameterAsText(2)
rivers_layer          = arcpy.GetParameterAsText(6)
buildings_layer       = arcpy.GetParameterAsText(7)
peak_al               = arcpy.GetParameterAsText(5)
minimal_plocha_tuhne  = arcpy.GetParameterAsText(8)
maximalni_plocha_tuhne= arcpy.GetParameterAsText(9)
delete_intermediate_data = arcpy.GetParameter(10)
# Převod velikosti buňky na float
try:
    cell_size_float = float(cell_size.replace(",", "."))
except ValueError:
    arcpy.AddError(f"Nelze převést '{cell_size}' na číslo.")
    raise

# Funkce pro konverzi .las na GeoPackage filtrem ground body
def convert_las_to_gpkg(las_path, out_gpkg):
    ext = os.path.splitext(las_path)[1].lower()
    if ext != ".las":
        raise ValueError("Očekávám .las soubor")
    with laspy.open(las_path) as las:
        pts = las.read()
        mask = pts.classification == 2
        x = pts.x[mask]
        y = pts.y[mask]
        z = pts.z[mask]
        if x is None or y is None or z is None or len(x) == 0:
            raise ValueError("Žádné ground body v LAS souboru")
        # zajistíme, že máme jednorozměrné pole i pro scalar
        x_arr = np.atleast_1d(x)
        y_arr = np.atleast_1d(y)
        z_arr = np.atleast_1d(z)
        df = pd.DataFrame({
            "X": x_arr,
            "Y": y_arr,
            "Elevation": z_arr
        })
    df["X"] = pd.to_numeric(df["X"], errors="coerce")
    df["Y"] = pd.to_numeric(df["Y"], errors="coerce")
    geom = [Point(xx, yy) for xx, yy in zip(df.X.tolist(), df.Y.tolist())]
    gdf = gpd.GeoDataFrame(df, geometry=geom, crs="EPSG:32633")
    gdf.to_file(out_gpkg, driver="GPKG")

# Kontrola a případná konverze inputu
if input_features.lower().endswith('.las'):
    arcpy.AddMessage("Vstup je .las soubor – konvertuji na GPKG ground bodů...")
    # získáme absolutní cestu k .las souboru
    las_desc = arcpy.Describe(input_features)
    las_path = las_desc.catalogPath
    gpkg = os.path.join(arcpy.env.scratchFolder, 'ground_points.gpkg')
    try:
        convert_las_to_gpkg(las_path, gpkg)
    except Exception as e:
        arcpy.AddError(f"Chyba při konverzi LAS: {e}")
        raise
    layer = arcpy.MakeFeatureLayer_management(gpkg, 'ground_lyr')
    input_features = layer
    value_field = 'Elevation'
else:
    arcpy.AddMessage("Vstup není .las – pokračuji s feature class...")

# Nastavení workspace a tvorba DMR rastru a tvorba DMR rastru
input_workspace = os.path.dirname(arcpy.Describe(input_features).catalogPath)
arcpy.AddMessage(f"Výstupní workspace: {input_workspace}")
DMR_raster = "DMR_RASTER"
out_raster = os.path.join(input_workspace, DMR_raster)

for i in range(3):
    arcpy.AddMessage(f"Parametr {i}: {arcpy.GetParameterAsText(i)}")
 
# Výpis parametrů 6 a 7 (řeky a budovy)
for i in [5, 6, 7, 8, 9]:
    arcpy.AddMessage(f"Parametr {i}: {arcpy.GetParameterAsText(i)}")
# Definice proměnných pro cesty a názvy polí
input_workspace = ""
output_raster = None
area_field_name = "PLOCHA_M2"
osa_field_name = "HLAVNI_OSA"
group_id_field = "GROUP_ID"  # Přidaná definice
DMR_raster = "DMR_RASTER"
# Nastavení scratch workspace pro dočasné soubory
arcpy.env.scratchWorkspace = arcpy.env.scratchGDB

# Definice názvů polí pro tvarové metriky
rectangularity_field = "RECTANGULARITY"
ellipticity_field = "ELLIPTICITY"
perimeter_area_ratio_field = "PERIM_AREA_RAT"
circularity_field = "CIRCULARITY"
convexity_field = "CONVEXITY"
elongation_field = "ELONGATION"
major_axis_field = "HLAVNI_OSA"
minor_axis_field = "VEDLEJSI_OSA"
shape_type_field = "TVAR_TUNE"


# Funkce pro klasifikaci tvaru tůně
def classify_shape(rectangularity, ellipticity, circularity, convexity, elongation):
    if rectangularity > 0.95 and ellipticity < 1.1 and circularity > 0.75 and convexity > 0.90:
        return "Čtverec"
    elif rectangularity > 0.4 and 1.1 <= ellipticity <= 3 and circularity < 0.5 and convexity > 0.6:
        return "Obdélník"
    elif rectangularity > 0.4 and ellipticity > 3 and circularity < 0.5 and convexity > 0.6:
        return "Protáhlý obdélník"
    elif circularity > 0.85 and convexity > 0.95 and elongation < 0.3 and ellipticity < 1.1:
        return "Kruh"
    elif circularity > 0.8 and convexity > 0.95 and elongation < 0.3 and ellipticity < 1.5:
        return "Elipsa"
    elif circularity > 0.6 and convexity > 0.95 and elongation > 0.3 and ellipticity >= 2.0:
        return "Protáhlá elipsa"
    elif ellipticity > 2.5 and circularity < 0.3 and elongation > 0.7 and convexity < 0.9:
        return "Protáhlý tvar"
    elif 1.0 < ellipticity < 1.5 and circularity < 0.35 and convexity < 0.6 and elongation < 0.35:
        return "Tvar písmene C"
    elif circularity < 0.4 and convexity < 0.8:
        return "Členitý tvar"
    elif ellipticity > 2.0 and 0.5 < circularity < 0.8 and convexity > 0.9:
        return "Protáhlá zaoblená tůň"
    else:
        return "Nepravidelný tvar"



try:
    # *** INICIALIZACE PROMENNÝCH HNED NA ZAČÁTKU TRY BLOKU ***
    output_tune_pred_sloucenim = "potencialni_tune_pred_sloucenim"
    output_final_tune_before_dissolve = None
    output_tune_pred_sloucenim_full = None
    rivers_buffer_distance = 10  # Vzdálenost bufferu kolem řek v metrech

    # Kontrola existence vstupního souboru
    if not arcpy.Exists(input_features):
        arcpy.AddError(f"Vstupní data '{input_features}' neexistují.")
        raise arcpy.ExecuteError

    # Získání popisu vstupních dat
    desc = arcpy.Describe(input_features)
    
    if not arcpy.Exists(input_features):
        arcpy.AddError(...)
        raise arcpy.ExecuteError

    # Získání popisu vstupních dat
    desc = arcpy.Describe(input_features)

    # --- DETEKCE A FILTRACE LAS DATASETU NA GROUND BODY ---
    if desc.dataType.lower() == "lasdataset" or input_features.lower().endswith(".las"):
        arcpy.AddMessage("Vstup je LAS dataset – filtrovat jen ground body...")

        # výstupní multipoint pro ground body
        filtered_fc = os.path.join(arcpy.env.scratchGDB, "filtered_ground_points")

        # konverze LAS → Multipoint (jen klasifikace 2 = ground)
        arcpy.ddd.LASToMultipoint(
            in_las_dataset=input_features,
            out_feature_class=filtered_fc,
            average_point_spacing=cell_size_float,
            class_code="2"
        )
        arcpy.AddMessage(f"Ground body vyfiltrovány do: {filtered_fc}")

        # teď budeme dál pracovat s touto feature class
        input_features = filtered_fc

        # ošetření pole pro nadmořskou výšku
        # pokud původní value_field neexistuje v nové vrstvě, vytvoříme pole 'Elevation'
        if not arcpy.ListFields(input_features, value_field):
            auto_field = "Elevation"
            arcpy.AddMessage(f"Pole '{value_field}' neexistuje – vytvořím pole '{auto_field}'.")
            arcpy.AddField_management(input_features, auto_field, "DOUBLE")
            arcpy.CalculateField_management(input_features, auto_field, "!SHAPE.Z!", "PYTHON3")
            value_field = auto_field
            arcpy.AddMessage(f"Pole '{auto_field}' připraveno jako nadmořská výška.")
    else:
        arcpy.AddMessage("Vstup je feature class – filtraci LAS datasetu přeskočím.")
    # --- KONEC FILTRACE LAS DATASETU ---


    # Získání pracovního prostoru (geodatabáze), kde jsou uložena vstupní data
    input_workspace = desc.catalogPath.split('\\')[:-1]
    input_workspace = '\\'.join(input_workspace)
    arcpy.AddMessage(f"Pracovní prostor vstupu: {input_workspace}")

    # Výsledná cesta k rastru
    output_raster_expected = os.path.join(input_workspace, DMR_raster)
    output_raster = output_raster_expected # Při úspěchu se tato hodnota přepíše

    # Kontrola, zda výstupní rastr již existuje
    if arcpy.Exists(output_raster):
        arcpy.AddWarning(f"Výstupní rastr '{DMR_raster}' již existuje a bude přepsán.")

    # Zpracování velikosti buňky s ohledem na lokalizaci
    try:
        cell_size_float = float(cell_size)
    except ValueError:
        cell_size = cell_size.replace(',', '.')
        try:
            cell_size_float = float(cell_size)
        except ValueError:
            arcpy.AddError(f"Nelze převést '{cell_size}' na číslo. Použijte formát jako '0.5' nebo '0,5'.")
            raise arcpy.ExecuteError

    arcpy.AddMessage(f"Vstupní feature class: {input_features}")
    arcpy.AddMessage(f"Použité hodnoty pole: {value_field}")
    arcpy.AddMessage(f"Velikost buňky: {cell_size_float}")
    arcpy.AddMessage(f"Výstupní rastr bude uložen do: {input_workspace} jako {DMR_raster}")

    arcpy.env.workspace = input_workspace
    arcpy.env.outputCoordinateSystem = arcpy.Describe(input_features).spatialReference

    arcpy.AddMessage("Konvertuji bodové mračno na rastr...")

    # Převod bodového mračna na rastr pomocí funkce PointToRaster
    try:
        arcpy.conversion.PointToRaster(
            in_features=input_features,
            value_field=value_field,
            out_rasterdataset=output_raster,
            cell_assignment="MEAN",
            cellsize=cell_size_float
        )
        arcpy.AddMessage("Konverze byla úspěšně dokončena.")

        # Nastavení vytvořeného rastru jako snap raster pro další operace
        arcpy.env.snapRaster = output_raster
    

        # --- PŘEVOD NoData HODNOT NA 1 POMOCÍ NÁSTROJE Con ---
        in_raster = output_raster
        out_raster_no_data = "nodata_raster"  # Jednodušší název
        output_raster_no_data_full = os.path.join(input_workspace, out_raster_no_data)

        arcpy.AddMessage("Převádím NoData pixely na hodnotu 1 pomocí nástroje Con...")

        # Vytvoříme podmínku: pokud je hodnota NoData, nastavíme 1, jinak 0
        out_con = arcpy.sa.Con(arcpy.sa.IsNull(in_raster), 1, 0)
        out_con.save(output_raster_no_data_full)
        output_raster = output_raster_no_data_full # Pro další kroky budeme pracovat s tímto rastrem
        arcpy.AddMessage("Převod NoData pixelů pomocí nástroje Con dokončen.")
        # --- KONEC PŘEVODU NoData HODNOT ---

        # --- SESKUPOVÁNÍ SOUVISLÝCH OBLASTÍ ---
        in_raster_grouped = output_raster  # Používáme rastr s hodnotou 1 pro NoData
        out_region_group_name = "seskupeny_raster"
        out_region_group_full = os.path.join(arcpy.env.scratchFolder, "grp_rst.tif")

        arcpy.AddMessage("Seskupuji souvislé oblasti NoData...")
        grouped_raster = arcpy.sa.RegionGroup(in_raster_grouped, "EIGHT", "WITHIN")
        grouped_raster.save(out_region_group_full)
        output_raster_grouped = out_region_group_full # Pro další kroky
        arcpy.AddMessage("Seskupení souvislých oblastí dokončeno.")
        # --- KONEC SESKUPOVÁNÍ ---

        # --- PŘEVOD RASTROVÝCH OBLASTÍ NA POLYGONY ---
        in_raster_for_polygon = output_raster_grouped
        out_polygon_features = "tune_polygony"
        out_polygon_features_full = os.path.join(input_workspace, "tune_polygony") # Jednoduchý název pro GDB

        arcpy.env.workspace = input_workspace # Explicitně nastavíme pracovní prostředí

        arcpy.conversion.RasterToPolygon(in_raster_for_polygon, out_polygon_features_full, "SIMPLIFY", "VALUE", False, None) # Explicitní parametry
        output_polygon_layer = out_polygon_features_full
        arcpy.AddMessage("Převod na polygony dokončen.")
        # --- KONEC PŘEVODU NA POLYGONY ---

        # --- VYHLAZENÍ POLYGONŮ ---
        smoothed_features = "hladke_polygony"
        smoothing_algorithm = "PAEK" # Nebo "BEZIER_INTERPOLATION"
        smoothing_tolerance = peak_al # Nastavte vhodnou toleranci

        arcpy.AddMessage(f"Vyhlazuji polygony pomocí algoritmu '{smoothing_algorithm}' s tolerancí '{smoothing_tolerance}'...")
        smoothed_features_full = os.path.join(input_workspace, smoothed_features)
        arcpy.cartography.SmoothPolygon(output_polygon_layer, smoothed_features_full, smoothing_algorithm, smoothing_tolerance, "NO_FIXED")
        output_polygon_layer = smoothed_features_full # Pro další kroky budeme pracovat s vyhlazenými polygony
        out_polygon_layer_with_area = output_polygon_layer # Ujistíme se, že odkazuje na správnou vrstvu
        arcpy.AddMessage("Vyhlazení polygonů dokončeno.")
        # --- KONEC VYHLAZENÍ ---

        # --- VÝPOČET PLOCHY POLYGONŮ (nyní pro vyhlazené polygony) ---
        arcpy.AddMessage("Přidávám pole pro plochu...")
        arcpy.AddField_management(out_polygon_layer_with_area, area_field_name, "DOUBLE")
        arcpy.AddMessage("Pole pro plochu přidáno.")

        arcpy.AddMessage("Počítám plochu polygonů...")
        with arcpy.da.UpdateCursor(out_polygon_layer_with_area, [area_field_name, "SHAPE@AREA"]) as cursor:
            for row in cursor:
                row[0] = row[1]
                cursor.updateRow(row)
        arcpy.AddMessage("Plocha polygonů spočítána.")
        # --- KONEC VÝPOČTU PLOCHY ---

        # --- FILTROVÁNÍ POLYGONŮ PODLE PLOCHY (PEVNÁ HODNOTA - MIN. A MAX.) ---
        arcpy.AddMessage(f"Filtruji polygony s plochou mezi {minimal_plocha_tuhne} a {maximalni_plocha_tuhne} m2...")

        # Vytvoření dočasné vrstvy
        temp_layer = "temp_polygons"
        arcpy.MakeFeatureLayer_management(out_polygon_layer_with_area, temp_layer)

        # Výběr polygonů podle atributu
        query = f"{area_field_name} ≥ {minimal_plocha_tuhne} AND {area_field_name} ≤ {maximalni_plocha_tuhne}"
        arcpy.SelectLayerByAttribute_management(temp_layer, "NEW_SELECTION", query)

        # Uložení vybraných prvků do nové feature class
        output_tune_pred_sloucenim_full = os.path.join(input_workspace, output_tune_pred_sloucenim)
        arcpy.CopyFeatures_management(temp_layer, output_tune_pred_sloucenim_full)
        output_final_tune_before_dissolve = output_tune_pred_sloucenim_full

        arcpy.AddMessage(f"Vytvořena vrstva potenciálních tůní před sloučením: {output_final_tune_before_dissolve}")
        # --- KONEC FILTROVÁNÍ ---