# -*- coding: utf-8 -*- """ Zadání: Vytvoř v arcpy pro ArcGIS Pro 3 shlukovací kód, který nashlukuje pověřené obce podle obslužnosti, kde pokud se obslužná vzdálenost obcí překrývá, patří do stejného shluku. Vstupem je shapefile polygonů pověřených obcí a shapefile silnic. Obslužná vzdálenost se určí síťovou analýzou silnic a centroidů pověřených obcí – procházíme každý centroid pověřené obce, zjistíme do jakých pověřených obcí dosahuje obslužnost, sjednotíme překrývající se obslužné vzdálenosti a do shapefile pověřených obcí přidáme sloupec s názvem podle obslužné vzdálenosti a s hodnotou příslušného shluku pro danou pověřenou obec. Hodnota obslužné vzdálenosti se zadá inputem – default hodnota je 5000 m. """ import arcpy import os # Nastavení vstupních parametrů # ------------------------------------ # Shapefile polygonů pověřených obcí obce_shp = r"C:\cesta\k\poverene_obce.shp" # Shapefile silnic silnice_shp = r"C:\cesta\k\silnice.shp" # Geodatabáze pro uložení Network Datasetu gdb_path = r"C:\cesta\k\silnice.gdb" # Název Network Datasetu network_name = "Silnice_ND" # Výstupní shapefile s přidaným sloupcem vystup_shp = r"C:\cesta\k\vystup_obce.shp" # Obslužná vzdálenost v metrech (default je 5000 m) obsluzna_vzdalenost = 5000 # Kontrola licence Network Analyst # ------------------------------------ # Zjistíme, zda je dostupná licence Network Analyst if arcpy.CheckExtension("Network") == "Available": arcpy.CheckOutExtension("Network") else: arcpy.AddError("Rozšíření Network Analyst není dostupné.") exit() # Krok 1: Převod silnic do File Geodatabase # ------------------------------------ # Pokud geodatabáze neexistuje, vytvoříme ji if not arcpy.Exists(gdb_path): arcpy.management.CreateFileGDB(os.path.dirname(gdb_path), os.path.basename(gdb_path)) # Cesta k feature class silnic v geodatabázi silnice_fc = os.path.join(gdb_path, "Silnice") # Pokud feature class silnic neexistuje, importujeme ji if not arcpy.Exists(silnice_fc): arcpy.conversion.FeatureClassToFeatureClass(silnice_shp, gdb_path, "Silnice") # Krok 2: Vytvoření Network Datasetu # ------------------------------------ # Nastavíme pracovní prostor na geodatabázi arcpy.env.workspace = gdb_path # Cesta k Network Datasetu silnice_nd = os.path.join(gdb_path, network_name) # Pokud Network Dataset neexistuje, vytvoříme jej if not arcpy.Exists(silnice_nd): # Vytvoříme novou síťovou datasetovou konfiguraci arcpy.na.CreateNetworkDataset(os.path.basename(gdb_path), network_name) # Přidáme Edge Source (silnice) do Network Datasetu arcpy.na.AddEdgeSource(silnice_nd, "Silnice") # Přidáme atribut délky (Length) jako náklad (Cost) arcpy.na.AddNetworkAttribute(silnice_nd, "Length", "Cost", units="Meters", use_by_default="True") # Sestavíme Network Dataset arcpy.na.BuildNetwork(silnice_nd) # Krok 3: Vytvoření centroidů pověřených obcí # ------------------------------------ # Vytvoříme centroidy z polygonů obcí centroidy_obci = r"in_memory\centroidy_obci" arcpy.FeatureToPoint_management(obce_shp, centroidy_obci, "INSIDE") # Přidání unikátního ID do obcí a centroidů arcpy.management.AddField(obce_shp, "ObecID", "LONG") arcpy.management.CalculateField(obce_shp, "ObecID", "!OBJECTID!", "PYTHON3") arcpy.management.AddField(centroidy_obci, "ObecID", "LONG") arcpy.management.CalculateField(centroidy_obci, "ObecID", "!OBJECTID!", "PYTHON3") # Krok 4: Vytvoření vrstvy pro analýzu obslužné oblasti # ------------------------------------ # Vytvoříme Service Area Layer pro síťovou analýzu outNALayer = arcpy.na.MakeServiceAreaLayer( in_network_dataset=silnice_nd, # Network Dataset silnic out_network_analysis_layer="ObsluznaOblast", # Název výstupní vrstvy impedance_attribute="Length", # Atribut nákladů (délka) travel_from_to="TRAVEL_FROM", # Směr cesty (z centroidů) default_break_values=[obsluzna_vzdalenost], # Obslužná vzdálenost polygon_type="DETAILED_POLYS", # Typ polygonů (detailní) merge="NO_MERGE", # Neslučovat výsledky nesting_type="RINGS", # Typ vnoření (kruhy) line_type="NO_LINES", # Nevytvářet linie overlap="OVERLAP", # Překrývání povoleno split="SPLIT" # Dělení polygonů ) # Získání objektu vrstvy a názvů podvrstev na_layer = outNALayer.getOutput(0) sub_layers = arcpy.na.GetNAClassNames(na_layer) facilities_layer_name = sub_layers["Facilities"] # Krok 5: Přidání centroidů jako zařízení (Facilities) # ------------------------------------ # Přidáme centroidy obcí do vrstvy síťové analýzy arcpy.na.AddLocations( in_network_analysis_layer=na_layer, sub_layer=facilities_layer_name, in_table=centroidy_obci, field_mappings="", search_tolerance="", append="CLEAR" # Vymažeme předchozí data ) # Krok 6: Vyřešení analýzy obslužné oblasti # ------------------------------------ # Spustíme síťovou analýzu pro výpočet obslužných oblastí arcpy.na.Solve(na_layer) # Krok 7: Získání vrstev polygonů obslužných oblastí # ------------------------------------ # Získáme vrstvu polygonů z výsledků analýzy polygons_layer_name = sub_layers["SAPolygons"] polygons_sublayer = na_layer.listLayers(polygons_layer_name)[0] # Kopírujeme polygony obslužných oblastí do in-memory vrstvy obsluzne_polygony = r"in_memory\obsluzne_polygony" arcpy.management.CopyFeatures(polygons_sublayer, obsluzne_polygony) # Krok 8: Sjednocení překrývajících se obslužných oblastí # ------------------------------------ # Sloučíme překrývající se polygony do jednotných shluků sloucene_oblasti = r"in_memory\sloucene_oblasti" arcpy.management.Dissolve( in_features=obsluzne_polygony, out_feature_class=sloucene_oblasti, dissolve_field="", multi_part="MULTI_PART" ) # Přidání pole pro shluk (ShlukID) arcpy.management.AddField(sloucene_oblasti, "ShlukID", "LONG") arcpy.management.CalculateField( in_table=sloucene_oblasti, field="ShlukID", expression="!OBJECTID!", expression_type="PYTHON3" ) # Krok 9: Prostorové spojení centroidů s shluky # ------------------------------------ # Přiřadíme centroidům obcí příslušný ShlukID na základě prostorového překrytí centroidy_s_shlukem = r"in_memory\centroidy_s_shlukem" arcpy.analysis.SpatialJoin( target_features=centroidy_obci, join_features=sloucene_oblasti, out_feature_class=centroidy_s_shlukem, join_operation="JOIN_ONE_TO_ONE", match_option="INTERSECT" ) # Krok 10: Vytvoření výstupního shapefile obcí s přidaným shlukem # ------------------------------------ # Kopírujeme původní shapefile obcí do nového výstupního souboru arcpy.management.CopyFeatures(obce_shp, vystup_shp) # Přidáme pole ShlukID do výstupního shapefile pomocí JoinField arcpy.management.JoinField( in_data=vystup_shp, in_field="ObecID", join_table=centroidy_s_shlukem, join_field="ObecID", fields=["ShlukID"] ) # Přejmenujeme pole ShlukID podle obslužné vzdálenosti nazev_pole = f"Shluk_{obsluzna_vzdalenost}m" arcpy.management.AlterField( in_table=vystup_shp, field="ShlukID", new_field_name=nazev_pole, new_field_alias=nazev_pole ) # Krok 11: Vyčištění dočasných dat # ------------------------------------ # Odstraníme dočasné vrstvy z paměti arcpy.management.Delete(centroidy_obci) arcpy.management.Delete(obsluzne_polygony) arcpy.management.Delete(sloucene_oblasti) arcpy.management.Delete(centroidy_s_shlukem) # Vracení licence Network Analyst arcpy.CheckInExtension("Network") # Informace o dokončení procesu print(f"Shlukování dokončeno. Výstup uložen v: {vystup_shp}")