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Í ---