------------------------------------------------------------------- import arcpy from arcpy.sa import * import os import sys import datetime import time import FindBraidedNetwork import BRAT_Braid_Handler from SupportingFunctions import make_layer, make_folder, getUUID, find_relative_path, write_xml_element_with_path import XMLBuilder import SupportingFunctions import importlib importlib.reload(XMLBuilder) XMLBuilder = XMLBuilder.XMLBuilder importlib.reload(FindBraidedNetwork) importlib.reload(BRAT_Braid_Handler) def main( proj_path, seg_network, in_DEM, flow_acc, coded_veg, coded_hist, valley_bottom, road, railroad, canal, landuse, ownership, perennial_network, out_name, description, find_clusters, should_segment_network, segment_by_ownership, is_verbose): """ Calculates, for each stream network segment, the attributes needed to trun the BRAT tools. :param proj_path: Path to the BRAT project folder. :param seg_network: The segmented (300m) network. :param in_DEM: The DEM for the entire project. :param flow_acc: The flow accumulation raster. :param coded_veg: The landfire EVT layer :param coded_hist: The landfire BPS layer :param valley_bottom: The valley bottom polygon that is associated with the input stream network. :param road: The shapefile containing all roads for the area :param railroad: The shapefile containing all railroads for the area :param canal: The shapefile containing all canals for the area :param landuse: The raster containing all land use data for the area :param ownership: The shapefile contating all land ownership data for the area :param perennial_network: The perennial network of streams :param out_name: The name for the output BRAT Table. :param description: A short description of the run that will be added to the XML :param find_clusters: If true, this option will create a ClusterID field and populate it :param should_segment_network: If true, this option divides reaches based on the roads input. :param segment_by_ownership: If true, this option divides reaches based on the land ownership input. :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :return: """ if flow_acc == "None": flow_acc = None if perennial_network == "None": perennial_network = None find_clusters = parse_input_bool(find_clusters) should_segment_network = parse_input_bool(should_segment_network) segment_by_ownership = parse_input_bool(segment_by_ownership) is_verbose = parse_input_bool(is_verbose) scratch = 'in_memory' #arcpy.env.workspace = scratch arcpy.env.overwriteOutput = True arcpy.env.outputZFlag = "Disabled" arcpy.env.outputMFlag = "Disabled" arcpy.CheckOutExtension("Spatial") # --check input projections-- validate_inputs(seg_network, road, railroad, canal, is_verbose) # name and create output folder new_output_folder, intermediate_folder, seg_network_copy = build_output_folder(proj_path, out_name, seg_network, road, should_segment_network, ownership, segment_by_ownership, is_verbose) # --check input network fields-- # add flowline reach id field ('ReachID') if it doens't already exist # this field allows for more for more 'stable' joining fields = [f.name for f in arcpy.ListFields(seg_network_copy)] if 'ReachID' not in fields: arcpy.AddField_management(seg_network_copy, 'ReachID', 'LONG') with arcpy.da.UpdateCursor(seg_network_copy, ['FID', 'ReachID']) as cursor: for row in cursor: row[1] = row[0] cursor.updateRow(row) # --create network buffers for analyses-- # create 'Buffers' folder if it doesn't exist buffers_folder = make_folder(intermediate_folder, "01_Buffers") # create network segment midpoints if is_verbose: arcpy.AddMessage("Finding network segment midpoints...") midpoints = arcpy.FeatureVerticesToPoints_management(seg_network_copy, scratch + "/midpoints", "MID") # remove unwanted fields from midpoints fields = arcpy.ListFields(midpoints) keep = ['ReachID'] drop = [] for field in fields: if not field.required and field.name not in keep and field.type != 'Geometry': drop.append(field.name) if len(drop) > 0: arcpy.DeleteField_management(midpoints, drop) if is_verbose: arcpy.AddMessage("Making buffers...") # create midpoint 100 m buffer midpoint_buffer = arcpy.Buffer_analysis(midpoints, scratch + "/midpoint_buffer", "100 Meters") # create network 30 m buffer buf_30m = os.path.join(buffers_folder, "buffer_30m.shp") arcpy.Buffer_analysis(seg_network_copy, buf_30m, "30 Meters", "", "ROUND") # create network 100 m buffer buf_100m = os.path.join(buffers_folder, "buffer_100m.shp") arcpy.Buffer_analysis(seg_network_copy, buf_100m, "100 Meters", "", "ROUND") # run geo attributes function arcpy.AddMessage('Adding "iGeo" attributes to network...') igeo_attributes(seg_network_copy, in_DEM, flow_acc, midpoint_buffer, scratch, is_verbose) # run vegetation attributes function arcpy.AddMessage('Adding "iVeg" attributes to network...') iveg_attributes(coded_veg, coded_hist, buf_100m, buf_30m, seg_network_copy, scratch, is_verbose) # find points of diversion if canals are defined if canal is not None: diversion_pts = find_points_of_diversion(canal, seg_network_copy, perennial_network, proj_path, is_verbose) else: diversion_pts = None # run ipc attributes function if conflict layers are defined by user if road is not None and valley_bottom is not None: arcpy.AddMessage('Adding "iPC" attributes to network...') ipc_attributes(seg_network_copy, road, railroad, canal, valley_bottom, ownership, diversion_pts, buf_30m, buf_100m, landuse, scratch, proj_path, is_verbose) if perennial_network is not None: find_is_perennial(seg_network_copy, perennial_network) handle_braids(seg_network_copy, canal, proj_path, find_clusters, perennial_network, is_verbose) # run write xml function arcpy.AddMessage('Writing project xml...') DrAr = find_dr_ar(flow_acc, in_DEM) trib_code_folder = os.path.dirname(os.path.abspath(__file__)) symbology_folder = os.path.join(trib_code_folder, 'BRATSymbology') flow_accumulation_sym_layer = os.path.join(symbology_folder, "Flow_Accumulation.lyr") make_layer(os.path.dirname(DrAr), DrAr, "Flow Accumulation", symbology_layer=flow_accumulation_sym_layer, is_raster=True) make_layers(seg_network_copy, diversion_pts) write_xml(new_output_folder, coded_veg, coded_hist, seg_network, in_DEM, valley_bottom, landuse, DrAr, road, railroad, canal, buf_30m, buf_100m, seg_network_copy, description) run_tests(seg_network_copy, is_verbose) arcpy.CheckInExtension("spatial") def find_is_perennial(seg_network_copy, perennial_network): """ Adds the IsPerennial attribute :param seg_network_copy: The BRAT Table output :param perennial_network: The input stream network that only contains perennial networks :return: """ arcpy.AddField_management(seg_network_copy, "IsPeren", "SHORT") arcpy.CalculateField_management(seg_network_copy, "IsPeren", 0, "PYTHON") seg_network_layer = "seg_network_lyr" perennial_network_layer = "perennial_network_layer" arcpy.MakeFeatureLayer_management(seg_network_copy, seg_network_layer) arcpy.MakeFeatureLayer_management(perennial_network, perennial_network_layer) arcpy.SelectLayerByLocation_management(seg_network_layer, "SHARE_A_LINE_SEGMENT_WITH", perennial_network_layer, '', "NEW_SELECTION") arcpy.CalculateField_management(seg_network_layer,"IsPeren",1,"PYTHON") def find_dr_ar(flow_acc, in_DEM): """ Finds the path to the drainage area. :param flow_acc: The flow accumulation raster. This may be empty :param in_DEM: The DEM for the entire area :return: The file path to the Drainage Area Raster """ if flow_acc is None: DrArea = os.path.join(os.path.join(os.path.dirname(in_DEM), "Flow", "DrainArea_sqkm.tif")) else: DrArea = os.path.join(os.path.join(os.path.dirname(in_DEM), "Flow"), os.path.basename(flow_acc)) return DrArea def build_output_folder(proj_path, out_name, seg_network, road, should_segment_network, ownership, segment_by_ownership, is_verbose): """ Builds the outputs folder where everything will be put :param proj_path: Path to the BRAT project folder. :param out_name: The name for the output BRAT Table. :param seg_network: The segmented (300m) network. :param road: The shapefile containing all roads for the area :param should_segment_network: If true, this option divides reaches based on the roads input. :param ownership: The shapefile contating all land ownership data for the area :param segment_by_ownership: The shapefile containing all land ownership data for the area :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :return: The created output folder, the created intermediate folder, and the potentially segmented network """ if is_verbose: arcpy.AddMessage("Building folder structure...") master_outputs_folder = os.path.join(proj_path, "Outputs") if not os.path.exists(master_outputs_folder): os.mkdir(master_outputs_folder) j = 1 str_num = '01' new_output_folder = os.path.join(master_outputs_folder, "Output_" + str_num) while os.path.exists(new_output_folder): j += 1 if j > 9: str_num = str(j) else: str_num = "0" + str(j) new_output_folder = os.path.join(master_outputs_folder, "Output_" + str_num) os.mkdir(new_output_folder) intermediate_folder = make_folder(new_output_folder, "01_Intermediates") # copy input segment network to output folder if out_name.endswith('.shp'): seg_network_copy = os.path.join(intermediate_folder, out_name) else: seg_network_copy = os.path.join(intermediate_folder, out_name + ".shp") # segment network by roads if should_segment_network is True if should_segment_network: segment_by_roads(seg_network, seg_network_copy, road, is_verbose) else: arcpy.CopyFeatures_management(seg_network, seg_network_copy) # segment network by roads if segment_by_ownership is True if segment_by_ownership: segment_network_by_ownership(seg_network_copy, ownership, is_verbose) else: pass return new_output_folder, intermediate_folder, seg_network_copy def segment_by_roads(seg_network, seg_network_copy, roads, is_verbose): """ Segments the seg_network by roads, and puts segmented network at seg_network_copy :param seg_network: Path to the seg_network that we want to segment further :param seg_network_copy: Path to where we want the new network to go :param roads: The shape file we use to segment :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :return: """ arcpy.AddMessage("Segmenting network by roads...") temp_network = os.path.join(os.path.dirname(seg_network_copy), "temp.shp") temp_layer = "temp_lyr" temp_seg_network_layer = "seg_network_lyr" arcpy.FeatureToLine_management([seg_network, roads], temp_network) arcpy.MakeFeatureLayer_management(temp_network, temp_layer) arcpy.MakeFeatureLayer_management(seg_network, temp_seg_network_layer) arcpy.SelectLayerByLocation_management(temp_layer, "WITHIN", temp_seg_network_layer) arcpy.CopyFeatures_management(temp_layer, seg_network_copy) delete_with_arcpy([temp_layer, temp_seg_network_layer, temp_network]) add_reach_dist(seg_network, seg_network_copy, is_verbose) def segment_network_by_ownership(seg_network_copy, ownership, is_verbose): """ Segments the seg_network by ownership, and puts segmented network at seg_network_copy :param seg_network_copy: Path to the seg_network that we want to segment further :param ownership: The shape file we use to segment :param is_verbose: Specifies whether to provide messages :return: """ arcpy.AddMessage("Segmenting network by ownership...") temp_seg_network_copy = os.path.join(os.path.dirname(seg_network_copy), "temp_seg_network_copy.shp") temp_network = os.path.join(os.path.dirname(seg_network_copy), "temp.shp") temp_layer = "temp_lyr" temp_seg_network_copy_layer = "seg_network_lyr" arcpy.CopyFeatures_management(seg_network_copy, temp_seg_network_copy) arcpy.FeatureToLine_management([seg_network_copy, ownership], temp_network) arcpy.MakeFeatureLayer_management(temp_network, temp_layer) arcpy.MakeFeatureLayer_management(seg_network_copy, temp_seg_network_copy_layer) arcpy.SelectLayerByLocation_management(temp_layer, "WITHIN", temp_seg_network_copy_layer) arcpy.CopyFeatures_management(temp_layer, seg_network_copy) add_reach_dist(temp_seg_network_copy, seg_network_copy, is_verbose) delete_with_arcpy([temp_layer, temp_seg_network_copy_layer, temp_network, temp_seg_network_copy]) def add_reach_dist(seg_network, seg_network_copy, is_verbose): """ Adds reach distance field to the network :param seg_network: The original segmented network :param seg_network_copy: The copy of the segmented network created by build_output_folder :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :return: """ if is_verbose: arcpy.AddMessage("Calculating ReachDist...") fields = [f.name for f in arcpy.ListFields(seg_network_copy)] if 'ReachID' in fields: with arcpy.da.UpdateCursor(seg_network_copy, ['FID', 'ReachID']) as cursor: for row in cursor: row[1] = row[0] cursor.updateRow(row) # get distance along route (LineID) for segment midpoints midpoints = arcpy.FeatureVerticesToPoints_management(seg_network_copy, 'in_memory/midpoints', "MID") seg_network_dissolve = arcpy.Dissolve_management(seg_network, 'in_memory/seg_network_dissolve', 'StreamID', '', 'SINGLE_PART', 'UNSPLIT_LINES') arcpy.AddField_management(seg_network_dissolve, 'From_', 'DOUBLE') arcpy.AddField_management(seg_network_dissolve, 'To_', 'DOUBLE') with arcpy.da.UpdateCursor(seg_network_dissolve, ['SHAPE@Length', 'From_', 'To_']) as cursor: for row in cursor: row[1] = 0.0 row[2] = row[0] cursor.updateRow(row) arcpy.CreateRoutes_lr(seg_network_dissolve, 'StreamID', 'in_memory/flowline_route', 'TWO_FIELDS', 'From_', 'To_') route_tbl = arcpy.LocateFeaturesAlongRoutes_lr(midpoints, 'in_memory/flowline_route', 'StreamID', 1.0, os.path.join(os.path.dirname(seg_network_copy), 'tbl_Routes.dbf'), 'RID POINT MEAS') dist_dict = {} # add reach id distance values to dictionary with arcpy.da.SearchCursor(route_tbl, ['ReachID', 'MEAS']) as cursor: for row in cursor: dist_dict[row[0]] = row[1] # populate dictionary value to output field by ReachID arcpy.AddField_management(seg_network_copy, 'ReachDist', 'DOUBLE') with arcpy.da.UpdateCursor(seg_network_copy, ['ReachID', 'ReachDist']) as cursor: for row in cursor: aKey = row[0] row[1] = dist_dict[aKey] cursor.updateRow(row) arcpy.Delete_management('in_memory') def zonalStatsWithinBuffer(buffer, ras, stat_type, stat_field, out_fc, out_FC_field, scratch): """ Calculate zonal statistics within buffer function :param buffer: The buffer around the stream :param ras: The DEM raster :param stat_type: The type of statistic to be calculated (MAXIMUM, MINIMUM, etc.) :param stat_field: The name of the field containing statistics to be calculated :param out_fc: The feature class to output to. :param out_FC_field: The field within the output feature class to output to. :param scratch: The current workspace :return: """ # get input raster stat value within each buffer # note: zonal stats as table does not support overlapping polygons so we will check which # reach buffers output was produced for and which we need to run tool on again stat_tbl = arcpy.sa.ZonalStatisticsAsTable(buffer, 'ReachID', ras, os.path.join(scratch, 'statTbl'), 'DATA', stat_type) # get list of segment buffers where zonal stats tool produced output have_stat_list = [row[0] for row in arcpy.da.SearchCursor(stat_tbl, 'ReachID')] # create dictionary to hold all reach buffer min dem z values stat_dict = {} # add buffer raster stat values to dictionary with arcpy.da.SearchCursor(stat_tbl, ['ReachID', stat_field]) as cursor: for row in cursor: stat_dict[row[0]] = row[1] # create list of overlapping buffer reaches (i.e., where zonal stats tool did not produce output) need_stat_list = [] with arcpy.da.SearchCursor(buffer, ['ReachID']) as cursor: for row in cursor: if row[0] not in have_stat_list: need_stat_list.append(row[0]) # run zonal stats until we have output for each overlapping buffer segment stat = None tmp_buff_lyr = None num_broken_repetitions = 0 BROKEN_REPS_ALLOWED = 5 while len(need_stat_list) > 0: # create tuple of segment ids where still need raster values need_stat = () for reach in need_stat_list: if reach not in need_stat: need_stat += (reach,) # use the segment id tuple to create selection query and run zonal stats tool if len(need_stat) == 1: quer = '"ReachID" = ' + str(need_stat[0]) else: quer = '"ReachID" IN ' + str(need_stat) tmp_buff_lyr = arcpy.MakeFeatureLayer_management(buffer, 'tmp_buff_lyr') arcpy.SelectLayerByAttribute_management(tmp_buff_lyr, 'NEW_SELECTION', quer) stat = arcpy.sa.ZonalStatisticsAsTable(tmp_buff_lyr, 'ReachID', ras, os.path.join(scratch, 'stat'), 'DATA', stat_type) # add segment stat values from zonal stats table to main dictionary with arcpy.da.SearchCursor(stat, ['ReachID', stat_field]) as cursor: for row in cursor: stat_dict[row[0]] = row[1] # create list of reaches that were run and remove from 'need to run' list have_stat_list2 = [row[0] for row in arcpy.da.SearchCursor(stat, 'ReachID')] if len(have_stat_list2) == 0: num_broken_repetitions += 1 if num_broken_repetitions >= BROKEN_REPS_ALLOWED: warning_message = "While calculating " + out_FC_field + ", the tool ran into an error. The following " warning_message += "ReachIDs did not recieve correct values:\n" for reach_id in need_stat_list: if reach_id == need_stat_list[-1]: warning_message += "and " warning_message += str(reach_id) if reach_id != need_stat_list[-1]: warning_message += ", " warning_message += "\n" arcpy.AddWarning(warning_message) for reach_id in need_stat_list: stat_dict[reach_id] = 0 need_stat_list = [] for reach in have_stat_list2: need_stat_list.remove(reach) # populate dictionary value to output field by ReachID with arcpy.da.UpdateCursor(out_fc, ['ReachID', out_FC_field]) as cursor: for row in cursor: try: aKey = row[0] row[1] = stat_dict[aKey] cursor.updateRow(row) except: pass stat_dict.clear() # delete temp fcs, tbls, etc. #items = [statTbl, haveStatList, haveStatList2, needStatList, stat, tmp_buff_lyr, needStat] items = [stat_tbl, stat, tmp_buff_lyr] for item in items: if item is not None: arcpy.Delete_management(item) def igeo_attributes(out_network, in_DEM, flow_acc, midpoint_buffer, scratch, is_verbose): """ calculates min and max elevation, length, slope, and drainage area for each flowline segment :param out_network: The output netwrok to add fields to. :param in_DEM: The DEM raster. :param flow_acc: Th eflow accumulation raster :param midpoint_buffer: The buffer created from midpoints :param scratch: The current workspace :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :return: Drainage Area """ # if fields already exist, delete them fields = [f.name for f in arcpy.ListFields(out_network)] drop = ["iGeo_ElMax", "iGeo_ElMin", "iGeo_Len", "iGeo_Slope", "iGeo_DA"] for field in fields: if field in drop: arcpy.DeleteField_management(out_network, field) # add flowline segment id field ('ReachID') for more 'stable' joining if 'ReachID' not in fields: arcpy.AddField_management(out_network, 'ReachID', 'LONG') with arcpy.da.UpdateCursor(out_network, ['FID', 'ReachID']) as cursor: for row in cursor: row[1] = row[0] cursor.updateRow(row) if is_verbose: arcpy.AddMessage("Preprocessing DEM...") # --smooth input dem by 3x3 cell window-- # define raster environment settings desc = arcpy.Describe(in_DEM) arcpy.env.extent = desc.Extent arcpy.env.outputCoordinateSystem = desc.SpatialReference arcpy.env.cellSize = desc.meanCellWidth # calculate mean z over 3x3 cell window neighborhood = NbrRectangle(3, 3, "CELL") tmp_dem = FocalStatistics(in_DEM, neighborhood, 'MEAN') # clip smoothed dem to input dem DEM = ExtractByMask(tmp_dem, in_DEM) # function to attribute start/end elevation (dem z) to each flowline segment def zSeg(vertex_type, out_field): if is_verbose: arcpy.AddMessage("Calculating values for " + out_field + "...") # create start/end points for each flowline reach segment tmp_pts = os.path.join(scratch, 'tmp_pts') arcpy.FeatureVerticesToPoints_management(out_network, tmp_pts, vertex_type) # create 30 meter buffer around each start/end point tmp_buff = os.path.join(scratch, 'tmp_buff') arcpy.Buffer_analysis(tmp_pts, tmp_buff, '30 Meters') # get min dem z value within each buffer arcpy.AddField_management(out_network, out_field, "DOUBLE") zonalStatsWithinBuffer(tmp_buff, DEM, 'MINIMUM', 'MIN', out_network, out_field, scratch) # delete temp fcs, tbls, etc. items = [tmp_pts, tmp_buff] for item in items: arcpy.Delete_management(item) # run zSeg function for start/end of each network segment zSeg('START', 'iGeo_ElMax') zSeg('END', 'iGeo_ElMin') # calculate network reach slope arcpy.AddField_management(out_network, "iGeo_Len", "DOUBLE") arcpy.CalculateField_management(out_network, "iGeo_Len", '!shape.length@meters!', "PYTHON_9.3") arcpy.AddField_management(out_network, "iGeo_Slope", "DOUBLE") with arcpy.da.UpdateCursor(out_network, ["iGeo_ElMax", "iGeo_ElMin", "iGeo_Len", "iGeo_Slope"]) as cursor: if is_verbose: arcpy.AddMessage("Calculating iGeo_Slope...") for row in cursor: try: row[3] = (abs(row[0] - row[1]))/row[2] except ZeroDivisionError: row[3] = 0 if row[3] == 0.0: row[3] = 0.0001 cursor.updateRow(row) # get DA values if flow_acc is None: arcpy.AddMessage("Calculating drainage area...") calc_drain_area(DEM, in_DEM) elif not os.path.exists(os.path.dirname(in_DEM) + "/Flow"): # if there's no folder for the flow accumulation, make one os.mkdir(os.path.dirname(in_DEM) + "/Flow") if is_verbose: arcpy.AddMessage("Copying drainage area raster...") arcpy.CopyRaster_management(flow_acc, os.path.dirname(in_DEM) + "/Flow/" + os.path.basename(flow_acc)) DrArea = find_dr_ar(flow_acc, in_DEM) # Todo: check this bc it seems wrong to pull from midpoint buffer # add drainage area 'iGeo_DA' field to flowline network arcpy.AddField_management(out_network, "iGeo_DA", "DOUBLE") # get max drainage area within 100 m midpoint buffer if is_verbose: arcpy.AddMessage("Calculating iGeo_DA...") zonalStatsWithinBuffer(midpoint_buffer, DrArea, "MAXIMUM", 'MAX', out_network, "iGeo_DA", scratch) # replace '0' drainage area values with tiny value with arcpy.da.UpdateCursor(out_network, ["iGeo_DA"]) as cursor: for row in cursor: if row[0] == 0: row[0] = 0.00000001 cursor.updateRow(row) return DrArea def iveg_attributes(coded_veg, coded_hist, buf_100m, buf_30m, out_network, scratch, is_verbose): """ Calculates both existing and potential mean vegetation value within 30 m and 100 m buffer of each stream segment :param coded_veg: The coded existing vegetation raster :param coded_hist: The coded historic vegetation raster :param buf_100m: The 100m stream buffer :param buf_30m: The 30m stream buffer :param out_network: The output network that data will be added to. :param scratch: The current workspace :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :return: """ # if fields already exist, delete them fields = [f.name for f in arcpy.ListFields(out_network)] drop = ["iVeg_100EX", "iVeg_30EX", "iVeg100Hpe", "iVeg_30Hpe"] for field in fields: if field in drop: arcpy.DeleteField_management(out_network, field) # --existing vegetation values-- if is_verbose: arcpy.AddMessage("Creating current veg lookup raster...") veg_lookup = Lookup(coded_veg, "VEG_CODE") # add mean veg value 'iVeg100EX' field to flowline network arcpy.AddField_management(out_network, "iVeg100EX", "DOUBLE") # get mean existing veg value within 100 m buffer if is_verbose: arcpy.AddMessage("Calculating iVeg100EX...") zonalStatsWithinBuffer(buf_100m, veg_lookup, 'MEAN', 'MEAN', out_network, "iVeg100EX", scratch) # add mean veg value 'iVeg_VT30EX' field to flowline network arcpy.AddField_management(out_network, "iVeg_30EX", "DOUBLE") # get mean existing veg value within 30 m buffer if is_verbose: arcpy.AddMessage("Calculating iVeg_30EX...") zonalStatsWithinBuffer(buf_30m, veg_lookup, 'MEAN', 'MEAN', out_network, "iVeg_30EX", scratch) # delete temp fcs, tbls, etc. items = [veg_lookup] for item in items: arcpy.Delete_management(item) # --historic (i.e., potential) vegetation values-- if is_verbose: arcpy.AddMessage("Creating historic veg lookup raster...") hist_veg_lookup = Lookup(coded_hist, "VEG_CODE") # add mean veg value 'iVeg100Hpe' field to flowline network arcpy.AddField_management(out_network, "iVeg100Hpe", "DOUBLE") # get mean potential veg value within 100 m buffer if is_verbose: arcpy.AddMessage("Calculating iVeg100Hpe...") zonalStatsWithinBuffer(buf_100m, hist_veg_lookup, 'MEAN', 'MEAN', out_network, "iVeg100Hpe", scratch) # add mean veg value 'iVeg_30Hpe' field to flowline network arcpy.AddField_management(out_network, "iVeg_30Hpe", "DOUBLE") # get mean potential veg value within 30 m buffer if is_verbose: arcpy.AddMessage("Calculating iVeg_30Hpe...") zonalStatsWithinBuffer(buf_30m, hist_veg_lookup, 'MEAN', 'MEAN', out_network, "iVeg_30Hpe", scratch) # delete temp fcs, tbls, etc. items = [hist_veg_lookup] for item in items: arcpy.Delete_management(item) def find_points_of_diversion(canal, network, perennial_network, proj_path, is_verbose): """ Finds potential points where water is being diverted into canals from the stream network by finding intersections between the canal and network :param: canal: input canals network :param: network: input stream network :param: perennial_network: input perennial stream network, will be used preferentially to the full network :param: is verbose: True/False for whether to return messages in ArcMap :return: points of diversion, the points where the canals and network intersect """ # find temp directory, or make if not present temp_dir = os.path.join(proj_path, 'Temp') if not os.path.exists(temp_dir): make_temp_dir(proj_path, is_verbose) # find points of diversion (intersection between perennial stream and canals if is_verbose: arcpy.AddMessage("Finding points of diversion...") # name output and place in canals folder canal_folder = os.path.dirname(canal) diversion_points = canal_folder + "\\points_of_diversion.shp" if not os.path.exists(diversion_points): # dissolve canals into single feature canal_dissolve = temp_dir + "\\canals_dissolved.shp" arcpy.Dissolve_management(canal, canal_dissolve, '', '', 'SINGLE_PART', 'UNSPLIT_LINES')# dissolve canals into single feature # intersect canals with perennial network if available if perennial_network is not None: arcpy.Intersect_analysis([perennial_network, canal_dissolve], diversion_points, "", 12, "POINT") # else intersect canals with full network minus the canals else: # select reaches that do NOT overlap with canals and save to temp_network_no_canals_shp temp_network_no_canals_lyr = arcpy.MakeFeatureLayer_management(network, 'temp_network_no_canals_lyr') arcpy.SelectLayerByLocation_management(in_layer=temp_network_no_canals_lyr, overlap_type='HAVE_THEIR_CENTER_IN', select_features=canal_dissolve, search_distance=5, selection_type='NEW_SELECTION') arcpy.SelectLayerByAttribute_management(temp_network_no_canals_lyr, 'SWITCH_SELECTION') # save temp network without canals to temp directory temp_network_no_canals_shp = os.path.join(temp_dir, "network_no_canals.shp") arcpy.CopyFeatures_management(temp_network_no_canals_lyr, temp_network_no_canals_shp) # find intersection between network without canals and canals arcpy.Intersect_analysis([temp_network_no_canals_shp, canal_dissolve], diversion_points, "", 12, "POINT") # return new diversion points shapefile return diversion_points def ipc_attributes(out_network, road, railroad, canal, valley_bottom, ownership, diversion_points, buf_30m, buf_100m, landuse, scratch, proj_path, is_verbose): """ Calculates distances from road intersections, adjacent roads, railroads and canals for each flowline segment :param out_network: The output network where fields will be added :param road: The roads shapefile for the entire area :param railroad: The railroads shapefile for the entire area :param canal: The canals shapefile for the entire area :param valley_bottom: The valley bottom shapefile for the entire area :param ownership: The land ownership shapefile for the entire area :param buf_30m: The 30m stream buffer :param buf_100m: The 100m stream buffer :param landuse: The landuse raster :param scratch: The current workspace :param projPath: The file path to the project folder :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :param perennial_network: The perennial network shapefile :return: """ # find temp directory, or make if not present temp_dir = os.path.join(proj_path, 'Temp') if not os.path.exists(temp_dir): make_temp_dir(proj_path, is_verbose) # if fields already exist, delete them fields = [f.name for f in arcpy.ListFields(out_network)] drop = ["iPC_RoadX", "iPC_Road", "iPC_RoadVB", "iPC_Rail", "iPC_RailVB", "iPC_Canal", "iPC_DivPts", "iPC_LU", "iPC_Privat", "ADMIN_AGEN"] for field in fields: if field in drop: arcpy.DeleteField_management(out_network, field) # calculate mean distance from road-stream crossings ('iPC_RoadX'), roads ('iPC_Road') and roads clipped to the valley bottom ('iPC_RoadVB') if road is not None: road_crossings = temp_dir + "\\roadx.shp" # create points at road-stream intersections arcpy.Intersect_analysis([out_network, road], road_crossings, "", "", "POINT") find_distance_from_feature(out_network, road_crossings, valley_bottom, temp_dir, buf_30m, "roadx", "iPC_RoadX", scratch, is_verbose, clip_feature = False) if road is not None: find_distance_from_feature(out_network, road, valley_bottom, temp_dir, buf_30m, "roadvb", "iPC_RoadVB", scratch, is_verbose, clip_feature = True) find_distance_from_feature(out_network, road, valley_bottom, temp_dir, buf_30m, "road", "iPC_Road", scratch, is_verbose, clip_feature = False) if railroad is not None: find_distance_from_feature(out_network, railroad, valley_bottom, temp_dir, buf_30m, "railroadvb", "iPC_RailVB", scratch, is_verbose, clip_feature = True) find_distance_from_feature(out_network, railroad, valley_bottom, temp_dir, buf_30m, "railroad", "iPC_Rail", scratch, is_verbose, clip_feature = False) if canal is not None: # find distance from canal find_distance_from_feature(out_network, canal, valley_bottom, temp_dir, buf_30m, "canal", "iPC_Canal", scratch, is_verbose, clip_feature=False) if diversion_points is not None: # calculate distance from points of diversion find_distance_from_feature(out_network, diversion_points, valley_bottom, temp_dir, buf_30m, "diversion", "iPC_DivPts", scratch, is_verbose, clip_feature = False) # assign land ownership agency to each reach if ownership is not None: if is_verbose: arcpy.AddMessage('Assigning land ownership to each reach...') spatial_join_temp = temp_dir + "\\ownership_network_join.shp" if os.path.exists(spatial_join_temp): arcpy.Delete_management(spatial_join_temp) ownership_fields = [f.name for f in arcpy.ListFields(ownership)] network_fields = [f.name for f in arcpy.ListFields(out_network)] if 'FID' in network_fields: network_fields.remove('FID') if 'Shape' in network_fields: network_fields.remove('Shape') for field in ownership_fields: if field in network_fields: arcpy.DeleteField_management(out_network, str(field)) arcpy.SpatialJoin_analysis(out_network, ownership, spatial_join_temp, 'JOIN_ONE_TO_ONE', 'KEEP_ALL', match_option = 'HAVE_THEIR_CENTER_IN') arcpy.JoinField_management(in_data=out_network, in_field='FID', join_table=spatial_join_temp, join_field='FID', fields='ADMIN_AGEN') with arcpy.da.UpdateCursor(out_network, 'ADMIN_AGEN') as cursor: for row in cursor: if row[0] == ' ': row[0] = 'None' cursor.updateRow(row) # calculate minimum distance from private or undetermined land ownership('iPC_Privat') private = temp_dir + "\\private_land.shp" private_lyr = arcpy.MakeFeatureLayer_management(ownership, "private_lyr") arcpy.SelectLayerByAttribute_management(private_lyr, 'NEW_SELECTION', """ "ADMIN_AGEN" = 'PVT' OR "ADMIN_AGEN" = 'UND' """) arcpy.CopyFeatures_management(private_lyr, private) find_distance_from_feature(out_network, private, valley_bottom, temp_dir, buf_30m, "private_land", "iPC_Privat", scratch, is_verbose, clip_feature=False) # calculate mean landuse value ('iPC_LU') if landuse is not None: add_landuse_to_table(out_network, landuse, buf_100m, scratch, is_verbose) add_min_distance(out_network) # clear the environment extent setting arcpy.ClearEnvironment("extent") def make_temp_dir(projPath, is_verbose): """ Creates a temporary directory :param projPath: The file path to the project folder :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :return: """ if is_verbose: arcpy.AddMessage("Deleting and remaking temp dir...") from shutil import rmtree temp_dir = os.path.join(projPath, 'Temp') if os.path.exists(temp_dir): rmtree(temp_dir) os.mkdir(temp_dir) def add_min_distance(out_network): """ Adds oPC_Dist field, and populated it with the minimum distance from multiple features :param out_network: The output network where fields will be added :return: """ arcpy.AddField_management(out_network, "oPC_Dist", 'DOUBLE') fields = [f.name for f in arcpy.ListFields(out_network)] all_dist_fields = ["oPC_Dist", "iPC_RoadX", "iPC_RoadVB", "iPC_RailVB", "iPC_Canal", "iPC_DivPts"] dist_fields = [] for field in all_dist_fields: if field in fields: dist_fields.append(field) with arcpy.da.UpdateCursor(out_network, dist_fields) as cursor: for row in cursor: row[0] = min(row[1:]) cursor.updateRow(row) def add_landuse_to_table(out_network, landuse, buf_100m, scratch, is_verbose): """ Adds landuse fields to the output network[iPC_LU, "iPC_VLowLU", "iPC_LowLU", "iPC_ModLU", "iPC_HighLU"] :param out_network: Output network to add fields to. :param landuse: The landuse raster. :param buf_100m: The 100m stream buffer :param scratch: Th current workspace :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :return: """ if is_verbose: arcpy.AddMessage("Calculating iPC_LU values...") arcpy.AddField_management(out_network, "iPC_LU", "DOUBLE") # create raster with just landuse code values lu_ras = Lookup(landuse, "LU_CODE") # calculate mean landuse value within 100 m buffer of each network segment zonalStatsWithinBuffer(buf_100m, lu_ras, 'MEAN', 'MEAN', out_network, "iPC_LU", scratch) # get percentage of each land use class in 100 m buffer of stream segment fields = [f.name.upper() for f in arcpy.ListFields(landuse)] if "LUI_CLASS" not in fields: arcpy.AddWarning("No field named \"LU_CLASS\" in the land use raster. Make sure that this field exists" + " with no typos if you wish to use the data from the land use raster") return buf_fields = [f.name for f in arcpy.ListFields(buf_100m)] if 'oArea' not in buf_fields: arcpy.AddField_management(buf_100m, 'oArea', 'DOUBLE') with arcpy.da.UpdateCursor(buf_100m, ['SHAPE@AREA', 'oArea']) as cursor: for row in cursor: row[1] = row[0] cursor.updateRow(row) landuse_poly = arcpy.RasterToPolygon_conversion(landuse, os.path.join(scratch, 'landuse_poly'), 'NO_SIMPLIFY', "LUI_Class") landuse_int = arcpy.Intersect_analysis([landuse_poly, buf_100m], os.path.join(scratch, 'landuse_int')) arcpy.AddField_management(landuse_int, 'propArea', 'DOUBLE') with arcpy.da.UpdateCursor(landuse_int, ['SHAPE@AREA', 'oArea', 'propArea']) as cursor: for row in cursor: row[2] = row[0]/row[1] cursor.updateRow(row) area_tbl = arcpy.Statistics_analysis(landuse_int, os.path.join(scratch, 'areaTbl'), [['propArea', 'SUM']], ['ReachID', 'LUI_CLASS']) area_piv_tbl = arcpy.PivotTable_management(area_tbl, ['ReachID'], 'LUI_CLASS', 'SUM_propArea', os.path.join(scratch, 'areaPivTbl')) sanitize_area_piv_tbl(area_piv_tbl) # create empty dictionary to hold input table field values tbl_dict = {} # add values to dictionary with arcpy.da.SearchCursor(area_piv_tbl, ['ReachID', 'VeryLow', 'Low', 'Moderate', 'High']) as cursor: for row in cursor: tbl_dict[row[0]] = [row[1], row[2], row[3], row[4]] # populate flowline network out fields arcpy.AddField_management(out_network, "iPC_VLowLU", 'DOUBLE') arcpy.AddField_management(out_network, "iPC_LowLU", 'DOUBLE') arcpy.AddField_management(out_network, "iPC_ModLU", 'DOUBLE') arcpy.AddField_management(out_network, "iPC_HighLU", 'DOUBLE') with arcpy.da.UpdateCursor(out_network, ['ReachID', 'iPC_VLowLU', 'iPC_LowLU', 'iPC_ModLU', 'iPC_HighLU']) as cursor: for row in cursor: try: aKey = row[0] row[1] = round(100*tbl_dict[aKey][0], 2) row[2] = round(100*tbl_dict[aKey][1], 2) row[3] = round(100*tbl_dict[aKey][2], 2) row[4] = round(100*tbl_dict[aKey][3], 2) cursor.updateRow(row) except: pass tbl_dict.clear() arcpy.Delete_management(lu_ras) def sanitize_area_piv_tbl(area_piv_tbl): """ Makes sure that the areaPivTbl has all the fields we need. If it doesn't, we'll add it. :param area_piv_tbl: :return: """ fields = [f.name for f in arcpy.ListFields(area_piv_tbl)] check_and_add_zero_fields(area_piv_tbl, fields, 'VeryLow') check_and_add_zero_fields(area_piv_tbl, fields, 'Low') check_and_add_zero_fields(area_piv_tbl, fields, 'Moderate') check_and_add_zero_fields(area_piv_tbl, fields, 'High') def check_and_add_zero_fields(table, fields, field_name): """ Checks that a field is in the table. If it isn't, we add it and populate it with zeros :param table: The table that we want to check :param fields: All the fields in the table (more efficient if doing multiple checks) :param field_name: The name of the field we want to check for :return: """ if field_name not in fields: arcpy.AddField_management(table, field_name, "DOUBLE") arcpy.CalculateField_management(table, field_name, 0, "PYTHON") def find_distance_from_feature(out_network, feature, valley_bottom, temp_dir, buf, temp_name, new_field_name, scratch, is_verbose, clip_feature = False): """ Finds the distance from a given feature to each stream segment and populates a new field :param out_network: The output network where new fields will be added :param feature: The feature that you want to calculate the distance from :param valley_bottom: The valley bottom shapefile :param temp_dir: The temporary folder directory :param buf: The 30m stream buffer shapefile :param temp_name: The name given to the temporary shapefile created :param new_field_name: The name of the new field to be added to the output :param scratch: The current workspace :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :param clip_feature: If true, the feature will be clipped to the valley bottom :return: """ if is_verbose: arcpy.AddMessage("Calculating " + new_field_name + " values...") arcpy.AddField_management(out_network, new_field_name, "DOUBLE") if clip_feature == True: # clip input feature to the valley bottom feature_subset = arcpy.Clip_analysis(feature, valley_bottom, os.path.join(temp_dir, temp_name + '_subset.shp')) # convert feature for count purposes feature_mts = arcpy.MultipartToSinglepart_management(feature_subset, os.path.join(temp_dir, temp_name + "_mts.shp")) else: feature_mts = arcpy.MultipartToSinglepart_management(feature, os.path.join(temp_dir, temp_name + "_mts.shp")) feature_subset = feature count = arcpy.GetCount_management(feature_mts) ct = int(count.getOutput(0)) # if there are features, then set the distance from to high value (10000 m) if ct < 1: with arcpy.da.UpdateCursor(out_network, new_field_name) as cursor: for row in cursor: row[0] = 10000.0 cursor.updateRow(row) # if there are features, calculate distance else: # set extent to the stream network arcpy.env.extent = out_network # calculate euclidean distance from input features ed_feature = EucDistance(feature_subset, cell_size = 5) # cell size of 5 m # get min distance from feature in the within 30 m buffer of each network segment if new_field_name == 'iPC_RoadX': zonalStatsWithinBuffer(buf, ed_feature, 'MINIMUM', 'MIN', out_network, new_field_name, scratch) else: zonalStatsWithinBuffer(buf, ed_feature, 'MEAN', 'MEAN', out_network, new_field_name, scratch) # delete temp fcs, tbls, etc. items = [] for item in items: arcpy.Delete_management(item) # calculate drainage area function def calc_drain_area(DEM, input_DEM): """ Calculate drainage area function :param DEM: Smoothed DEM :param input_DEM: The original input DEM :return: """ # define raster environment settings desc = arcpy.Describe(DEM) arcpy.env.extent = desc.Extent arcpy.env.outputCoordinateSystem = desc.SpatialReference arcpy.env.cellSize = desc.meanCellWidth # calculate cell area for use in drainage area calcultion height = desc.meanCellHeight width = desc.meanCellWidth cell_area = height * width # derive drainage area raster (in square km) from input DEM # note: draiange area calculation assumes input dem is in meters filled_DEM = Fill(DEM) # fill sinks in dem flow_direction = FlowDirection(filled_DEM) # calculate flow direction flow_accumulation = FlowAccumulation(flow_direction) # calculate flow accumulation drain_area = flow_accumulation * cell_area / 1000000 # calculate drainage area in square kilometers # save drainage area raster if os.path.exists(os.path.dirname(input_DEM) + "/Flow/DrainArea_sqkm.tif"): arcpy.Delete_management(os.path.dirname(input_DEM) + "/Flow/DrainArea_sqkm.tif") arcpy.CopyRaster_management(drain_area, os.path.dirname(input_DEM) + "/Flow/DrainArea_sqkm.tif") else: os.mkdir(os.path.dirname(input_DEM) + "/Flow") arcpy.CopyRaster_management(drain_area, os.path.dirname(input_DEM) + "/Flow/DrainArea_sqkm.tif") def write_xml(output_folder, coded_veg, coded_hist, seg_network, inDEM, valley_bottom, landuse, DrAr, road, railroad, canal, buf_30m, buf_100m, out_network, description): """ Writes the xml file for the project :param output_folder: The folder where everything will be output :param coded_veg: The landfire EVT layer :param coded_hist: The landfire BPS layer :param seg_network: The segmented (300m) network :param inDEM: The DEM for the entire project :param valley_bottom: The valley bottom shapefile :param landuse: The landuse raster :param DrAr: The drainage area raster :param road: The roads shapefile :param railroad: The railroads shapefile :param canal: The canals shapefile :param buf_30m: The 30m buffer shapefile :param buf_100m: The 100m buffer shapefile :param out_network: The network that new data has been added to :param description: A short description of the run that will be added to the XML :return: """ proj_path = os.path.dirname(os.path.dirname(output_folder)) output_folder_num = str(int(output_folder[-2:])) xml_file_path = proj_path + "/project.rs.xml" if not os.path.exists(xml_file_path): arcpy.AddWarning("XML file not found. Could not update XML file") return xml_file = XMLBuilder(xml_file_path) add_drain_area_to_inputs_xml(xml_file, DrAr, proj_path) realizations_element = xml_file.find("Realizations") if realizations_element is None: realizations_element = xml_file.add_sub_element(xml_file.root, "Realizations") creation_time = datetime.datetime.today().isoformat() brat_element = xml_file.add_sub_element(realizations_element, "BRAT", tags=[("dateCreated", creation_time), ("guid", getUUID()), ("id", "RZ" + output_folder_num), ("ProductVersion", "3.0.21")]) xml_file.add_sub_element(brat_element, "Name", "BRAT Realization " + output_folder_num) meta_element = xml_file.add_sub_element(brat_element, "MetaData") write_description(xml_file, meta_element, description) write_input_xml(xml_file, brat_element, proj_path, coded_veg, coded_hist, landuse, valley_bottom, road, railroad, canal, inDEM, DrAr, seg_network, buf_30m, buf_100m) write_intermediate_xml(xml_file, brat_element, proj_path, out_network) xml_file.write() def write_description(xml_file, meta_element, description): """ Writes the description to the XML :param xml_file: The XML file to write to. :param meta_element: The current meta element :param description: A short description of the run that will be added to the XML :return: """ if description is None: xml_file.add_sub_element(meta_element, "Description") elif len(description) <= 100: xml_file.add_sub_element(meta_element, "Meta", description, tags=[("name", "description")]) elif len(description) > 100: raise Exception("Description must be less than 100 characters") def write_intermediate_xml(xml_file, brat_element, proj_path, out_network): """ Writes part of the XML file for all intermediate data. :param xml_file: The XML file to write to. :param brat_element: The XML element containing BRAT data. :param proj_path: The file path to the project folder. :param out_network: The output network that new data has been added to. :return: """ intermediates_element = xml_file.add_sub_element(brat_element, "Intermediates") intermediate_element = xml_file.add_sub_element(intermediates_element, "Intermediate") xml_file.add_sub_element(intermediate_element, "Name", "BRAT Intermediate") write_xml_element_with_path(xml_file, intermediate_element, "Vector", "BRAT Input Table", out_network, proj_path) def write_input_xml(xml_file, brat_element, proj_path, coded_veg, coded_hist, landuse, valley_bottom, road, railroad, canal, inDEM, DrAr, seg_network, buf_30m, buf_100m): """ Writes part of the XML file for all input data. :param xml_file: The XML file to write to. :param brat_element: The XML element containing BRAT data. :param proj_path: The file path to the project folder. :param coded_veg: The landfire EVT layer :param coded_hist: The landfire BPS layer :param landuse: The landuse raster :param valley_bottom: The valley bottom shapefile :param road: The roads shapefile :param railroad: The railroads shapefile :param canal: The canals shapefile :param inDEM: The DEM for the entire project :param DrAr: The drainage area raster :param seg_network: The segmented (300m) network :param buf_30m: The 30m buffer shapefile :param buf_100m: The 100m buffer shapefile :return: """ inputs_element = xml_file.add_sub_element(brat_element, "Inputs") add_input_ref_element(xml_file, proj_path, inputs_element, coded_veg, "ExistingVegetation") add_input_ref_element(xml_file, proj_path, inputs_element, coded_hist, "HistoricVegetation") add_input_ref_element(xml_file, proj_path, inputs_element, landuse, "LandUse") add_input_ref_element(xml_file, proj_path, inputs_element, valley_bottom, "ValleyBottom") add_input_ref_element(xml_file, proj_path, inputs_element, road, "Roads") add_input_ref_element(xml_file, proj_path, inputs_element, railroad, "Railroads") add_input_ref_element(xml_file, proj_path, inputs_element, canal, "Canals") topo_element = xml_file.add_sub_element(inputs_element, "Topography") add_input_ref_element(xml_file, proj_path, topo_element, inDEM, "DEM") add_input_ref_element(xml_file, proj_path, topo_element, DrAr, "Flow") drain_network_element = xml_file.add_sub_element(inputs_element, "DrainageNetworks") network_element = add_input_ref_element(xml_file, proj_path, drain_network_element, seg_network, "Network") buffers_element = xml_file.add_sub_element(network_element, "Buffers") write_xml_element_with_path(xml_file, buffers_element, "Buffer", "30m Buffer", buf_30m, proj_path) write_xml_element_with_path(xml_file, buffers_element, "Buffer", "100m Buffer", buf_100m, proj_path) def add_drain_area_to_inputs_xml(xml_file, drainage_area, proj_path): """ Adds drainage area data to the input XML :param xml_file: The XML file to write to. :param drainage_area: The drainage area raster :param proj_path: The file path to the project folder. :return: """ element = xml_file.find_by_text(find_relative_path(drainage_area, proj_path)) if element is not None: # if the flow acc is already in the xml file, we don't need to do anything return inputs_element = xml_file.find("Inputs") id = find_next_available_id(xml_file, "DR") write_xml_element_with_path(xml_file, inputs_element, "Raster", "Drainage Area", drainage_area, proj_path, xml_id=id) def find_next_available_id(xml_file, id_base): """ Finds the next ID available in the current XML File :param xml_file: The XML file to write to. :param id_base: The ID to start with, then count up :return: The next available ID """ i = 1 element = xml_file.find_by_id(id_base + str(i)) while element is not None: i += 1 element = xml_file.find_by_id(id_base + str(i)) return id_base + str(i) def add_input_ref_element(xml_file, proj_path, inputs_element, input_path, new_element_name): """ :param xml_file: The XML file to write to. :param proj_path: The file path to the project folder. :param inputs_element: The element that holds Input data :param input_path: The path to the specific Input data :param new_element_name: The ame of the newly added input element. :return: """ if input_path is None: return ref_id = find_element_id_with_path(xml_file, input_path, proj_path) if ref_id is not None: return xml_file.add_sub_element(inputs_element, new_element_name, tags=[('ref', ref_id)]) else: arcpy.AddMessage(new_element_name + " could not be found in the Inputs XML, and so could not be added to the new Realization in the XML") return None def find_element_id_with_path(xml_file, path, proj_path): """ Returns the ID of the input element that has the relative path given to it :param xml_file: The XMLBuilder object :param path: The path we want to find :param proj_path: Path to the root folder of the project :return: """ relative_path = find_relative_path(path, proj_path) element = xml_file.find_by_text(relative_path) if element is not None: parent = xml_file.find_element_parent(element) return parent.attrib['id'] else: return None def validate_inputs(seg_network, road, railroad, canal, is_verbose): """ Checks if the spatial references are correct and that the inputs are what we want :param seg_network: The stream network shape file :param road: The roads shapefile :param railroad: The railroads shape file :param canal: The canals shapefile :param is_verbose: Tells us if we should print out extra debug messages :return: """ if is_verbose: arcpy.AddMessage("Validating inputs...") try: network_sr = arcpy.Describe(seg_network).spatialReference except: raise Exception("There was a problem finding the spatial reference of the stream network. " + "This is commonly caused by trying to run the Table tool directly after running the project " + "builder. Restarting ArcGIS fixes this problem most of the time.") if not network_sr.type == "Projected": raise Exception("Input stream network must have a projected coordinate system") if road is not None: if not arcpy.Describe(road).spatialReference.type == "Projected": raise Exception("Input roads must have a projected coordinate system") if railroad is not None: if not arcpy.Describe(railroad).spatialReference.type == "Projected": raise Exception("Input railroads must have a projected coordinate system") if canal is not None: if not arcpy.Describe(canal).spatialReference.type == "Projected": raise Exception("Input canals must have projected coordinate system") # --check that input network is shapefile-- if not seg_network.endswith(".shp"): raise Exception("Input network must be a shapefile (.shp)") def add_mainstem_attribute(out_network): """ Adds the mainstem attribute to our output network :param out_network: The network that we want to work with :return: None """ list_fields = arcpy.ListFields(out_network,"IsMainCh") if len(list_fields) is not 1: arcpy.AddField_management(out_network, "IsMainCh", "SHORT", "", "", "", "", "NULLABLE") arcpy.CalculateField_management(out_network,"IsMainCh",1,"PYTHON") def make_layers(out_network, diversion_pts): """ Writes the layers :param out_network: The output network, which we want to make into a layer :param diversion_pts: The diversion points shapefile :return: """ arcpy.AddMessage("Making layers...") intermediates_folder = os.path.dirname(out_network) buffers_folder = os.path.join(intermediates_folder, "01_Buffers") topo_folder = make_folder(intermediates_folder, "02_TopographicMetrics") anthropogenic_metrics_folder = make_folder(intermediates_folder, "03_AnthropogenicMetrics") perennial_folder = make_folder(intermediates_folder, "04_Perennial") trib_code_folder = os.path.dirname(os.path.abspath(__file__)) symbology_folder = os.path.join(trib_code_folder, 'BRATSymbology') dist_to_infrastructure_symbology = os.path.join(symbology_folder, "DistancetoClosestInfrastructure.lyr") dist_to_road_in_valley_bottom_symbology = os.path.join(symbology_folder, "DistancetoRoadinValleyBottom.lyr") dist_to_road_crossing_symbology = os.path.join(symbology_folder, "DistancetoRoadCrossing.lyr") dist_to_road_symbology = os.path.join(symbology_folder, "DistancetoRoad.lyr") dist_to_railroad_in_valley_bottom_symbology = os.path.join(symbology_folder, "DistancetoRailroadinValleyBottom.lyr") dist_to_railroad_symbology = os.path.join(symbology_folder, "DistancetoRailroad.lyr") dist_to_canal_symbology = os.path.join(symbology_folder, "DistancetoCanal.lyr") pts_diversion_symbology = os.path.join(symbology_folder, "PointsofDiversion.lyr") dist_to_pts_diversion_symbology = os.path.join(symbology_folder, "DistancetoPointsofDiversion.lyr") land_use_symbology = os.path.join(symbology_folder, "LandUseIntensity.lyr") land_ownership_per_reach_symbology = os.path.join(symbology_folder, "LandOwnershipperReach.lyr") priority_translocations_symbology = os.path.join(symbology_folder, "PriorityBeaverTranslocationAreas.lyr") slope_symbology = os.path.join(symbology_folder, "ReachSlope.lyr") drain_area_symbology = os.path.join(symbology_folder, "UpstreamDrainageArea.lyr") buffer_30m_symbology = os.path.join(symbology_folder, "buffer_30m.lyr") buffer_100m_symbology = os.path.join(symbology_folder, "buffer_100m.lyr") perennial_symbology = os.path.join(symbology_folder, "Perennial.lyr") make_buffer_layers(buffers_folder, buffer_30m_symbology, buffer_100m_symbology) make_layer(topo_folder, out_network, "Reach Slope", slope_symbology, is_raster=False) make_layer(topo_folder, out_network, "Upstream Drainage Area", drain_area_symbology, is_raster=False) if diversion_pts is not None: try: make_layer(os.path.dirname(diversion_pts), diversion_pts, "Provisional Points of Diversion", pts_diversion_symbology, is_raster=False) except Exception as err: print(err) fields = [f.name for f in arcpy.ListFields(out_network)] if 'iPC_LU' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Land Use Intensity", land_use_symbology, is_raster=False) if 'iPC_RoadX' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Distance to Road Crossing", dist_to_road_crossing_symbology, is_raster=False, symbology_field ='iPC_RoadX') if 'iPC_Road' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Distance to Road", dist_to_road_symbology, is_raster=False, symbology_field ='iPC_Road') if 'iPC_RoadVB' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Distance to Road in Valley Bottom", dist_to_road_in_valley_bottom_symbology, is_raster=False, symbology_field ='iPC_RoadVB') if 'iPC_Rail' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Distance to Railroad", dist_to_railroad_symbology, is_raster=False, symbology_field ='iPC_Rail') if 'iPC_RailVB' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Distance to Railroad in Valley Bottom", dist_to_railroad_in_valley_bottom_symbology, is_raster=False, symbology_field ='iPC_RailVB') if 'iPC_Canal' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Distance to Canal", dist_to_canal_symbology, is_raster=False, symbology_field ='iPC_Canal') if 'iPC_DivPts' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Distance to Points of Diversion", dist_to_pts_diversion_symbology, is_raster=False, symbology_field='iPC_DivPts') if 'ADMIN_AGEN' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Land Ownership per Reach", land_ownership_per_reach_symbology, is_raster=False, symbology_field='ADMIN_AGEN') if 'iPC_Privat' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Priority Beaver Translocation Areas", priority_translocations_symbology, is_raster=False, symbology_field='iPC_Privat') if 'oPC_Dist' in fields: make_layer(anthropogenic_metrics_folder, out_network, "Distance to Closest Infrastructure", dist_to_infrastructure_symbology, is_raster=False, symbology_field ='oPC_Dist') if 'IsPeren' in fields: make_layer(perennial_folder, out_network, "Perennial", perennial_symbology, is_raster=False, symbology_field="IsPeren") def handle_braids(seg_network_copy, canal, proj_path, find_clusters, perennial_network, is_verbose): """ Finds multi-threaded attributes in the network, and optionally finds clusters. :param seg_network_copy: The network where data will be updated :param canal: The canals shapefile :param proj_path: The file path to the project folder. :param find_clusters: If true, clusters will be found via the braid handler script :param perennial_network: The perennial network shapefile :param is_verbose: If true, this option enables ArcMap to provide messages for each step conducted by the tool. :return: """ if is_verbose: arcpy.AddMessage("Finding multi-threaded attributes...") add_mainstem_attribute(seg_network_copy) # find braided reaches temp_dir = os.path.join(proj_path, 'Temp') if not os.path.exists(temp_dir): os.mkdir(temp_dir) FindBraidedNetwork.main(seg_network_copy, canal, temp_dir, perennial_network, is_verbose) if find_clusters: arcpy.AddMessage("Finding Clusters...") clusters = BRAT_Braid_Handler.find_clusters(seg_network_copy) BRAT_Braid_Handler.add_cluster_id(seg_network_copy, clusters) # if 'StreamName' is field then run the update_multiCh function fields = [f.name for f in arcpy.ListFields(seg_network_copy)] if 'StreamName' in fields: BRAT_Braid_Handler.update_multi_ch(seg_network_copy) def make_buffer_layers(buffers_folder, buffer_30m_symbology, buffer_100m_symbology): """ Makes a layer for each buffer :param buffers_folder: The path to the buffers folder :return: Nothing """ for file_name in os.listdir(buffers_folder): if file_name.endswith(".shp"): file_path = os.path.join(buffers_folder, file_name) given_symbology = None if "30m" in file_name: new_layer_name = "30 m Buffer" given_symbology = buffer_30m_symbology elif "100m" in file_name: new_layer_name = "100 m Buffer" given_symbology = buffer_100m_symbology make_layer(buffers_folder, file_path, new_layer_name, given_symbology) def parse_input_bool(given_input): """ Takes an ArcMap bool input and coverts it to be usable in Python :param given_input: The given ArcMap input :return: Converted Bool """ if given_input == 'false' or given_input is None: return False else: return True def delete_with_arcpy(stuffToDelete): """ Deletes everything in a list with arcpy.Delete_management() :param stuffToDelete: A list of stuff to delete :return: """ for thing in stuffToDelete: arcpy.Delete_management(thing) def run_tests(seg_network_copy, is_verbose): """ Runs tests on the tool's output :param seg_network_copy: The network that we want to test :return: """ if is_verbose: arcpy.AddMessage("Running tests...") run_tests = True if not run_tests: # don't run tests in execution return from Tests import test_reach_id_is_unique, report_exceptions, TestException test_exceptions = [] try: test_reach_id_is_unique(seg_network_copy) except TestException as e: test_exceptions.append(str(e)) report_exceptions(test_exceptions) if __name__ == '__main__': main( sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7], sys.argv[8], sys.argv[9], sys.argv[10], sys.argv[11], sys.argv[12], sys.argv[13], sys.argv[14], sys.argv[15], sys.argv[16], sys.argv[17], sys.argv[18], sys.argv[19])