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_tuny_pred_sloucenim = "potencialni_tune_pred_sloucenim" output_final_tuny_before_dissolve = None output_tuny_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_tuny_pred_sloucenim_full = os.path.join(input_workspace, output_tuny_pred_sloucenim) arcpy.CopyFeatures_management(temp_layer, output_tuny_pred_sloucenim_full) output_final_tuny_before_dissolve = output_tuny_pred_sloucenim_full arcpy.AddMessage(f"Vytvořena vrstva potenciálních tůní před sloučením: {output_final_tuny_before_dissolve}") # --- KONEC FILTROVÁNÍ --- # --- PŘIDÁNÍ UNIKÁTNÍHO ID PRO KAŽDÝ POLYGON --- arcpy.AddMessage("Přidávám unikátní ID pro každý polygon...") orig_id_field = "ORIG_ID" arcpy.AddField_management(output_final_tuny_before_dissolve, orig_id_field, "LONG") with arcpy.da.UpdateCursor(output_final_tuny_before_dissolve, [orig_id_field, "OID@"]) as cursor: for row in cursor: row[0] = row[1] # Použijeme OID jako unikátní ID cursor.updateRow(row) arcpy.AddMessage("Unikátní ID přidáno.") # --- KONEC PŘIDÁNÍ UNIKÁTNÍHO ID --- # --- VYTVOŘENÍ BUFFERU PRO IDENTIFIKACI PŘEKRÝVAJÍCÍCH SE TŮNÍ --- arcpy.AddMessage("Vytvářím buffer kolem každé tůně...") buffered_polygons = os.path.join(input_workspace, "tune_buffer") # Vytvoření bufferu kolem každé tůně s konstantní hodnotou 4 metry arcpy.Buffer_analysis( in_features=output_final_tuny_before_dissolve, out_feature_class=buffered_polygons, buffer_distance_or_field="2 Meters", # Změněno z 5 na 4 metry line_side="FULL", line_end_type="ROUND", dissolve_option="NONE" # Zachování jednotlivých záznamů ) arcpy.AddMessage("Buffery vytvořeny.") # --- KONEC VYTVOŘENÍ BUFFERU --- # Přidání pole GROUP_ID do buffered_polygons arcpy.AddMessage("Přidávám pole GROUP_ID do vrstvy bufferů...") arcpy.AddField_management(buffered_polygons, "GROUP_ID", "LONG") # Inicializace GROUP_ID na -1 with arcpy.da.UpdateCursor(buffered_polygons, ["GROUP_ID"]) as cursor: for row in cursor: row[0] = -1 cursor.updateRow(row) # --- ITERATIVNÍ IDENTIFIKACE PŘEKRÝVAJÍCÍCH SE BUFFERŮ --- arcpy.AddMessage("Začínám iterativní slučování překrývajících se tůní...") # Funkce pro provedení jednoho cyklu slučování def perform_merge_iteration(input_features, workspace, iteration_number): arcpy.AddMessage(f"Začínám {iteration_number}. cyklus slučování...") # Vytvoření dočasné vrstvy pro buffery buffer_layer = f"buffer_layer_{iteration_number}" if arcpy.Exists(buffer_layer): arcpy.Delete_management(buffer_layer) arcpy.MakeFeatureLayer_management(input_features, buffer_layer) # Připravení pole pro skupiny group_id_field = "GROUP_ID" field_names = [f.name for f in arcpy.ListFields(input_features)] if group_id_field not in field_names: arcpy.AddField_management(input_features, group_id_field, "LONG") # Nastavení výchozí hodnoty GROUP_ID na -1 (samostatné tůně) with arcpy.da.UpdateCursor(input_features, [group_id_field]) as cursor: for row in cursor: row[0] = -1 cursor.updateRow(row) # Inicializace čítače skupin group_counter = 1 changes_made = False processed_ids = [] # Procházíme všechny buffery a hledáme překryvy with arcpy.da.SearchCursor(input_features, ["OID@", "SHAPE@", group_id_field]) as cursor: for row in cursor: oid = row[0] geom = row[1] current_group = row[2] if oid in processed_ids: continue # Vytvoříme výběr překrývajících se bufferů arcpy.SelectLayerByLocation_management( buffer_layer, "INTERSECT", geom, selection_type="NEW_SELECTION" ) # Získáme seznam ID a aktuálních skupin vybraných bufferů overlapping_data = [] with arcpy.da.SearchCursor(buffer_layer, ["OID@", group_id_field]) as overlap_cursor: for overlap_row in overlap_cursor: overlapping_data.append((overlap_row[0], overlap_row[1])) # Máme překryv s více než jedním bufferem if len(overlapping_data) > 1: changes_made = True # Zjistíme, jestli některé z překrývajících se polygonů již mají přiřazenou skupinu assigned_group_ids = [gid for _, gid in overlapping_data if gid > 0] # Určíme jakou skupinu použít: existující nebo novou target_group_id = min(assigned_group_ids) if assigned_group_ids else group_counter if not assigned_group_ids: group_counter += 1 # Přiřadíme vybranou skupinu všem překrývajícím se bufferům for overlap_id, overlap_group in overlapping_data: processed_ids.append(overlap_id) # Aktualizace v aktuální vrstvě bufferů where_clause = f"OBJECTID = {overlap_id}" with arcpy.da.UpdateCursor(buffer_layer, [group_id_field], where_clause) as update_cursor: for update_row in update_cursor: update_row[0] = target_group_id update_cursor.updateRow(update_row) if changes_made: # Vytvoření výstupní vrstvy pro tento cyklus output_merged = os.path.join(workspace, f"merged_polygons_{iteration_number}") # Nejprve vybereme a rozpustíme polygony se skupinou arcpy.SelectLayerByAttribute_management(buffer_layer, "NEW_SELECTION", f"{group_id_field} > 0") temp_dissolved = os.path.join(workspace, f"temp_dissolved_{iteration_number}") arcpy.Dissolve_management( in_features=buffer_layer, out_feature_class=temp_dissolved, dissolve_field=[group_id_field], multi_part="MULTI_PART" ) # Pak vybereme samostatné polygony arcpy.SelectLayerByAttribute_management(buffer_layer, "NEW_SELECTION", f"{group_id_field} = -1") temp_standalone = os.path.join(workspace, f"temp_standalone_{iteration_number}") arcpy.CopyFeatures_management(buffer_layer, temp_standalone) # Sloučíme rozpuštěné a samostatné do výstupní vrstvy arcpy.Merge_management([temp_dissolved, temp_standalone], output_merged) # Vyčistíme for temp_fc in [temp_dissolved, temp_standalone]: if arcpy.Exists(temp_fc): arcpy.Delete_management(temp_fc) return output_merged, changes_made return input_features, changes_made # První cyklus slučování merged_features_1, changes_1 = perform_merge_iteration(buffered_polygons, input_workspace, 1) arcpy.AddMessage("První cyklus slučování dokončen.") # Druhý cyklus slučování if changes_1: merged_features_2, changes_2 = perform_merge_iteration(merged_features_1, input_workspace, 2) final_merged_features = merged_features_2 arcpy.AddMessage("Druhý cyklus slučování dokončen.") else: final_merged_features = merged_features_1 arcpy.AddMessage("Druhý cyklus slučování nebyl potřeba.") # Vytvoříme prostorové spojení mezi původními polygony a finálními buffery spatial_join = "spatial_join" spatial_join_full = os.path.join(input_workspace, spatial_join) # Přidání pole GROUP_ID do původních polygonů před spojením if not arcpy.ListFields(output_final_tuny_before_dissolve, group_id_field): arcpy.AddField_management(output_final_tuny_before_dissolve, group_id_field, "LONG") arcpy.SpatialJoin_analysis( target_features=output_final_tuny_before_dissolve, join_features=final_merged_features, out_feature_class=spatial_join_full, join_operation="JOIN_ONE_TO_ONE", join_type="KEEP_ALL", match_option="INTERSECT" # Změněno z HAVE_THEIR_CENTER_IN na INTERSECT ) # Aktualizace GROUP_ID v původních polygonech podle prostorového spojení arcpy.MakeFeatureLayer_management(output_final_tuny_before_dissolve, "original_layer") # Kontrola, zda pole GROUP_ID existuje ve výsledku spojení field_names = [f.name for f in arcpy.ListFields(spatial_join_full)] join_field_name = f"JOIN_{group_id_field}" if join_field_name not in field_names: if group_id_field in field_names: join_field_name = group_id_field else: arcpy.AddError(f"Nelze najít pole '{join_field_name}' ani '{group_id_field}' ve výsledku prostorového spojení.") arcpy.AddMessage("Dostupná pole: " + ", ".join(field_names)) raise arcpy.ExecuteError arcpy.AddMessage(f"Používám pole {join_field_name} pro přenos GROUP_ID.") with arcpy.da.SearchCursor(spatial_join_full, ["ORIG_ID", join_field_name]) as cursor: for row in cursor: orig_id = row[0] group_id = row[1] where_clause = f"ORIG_ID = {orig_id}" with arcpy.da.UpdateCursor("original_layer", [group_id_field], where_clause) as update_cursor: for update_row in update_cursor: update_row[0] = group_id update_cursor.updateRow(update_row) arcpy.AddMessage("Iterativní slučování dokončeno.") # --- KONEC ITERATIVNÍHO SLUČOVÁNÍ --- # --- VYTVOŘENÍ VÝSLEDNÉ VRSTVY S PŘEKRÝVAJÍCÍMI SE A SAMOSTATNÝMI TŮNĚMI --- arcpy.AddMessage("Vytvářím výslednou vrstvu...") # Sloučení překrývajících se bufferů podle GROUP_ID overlapping_pools = os.path.join(input_workspace, "prekryvajici_se_tune") arcpy.MakeFeatureLayer_management(buffered_polygons, "tune_buffer_layer") # Výběr tůní, které mají překryv (GROUP_ID > 0) arcpy.SelectLayerByAttribute_management("tune_buffer_layer", "NEW_SELECTION", f"{group_id_field} > 0") # Sloučení překrývajících se tůní podle GROUP_ID arcpy.Dissolve_management( in_features="tune_buffer_layer", out_feature_class=overlapping_pools, dissolve_field=[group_id_field], multi_part="MULTI_PART" ) # Výběr samostatných tůní (GROUP_ID = -1) arcpy.SelectLayerByAttribute_management("tune_buffer_layer", "NEW_SELECTION", f"{group_id_field} = -1") # Uložení samostatných tůní do dočasné vrstvy standalone_pools = os.path.join(input_workspace, "samostatne_tune") arcpy.CopyFeatures_management("tune_buffer_layer", standalone_pools) # Sloučení obou typů tůní do jedné výsledné vrstvy final_output = "vysledne_tune" final_output_full = os.path.join(input_workspace, final_output) arcpy.Merge_management([overlapping_pools, standalone_pools], final_output_full) # --- ZPRACOVÁNÍ VOLITELNÝCH VRSTEV (ŘEKY A BUDOVY) --- if rivers_layer and arcpy.Exists(rivers_layer): arcpy.AddMessage("Zpracovávám vrstvu řek...") # Vytvoření bufferu kolem řek rivers_buffer = os.path.join(input_workspace, "rivers_buffer") arcpy.Buffer_analysis( in_features=rivers_layer, out_feature_class=rivers_buffer, buffer_distance_or_field="3 Meters", line_side="FULL", line_end_type="ROUND", dissolve_option="ALL" ) # Odstranění tůní, které se překrývají s řekami arcpy.MakeFeatureLayer_management(final_output_full, "final_tune_layer") # Počítání tůní před odstraněním result = arcpy.GetCount_management("final_tune_layer") count_before = int(result[0]) # Výběr tůní, které jsou blízko řek arcpy.SelectLayerByLocation_management("final_tune_layer", "WITHIN_A_DISTANCE", rivers_buffer, "3 Meters") # Počítání vybraných tůní result = arcpy.GetCount_management("final_tune_layer") count_selected = int(result[0]) # Odstranění vybraných tůní arcpy.DeleteFeatures_management("final_tune_layer") # Počítání tůní po odstranění result = arcpy.GetCount_management("final_tune_layer") count_after = int(result[0]) arcpy.AddMessage(f"Bylo odstraněno {count_before - count_after} tůní, které procházely řekami nebo byly příliš blízko řek.") if buildings_layer and arcpy.Exists(buildings_layer): arcpy.AddMessage("Zpracovávám vrstvu budov...") # Vytvoření bufferu kolem budov buildings_buffer = os.path.join(input_workspace, "buildings_buffer") arcpy.Buffer_analysis( in_features=buildings_layer, out_feature_class=buildings_buffer, buffer_distance_or_field="3 Meters", line_side="FULL", line_end_type="ROUND", dissolve_option="ALL" ) # Odstranění tůní, které se překrývají s budovami arcpy.MakeFeatureLayer_management(final_output_full, "final_tune_layer") # Počítání tůní před odstraněním result = arcpy.GetCount_management("final_tune_layer") count_before = int(result[0]) # Výběr tůní, které jsou blízko budov arcpy.SelectLayerByLocation_management("final_tune_layer", "WITHIN_A_DISTANCE", buildings_buffer, "3 Meters") # Počítání vybraných tůní result = arcpy.GetCount_management("final_tune_layer") count_selected = int(result[0]) # Odstranění vybraných tůní arcpy.DeleteFeatures_management("final_tune_layer") # Počítání tůní po odstranění result = arcpy.GetCount_management("final_tune_layer") count_after = int(result[0]) arcpy.AddMessage(f"Bylo odstraněno {count_before - count_after} tůní, které procházely budovami nebo byly příliš blízko budov.") # --- PŘEPOČET PLOCHY A TVAROVÝCH METRIK PRO VÝSLEDNOU VRSTVU --- arcpy.AddMessage("Přidávám a počítám tvarové metriky...") try: # Nejprve přidáme všechna potřebná pole arcpy.AddField_management(final_output_full, rectangularity_field, "DOUBLE", "", "", "", "Rectangularity") arcpy.AddField_management(final_output_full, ellipticity_field, "DOUBLE", "", "", "", "Ellipticity") arcpy.AddField_management(final_output_full, perimeter_area_ratio_field, "DOUBLE", "", "", "", "Perimeter/Area Ratio") arcpy.AddField_management(final_output_full, circularity_field, "DOUBLE", "", "", "", "Circularity") arcpy.AddField_management(final_output_full, convexity_field, "DOUBLE", "", "", "", "Convexity") arcpy.AddField_management(final_output_full, elongation_field, "DOUBLE", "", "", "", "Elongation") arcpy.AddField_management(final_output_full, major_axis_field, "DOUBLE", "", "", "", "Hlavní osa") arcpy.AddField_management(final_output_full, minor_axis_field, "DOUBLE", "", "", "", "Vedlejší osa") arcpy.AddField_management(final_output_full, shape_type_field, "TEXT", "", "", 50, "Tvar tůně") # Přidáno pole pro tvar arcpy.AddMessage("Pole pro tvarové metriky byla úspěšně přidána.") # Kontrola, zda byla pole skutečně přidána existing_fields = [f.name for f in arcpy.ListFields(final_output_full)] arcpy.AddMessage(f"Existující pole: {', '.join(existing_fields)}") # Výpočet metrik pro každý polygon arcpy.AddMessage("Počítám tvarové metriky pro každý polygon...") # Seznam všech polí pro výpočet fields = ["SHAPE@", rectangularity_field, ellipticity_field, perimeter_area_ratio_field, circularity_field, convexity_field, elongation_field, major_axis_field, minor_axis_field, shape_type_field] with arcpy.da.UpdateCursor(final_output_full, fields) as cursor: for row in cursor: # Získání geometrie geometry = row[0] # Základní měření area = geometry.area perimeter = geometry.length extent = geometry.extent width = extent.width height = extent.height # Výpočet metrik try: # Rectangularity rectangularity = area / (width * height) if (width * height) != 0 else 0 row[1] = rectangularity # Ellipticity major_axis = max(width, height) minor_axis = min(width, height) ellipticity = major_axis / minor_axis if minor_axis != 0 else 0 row[2] = ellipticity # Perimeter/Area Ratio row[3] = perimeter / area if area != 0 else 0 # Circularity circularity = (4 * math.pi * area) / (perimeter**2) if perimeter != 0 else 0 row[4] = circularity # Convexity convex_hull = geometry.convexHull() convexity = area / convex_hull.area if convex_hull.area != 0 else 0 row[5] = convexity # Elongation elongation = 1 - (minor_axis / major_axis) if major_axis != 0 else 0 row[6] = elongation # Hlavní osa row[7] = major_axis # Vedlejší osa row[8] = minor_axis # Klasifikace tvaru row[9] = classify_shape(rectangularity, ellipticity, circularity, convexity, elongation) cursor.updateRow(row) except Exception as e: arcpy.AddWarning(f"Chyba při výpočtu metrik pro polygon: {str(e)}") continue arcpy.AddMessage("Výpočet tvarových metrik a klasifikace tvarů dokončeny.") except Exception as e: arcpy.AddError(f"Chyba při přidávání nebo výpočtu tvarových metrik: {str(e)}") raise # --- VYTVOŘENÍ NOVÉ VRSTVY S FILTRACÍ PODLE KONVEXITY A HLAVNÍ OSY --- arcpy.AddMessage("Vytvářím novou vrstvu filtrovanou podle konvexity a hlavní osy...") # Vytvoření dočasné vrstvy filtered_output = "vysledne_tune_konvexita" filtered_output_full = os.path.join(input_workspace, filtered_output) # Nejprve zkopírujeme původní vrstvu arcpy.CopyFeatures_management(final_output_full, filtered_output_full) # Vytvoření layer file pro výběr arcpy.MakeFeatureLayer_management(filtered_output_full, "konvexita_layer") # Výběr polygonů s konvexitou >= 0.6 # odstranění všech tůní s hlavní osou > 80 where_clause = f"{convexity_field} >= 0.6 AND {major_axis_field} <= 80" arcpy.SelectLayerByAttribute_management("konvexita_layer", "NEW_SELECTION", where_clause) # Vytvoření nové feature class pouze s vybranými prvky final_filtered_output = "vysledne_tune_konvexita_filtered" final_filtered_output_full = os.path.join(input_workspace, final_filtered_output) arcpy.CopyFeatures_management("konvexita_layer", final_filtered_output_full) # Přidání informace o počtu odfiltrovaných tůní result = arcpy.GetCount_management(final_output_full) total_count = int(result[0]) result = arcpy.GetCount_management(final_filtered_output_full) filtered_count = int(result[0]) removed_count = total_count - filtered_count arcpy.AddMessage(f"Celkový počet tůní před filtrací: {total_count}") arcpy.AddMessage(f"Počet tůní po filtraci (konvexita >= 0.6 a hlavní osa <= 80): {filtered_count}") arcpy.AddMessage(f"Počet odstraněných tůní: {removed_count}") # --- VYTVOŘENÍ BUFFERU PRO OŘÍZNUTÍ BODŮ --- arcpy.AddMessage("Vytvářím buffer pro oříznutí bodů...") # Odstranění nepotřebných polí arcpy.AddMessage("Odstraňuji nepotřebná pole...") fields_to_delete = ["Id", "gridcode", "InPoly_FID", "PLOCHA_M2", "ORIG_ID", "BUFF_DIST", "ORIG_FID"] # Získání seznamu existujících polí existing_fields = [f.name for f in arcpy.ListFields(final_output_full)] # Odstranění pouze těch polí, která existují fields_to_actually_delete = [field for field in fields_to_delete if field in existing_fields] if fields_to_actually_delete: arcpy.DeleteField_management(final_output_full, fields_to_actually_delete) arcpy.AddMessage(f"Odstraněna pole: {', '.join(fields_to_actually_delete)}") else: arcpy.AddMessage("Žádná z určených polí nebyla nalezena k odstranění.") buffer_pro_orez = "buffer_pro_orez" buffer_pro_orez_full = os.path.join(input_workspace, buffer_pro_orez) arcpy.Buffer_analysis( in_features=final_filtered_output_full, out_feature_class=buffer_pro_orez_full, buffer_distance_or_field= "3 Meters", # Konstantní hodnota line_side="FULL", line_end_type="ROUND", dissolve_option="NONE" # Sloučíme všechny buffery do jednoho ) arcpy.AddMessage("Buffer pro oříznutí bodů vytvořen.") # --- KONEC VYTVOŘENÍ BUFFERU --- # --- OŘEZÁNÍ BODOVÉHO MRAČNA POMOCÍ JEDNOTLIVÝCH BUFFERŮ --- arcpy.AddMessage("Začínám ořezávat bodové mračno pomocí jednotlivých bufferů...") # Vytvoření seznamu pro ukládání ORIG_FID orig_fid_list = [] # Vytvoření FeatureLayer z vstupního bodového mračna pro výběr in_points_layer = "vstupni_body_pro_orez" arcpy.MakeFeatureLayer_management(input_features, in_points_layer) # Vytvoření FeatureLayer z bufferů pro iteraci buffer_layer = "buffery_pro_iteraci" arcpy.MakeFeatureLayer_management(buffer_pro_orez_full, buffer_layer) # Získání počtu bufferů pro informaci result = arcpy.GetCount_management(buffer_layer) num_buffers = int(result[0]) arcpy.AddMessage(f"Počet bufferů k ořezání: {num_buffers}") # Proměnná pro počítání vytvořených výstupních vrstev output_counter = 1 # Iterace přes každý buffer with arcpy.da.SearchCursor(buffer_layer, ["OID@", "SHAPE@"]) as cursor: for row in cursor: buffer_oid = row[0] buffer_geometry = row[1] # Provedení prostorového výběru bodů uvnitř aktuálního bufferu arcpy.SelectLayerByLocation_management(in_points_layer, "WITHIN", buffer_geometry, None, "NEW_SELECTION") # Vytvoření názvu pro výstupní vrstvu output_points_name = f"tune_body_{output_counter}" output_points_full = os.path.join(input_workspace, output_points_name) # Vytvoření nové feature class arcpy.CreateFeatureclass_management( out_path=os.path.dirname(output_points_full), out_name=os.path.basename(output_points_full), geometry_type="POINT", template=in_points_layer, has_m="ENABLED", has_z="ENABLED" ) # Přidání pole "TVAR_TUNE", pokud neexistuje if not any(f.name.lower() == "tvar_tune" for f in arcpy.ListFields(output_points_full)): arcpy.AddField_management(output_points_full, "TVAR_TUNE", "TEXT", field_length=50) # Uložení ORIG_FID pro aktuální buffer with arcpy.da.SearchCursor(buffer_pro_orez_full, ["ORIG_FID"], f"OBJECTID = {buffer_oid}") as search_cursor_buffer: for search_row_buffer in search_cursor_buffer: orig_fid = search_row_buffer[0] orig_fid_list.append(orig_fid) break # Získání názvů polí X, Y a Elevation x_field_name = None y_field_name = None elevation_field_name = None for field in arcpy.ListFields(in_points_layer): if field.type.lower() in ["double", "single", "integer", "smallinteger", "oid"]: if field.name.lower() in ["x", "point_x", "shape.pointx"]: x_field_name = field.name elif field.name.lower() in ["y", "point_y", "shape.pointy"]: y_field_name = field.name elif field.name.lower() in ["elevation", "z", "point_z", "shape.pointz"]: elevation_field_name = field.name # Vložení dat fields_to_insert = ["SHAPE@"] fields_to_search = ["SHAPE@"] if x_field_name: fields_to_insert.append("X") fields_to_search.append(x_field_name) if y_field_name: fields_to_insert.append("Y") fields_to_search.append(y_field_name) if elevation_field_name: fields_to_insert.append("Elevation") fields_to_search.append(elevation_field_name) with arcpy.da.InsertCursor(output_points_full, fields_to_insert) as insert_cursor: with arcpy.da.SearchCursor(in_points_layer, fields_to_search) as search_cursor: for search_row in search_cursor: insert_cursor.insertRow(search_row) arcpy.AddMessage(f"Vytvořena vrstva bodů pro buffer {buffer_oid}: {output_points_full}") # Přenos "TVAR_TUNE" z bufferu with arcpy.da.SearchCursor(buffer_pro_orez_full, ["OID@", "TVAR_TUNE"], f"OBJECTID = {buffer_oid}") as search_cursor_buffer: for search_row_buffer in search_cursor_buffer: tvar_tune_value = search_row_buffer[1] break with arcpy.da.UpdateCursor(output_points_full, ["TVAR_TUNE"]) as update_cursor_points: for update_row_points in update_cursor_points: update_row_points[0] = tvar_tune_value update_cursor_points.updateRow(update_row_points) arcpy.SelectLayerByAttribute_management(in_points_layer, "CLEAR_SELECTION") output_counter += 1 if arcpy.Exists(in_points_layer): arcpy.Delete_management(in_points_layer) if arcpy.Exists(buffer_layer): arcpy.Delete_management(buffer_layer) arcpy.AddMessage("Ořezání bodového mračna pomocí bufferů dokončeno. Vytvořeno " + str(output_counter - 1) + " vrstev bodových mračen.") # --- KONEC OŘEZÁNÍ BODOVÉHO MRAČNA --- # Odstranění děr v polygonech arcpy.AddMessage("Odstraňuji všechny díry v polygonech...") # Vytvoření dočasné vrstvy bez děr no_holes_output = "vysledne_tune_no_holes" no_holes_output_full = os.path.join(input_workspace, no_holes_output) # Použití nástroje EliminatePolygonPart pro odstranění všech děr arcpy.EliminatePolygonPart_management( in_features=final_output_full, out_feature_class=no_holes_output_full, condition="PERCENT", part_area="0", part_area_percent="0", part_option="CONTAINED_ONLY" ) # Aktualizace proměnné pro další zpracování final_output_full = no_holes_output_full # Vyčištění pouze dočasných vrstev temp_layers = ["grp_rst", "temp_dissolved", "temp_standalone", "spatial_join", "konvexita_layer"] # Odstranění polygonů s plochou větší než 5000 m2 arcpy.AddMessage("Odstraňuji polygony s plochou větší než 5000 m2...") # Vytvoření dočasné vrstvy pro výběr arcpy.MakeFeatureLayer_management(final_output_full, "area_filter_layer") # Výběr polygonů s plochou menší nebo rovnou 5000 m2 where_clause = f"Shape_Area <= 5000" arcpy.SelectLayerByAttribute_management("area_filter_layer", "NEW_SELECTION", where_clause) # Vytvoření nové feature class pouze s vybranými prvky final_filtered_by_area = "vysledne_tune_filtered" final_filtered_by_area_full = os.path.join(input_workspace, final_filtered_by_area) arcpy.CopyFeatures_management("area_filter_layer", final_filtered_by_area_full) # Přidání informace o počtu odstraněných polygonů result = arcpy.GetCount_management(final_filtered_output_full) total_count = int(result[0]) result = arcpy.GetCount_management(final_filtered_output_full) filtered_count = int(result[0]) removed_count = total_count - filtered_count arcpy.AddMessage(f"Celkový počet tůní před filtrací dle plochy: {total_count}") arcpy.AddMessage(f"Počet tůní po filtraci (plocha <= 5000 m2): {filtered_count}") arcpy.AddMessage(f"Počet odstraněných tůní: {removed_count}") # --- EBK INTERPOLACE PRO TŮNĚ S PEVNOU PARAMETRIZACÍ --- arcpy.AddMessage("Provádím EBK interpolaci všech tůní s jednotnými parametry...") arcpy.env.overwriteOutput = True arcpy.CheckOutExtension("GeoStats") # Dynamické nastavení SnapRaster podle zadaného referenčního rastru snap_raster = os.path.join(input_workspace, DMR_raster) arcpy.env.snapRaster = snap_raster arcpy.AddMessage(f"Nastaven snapRaster: {snap_raster}") # Získání všech bodových vrstev tune_body_* arcpy.env.workspace = input_workspace point_layers = arcpy.ListFeatureClasses("tune_body_*") for point_layer in point_layers: point_layer_full = os.path.join(input_workspace, point_layer) arcpy.AddMessage(f"Zpracovávám vrstvu: {point_layer}") arcpy.AddMessage(f"Cesta k vrstvě: {point_layer_full}") try: # Kontrola, jestli vrstva má data num_points = int(arcpy.GetCount_management(point_layer_full).getOutput(0)) if num_points == 0: arcpy.AddWarning(f"Vrstva {point_layer} je prázdná. Přeskočeno.") continue # Nastavení parametrů interpolace - fixně radius = 60 nbr_min = 16 nbr_max = 16 sector_type = "FOUR_SECTOR" #FOUR_SECTORS search_neighborhood = ( f"NBRTYPE=StandardCircular RADIUS={radius} ANGLE=0 " f"NBR_MAX={nbr_max} NBR_MIN={nbr_min} SECTOR_TYPE={sector_type}" ) # Výstupní název rastru output_raster_name = f"EBK_{point_layer}" output_raster_path = os.path.join(input_workspace, output_raster_name) arcpy.AddMessage(f"Výstupní raster bude: {output_raster_path}") # Volání EBK interpolace arcpy.ga.EmpiricalBayesianKriging( in_features=point_layer_full, z_field="Elevation", out_ga_layer="", out_raster=output_raster_path, cell_size="0,25", search_neighborhood=search_neighborhood ) arcpy.AddMessage(f" ✅ Raster byl úspěšně uložen: {output_raster_path}") except arcpy.ExecuteError: error_message = arcpy.GetMessages(2) arcpy.AddError(f"Chyba při EBK interpolaci vrstvy {point_layer}: {error_message}") with open("ebk_errors.txt", "a") as f: f.write(f"Chyba EBK ({point_layer}): {error_message}\n") except Exception as e: arcpy.AddError(f"Neočekávaná chyba při EBK interpolaci vrstvy {point_layer}: {str(e)}") with open("ebk_errors.txt", "a") as f: f.write(f"Neočekávaná chyba EBK ({point_layer}): {str(e)}\n") arcpy.AddMessage("EBK interpolace s pevnými parametry byla dokončena.") # --- KONEC EBK --- # --- OŘEZÁVÁNÍ INTERPOLACÍ POMOCÍ buffer_pro_orez --- arcpy.AddMessage("Začínám ořezávat interpolace pomocí buffer_pro_orez...") # Workspace je stále input_workspace arcpy.env.workspace = input_workspace # Cesta k vrstvě buffer_pro_orez buffer_layer = os.path.join(input_workspace, "buffer_pro_orez") # Vytvoření vrstvy z buffer_pro_orez pro výběr arcpy.MakeFeatureLayer_management(buffer_layer, "buffer_layer_tmp") # Vyhledání všech interpolovaných rasterů ebk_rasters = arcpy.ListRasters("EBK_tune_body_*") for ebk_raster in ebk_rasters: # Přeskakujeme rastry, které už obsahují "_mask" if "_mask" in ebk_raster: continue ebk_raster_full = os.path.join(input_workspace, ebk_raster) # Získání čísla tůně try: index = int(ebk_raster.replace("EBK_tune_body_", "")) - 1 # protože tune_body_1 je první prvek seznamu except ValueError: arcpy.AddWarning(f"Nesprávný název rastru: {ebk_raster}, přeskočeno.") continue # Kontrola, zda index existuje v orig_fid_list if index < 0 or index >= len(orig_fid_list): arcpy.AddWarning(f"Index {index} mimo rozsah orig_fid_list pro raster {ebk_raster}, přeskočeno.") continue orig_fid = orig_fid_list[index] selection_query = f"ORIG_FID = {orig_fid}" # Výběr správného polygonu v buffer_layer_tmp podle ORIG_FID arcpy.SelectLayerByAttribute_management("buffer_layer_tmp", "NEW_SELECTION", selection_query) # Kontrola, jestli existuje nějaký výběr selection_count = int(arcpy.GetCount_management("buffer_layer_tmp").getOutput(0)) if selection_count == 0: arcpy.AddWarning(f"Nebyl nalezen žádný buffer s ORIG_FID = {orig_fid}. Raster {ebk_raster} bude přeskočen.") continue # Výstupní název output_clipped_raster_name = f"{ebk_raster}_mask" output_clipped_raster_path = os.path.join(input_workspace, output_clipped_raster_name) arcpy.AddMessage(f"Ořezávám raster {ebk_raster} podle bufferu s ORIG_FID = {orig_fid}.") # Spuštění ExtractByMask clipped_raster = arcpy.sa.ExtractByMask(ebk_raster_full, "buffer_layer_tmp") clipped_raster.save(output_clipped_raster_path) arcpy.AddMessage(f"Ořezaný raster uložen: {output_clipped_raster_path}") # Vyčištění vrstvy arcpy.Delete_management("buffer_layer_tmp") arcpy.AddMessage("Ořezávání všech interpolací bylo úspěšně dokončeno.") # --- KONEC OŘEZÁVÁNÍ --- # --- VÝSLEDNÝ RASTR (REFERENČNÍ + INTERPOLACE) --- arcpy.AddMessage("Začínám vytvářet novou kompletní vrstvu založenou na referenčním rastru a interpolacích...") arcpy.env.workspace = input_workspace mask_rasters = arcpy.ListRasters("*_mask") if not mask_rasters: arcpy.AddWarning("Nenalezeny žádné ořezané rastry.") else: reference_raster_path = os.path.join(input_workspace, DMR_raster) # Nejprve kopíruj referenční rastr final_raster_name = "final_combined" final_raster_path = os.path.join(input_workspace, final_raster_name) arcpy.CopyRaster_management( in_raster=reference_raster_path, out_rasterdataset=final_raster_path ) arcpy.AddMessage("Kopírování referenčního rastru dokončeno.") # Připrav seznam všech rastrů (nejprve interpolace, pak podklad) rasters_to_mosaic = mask_rasters + [final_raster_name] # Výstupní mozaika mosaic_output_name = "rastr_m" mosaic_output_path = os.path.join(input_workspace, mosaic_output_name) arcpy.MosaicToNewRaster_management( input_rasters=rasters_to_mosaic, output_location=input_workspace, raster_dataset_name_with_extension=mosaic_output_name, coordinate_system_for_the_raster="", pixel_type="32_BIT_FLOAT", cellsize = cell_size, # stejná velikost buňky number_of_bands=1, mosaic_method="FIRST", mosaic_colormap_mode="FIRST" ) arcpy.AddMessage(f"Nová kompletní mozaikovaná vrstva byla úspěšně vytvořena: {mosaic_output_path}") # --- KÓD PRO DOPLNĚNÍ NODATA HODNOT --- arcpy.AddMessage("Doplňuji NoData hodnot ve výsledném rastru...") # Zkontrolujeme, že je dostupná licence Spatial Analyst arcpy.CheckOutExtension("Spatial") try: # Vstupní rastr s NoData hodnotami input_raster = mosaic_output_path # "Vysledny_rastr_s_interpolacemi" # Výstupní rastr bez NoData hodnot output_raster_path = os.path.join(input_workspace, "VYSLEDNA_INTERPOLACE") # Výchozí rastr s NoData in_raster = arcpy.Raster(input_raster) # Iterativní přístup k vyplňování NoData hodnot # Vytvoříme dočasný rastr temp_raster = in_raster # Opakujeme proces několikrát pro postupné vyplňování větších mezer for i in range(3): # 3 iterace obvykle stačí pro menší mezery arcpy.AddMessage(f"Provádím {i+1}. iteraci vyplňování NoData hodnot...") # Použijeme nástroj FocalStatistics pro vyplnění NoData hodnot # s kruhovou sousedností o poloměru 3 buňky neighborhood = arcpy.sa.NbrCircle(3, "CELL") result = arcpy.sa.Con( arcpy.sa.IsNull(temp_raster), arcpy.sa.FocalStatistics(temp_raster, neighborhood, "MEAN", "DATA"), temp_raster ) # Aktualizujeme dočasný rastr pro další iteraci temp_raster = result # Uložení výsledku result.save(output_raster_path) arcpy.AddMessage(f"NoData hodnoty byly úspěšně doplněny: {output_raster_path}") except arcpy.ExecuteError: arcpy.AddError(f"Chyba při doplňování NoData hodnot: {arcpy.GetMessages(2)}") except Exception as e: arcpy.AddError(f"Neočekávaná chyba: {str(e)}") finally: # Vrátíme licenci arcpy.CheckInExtension("Spatial") # *** PŘIDÁNO: PŘEVOD POLYGONŮ NA CENTROIDY *** arcpy.AddMessage("Převádím polygony na centroidy...") centroid_layer = "vysledne_tune_centroidy" centroid_layer_full = os.path.join(input_workspace, centroid_layer) arcpy.FeatureToPoint_management( in_features=final_filtered_output, out_feature_class=centroid_layer_full, point_location="INSIDE" # Umístění centroidu uvnitř polygonu ) arcpy.AddMessage(f"Centroidy vytvořeny: {centroid_layer_full}") # Aktualizace výstupního parametru na filtrovanou vrstvu #arcpy.SetParameterAsText(4, output_raster_path) # Hlavní výstup - finální rastr #arcpy.SetParameterAsText(5, final_filtered_output) # tůně filtrované arcpy.AddMessage(f"Proces úspěšně dokončen!") feature_dataset_name = "vysledne_vrstvy" feature_dataset_path = f"{input_workspace}/{feature_dataset_name}" # Nejprve potřebujeme určit souřadnicový systém - použijeme ten, který má jedna z vrstev spatial_reference = arcpy.Describe(centroid_layer_full).spatialReference # Vytvoření feature datasetu arcpy.management.CreateFeatureDataset(input_workspace, feature_dataset_name, spatial_reference) # Kopírování vektorových vrstev do feature datasetu arcpy.management.CopyFeatures(centroid_layer_full, f"{feature_dataset_path}/Polygonova_vrstva_tuni_centroidy") arcpy.management.CopyFeatures(final_filtered_output, f"{feature_dataset_path}/Polygonova_vrstva_tuni") print(f"Feature dataset '{feature_dataset_name}' byl vytvořen a vrstvy byly přidány.") # Funkce pro mazání mezikrokových vrstev def delete_intermediate_layers(workspace, delete_flag): if delete_flag: arcpy.env.workspace = workspace # Explicitní nastavení pracovního prostoru arcpy.AddMessage(f"Pokouším se smazat mezikrokové vrstvy z: {workspace}") layers_to_delete = [ "nodata_raster", "seskupeny_raster", "tune_polygony", "hladke_polygony", "tune_buffer", "merged_polygons_1", "merged_polygons_2", "spatial_join", "prekryvajici_se_tune", "samostatne_tune", "vysledne_tune_konvexita", "buffer_pro_orez", "vysledne_tune_no_holes", "vysledne_tune_filtered", "final_combined", "vysledne_tune", "potencialni_tune_pred_sloucenim", "rastr_m" # Přidejte další názvy vrstev, které chcete smazat ] for layer in layers_to_delete: full_path = os.path.join(workspace, layer) if arcpy.Exists(full_path): try: arcpy.Delete_management(full_path) arcpy.AddMessage(f"Smazána mezikroková vrstva: {full_path}") except Exception as e: arcpy.AddWarning(f"Nepodařilo se smazat vrstvu {full_path}: {e}") else: arcpy.AddWarning(f"Vrstva nenalezena: {full_path}") else: arcpy.AddMessage("mezikrokové vrstvy NEBYLY smazány.") # Volání funkce pro mazání na konci skriptu delete_intermediate_layers(input_workspace, delete_intermediate_data) for layer in temp_layers: temp_path = os.path.join(arcpy.env.scratchFolder, layer) if arcpy.Exists(temp_path): try: arcpy.Delete_management(temp_path) arcpy.AddMessage(f"Dočasná vrstva {layer} byla odstraněna.") except: arcpy.AddWarning(f"Nepodařilo se odstranit dočasnou vrstvu {layer}.") except arcpy.ExecuteError: arcpy.AddError("Chyba při provádění nástroje.") arcpy.AddError(arcpy.GetMessages(2)) raise except Exception as e: arcpy.AddError(f"Neočekávaná chyba: {str(e)}") raise except arcpy.ExecuteError: arcpy.AddError("Během zpracování došlo k chybě v rámci ArcPy.") arcpy.AddError(arcpy.GetMessages(2)) except Exception as e: arcpy.AddError(f"Během zpracování došlo k neočekávané chybě: {str(e)}")