# Knickzone Extraction Tool (KET) Part2: Extraction of knickzones # Import arcpy module import arcpy from arcpy.sa import * from arcpy import env # Check out any necessary licenses arcpy.CheckOutExtension("spatial") arcpy.env.overwriteOutput = True env.qualifiedFieldNames = False #Input env.workspace = arcpy.GetParameterAsText(0) ##"C:/Dropbox/Files/py" in_dem = arcpy.GetParameterAsText(1) ##"C:/Dropbox/Files/py/wdem1.tif" aRas = arcpy.GetParameterAsText(2) ##"C:/Dropbox/Files/py/wcurv" #the final raster from previous tool kz_th = float(arcpy.GetParameterAsText(3)) ##-8 kz_ln_th = int(arcpy.GetParameterAsText(4)) ##20 #Minimum length of thresholds to keep in meter rvr_th = int(arcpy.GetParameterAsText(5)) ##5000 rvr = arcpy.GetParameterAsText(6) knickzones = arcpy.GetParameterAsText(7) #Output zonal_table = "zonal_" river_shp = rvr filldem = Fill(in_dem) if not rvr: fdr = FlowDirection(filldem, "NORMAL") fla = FlowAccumulation(fdr) rvr = "river_"+str(rvr_th) rvr = Con(fla > rvr_th,1) rvr.save("rvr_ras") river_shp = "river_"+str(rvr_th)+".shp" StreamToFeature(rvr, fdr, river_shp, "NO_SIMPLIFY") # Process: Extract by Mask rvr_aRas = ExtractByMask(aRas,rvr) # Process: KZ kz_Ras = Con(rvr_aRas < kz_th,1) # Process: INT kz_Ras_INT = Int(kz_Ras) # Process: Region Group kz_group = RegionGroup(kz_Ras_INT, "EIGHT", "WITHIN", "ADD_LINK", "") # Process: Raster to Polyline kz_polyline = arcpy.RasterToPolyline_conversion(kz_group, "kz_polyline", "NODATA", "0", "NO_SIMPLIFY", "VALUE") # Process: Add Field 'Length' arcpy.AddField_management(kz_polyline, "Length", "LONG", "", "", "", "", "NULLABLE", "NON_REQUIRED", "") # Process: Calculate length arcpy.CalculateField_management(kz_polyline, "Length", "!shape.length!", "PYTHON", "") # Process: Removing location smaller than "kz_ln_th" in length rows = arcpy.UpdateCursor(kz_polyline) for row in rows: if row.Length <= kz_ln_th: rows.deleteRow(row) # Process: Zonal statistics using the kz_polyline and filldem ZonalStatisticsAsTable(kz_polyline,"FID",filldem,zonal_table,"DATA","ALL") # Process: Join #Feature layer can be used in join but not a feature class or a shapefile arcpy.MakeFeatureLayer_management(kz_polyline, "tempLayer") arcpy.AddJoin_management("tempLayer","FID",zonal_table,"FID","KEEP_COMMON") kz_final = arcpy.CopyFeatures_management ("tempLayer",knickzones) ##arcpy.FeatureClassToFeatureClass_conversion ("tempLayer","/","kz_final.shp") arcpy.DeleteField_management(kz_final,drop_field="AREA;SUM;ARCID;GRID_CODE;FROM_NODE;TO_NODE;Rowid_;FID_1") arcpy.ExportXYv_stats(kz_final,"FID;Length;COUNT;MIN;MAX;RANGE;MEAN;STD", "COMMA","Final_table.csv","ADD_FIELD_NAMES") print("Done")