# Wave Attenuation Toolbox (WATTE) # Copyright (C) 2020 M. Foster-Martinez (madeline@berkeley.edu) # Developed at Louisiana State University (LSU) # Collaborators: Scott Hagen, LSU; Karim Alizad, University of South Carolina # Utilizes code from http://gis4geomorphology.com/stream-transects-partial/ and http://nodedangles.wordpress.com/2011/05/01/quick-dirty-arcpy-batch-splitting-polylines-to-a-specific-length/ # Version 2.0 # WATTE is an ArcGIS toolbox. It estimates and maps wave attenuation along coastlines following an exponential decay. # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # You should have received a copy of the GNU General Public License # along with this program (LICENSE.txt file). If not, see . ########### import arcpy from arcpy import env, Raster from arcpy.sa import * import math, numpy arcpy.env.overwriteOutput = True arcpy.env.XYResolution = "0.00001 Meters" arcpy.env.XYTolerance = "0.0001 Meters" class Toolbox(object): def __init__(self): """Define the toolbox (the name of the toolbox is the name of the .pyt file).""" self.label = "Wave Attenuation Mapping" self.alias = "Wave Attenuation Mapping" # List of tool classes associated with this toolbox self.tools = [WA_map] class WA_map(object): def __init__(self): """Define the tool (tool name is the name of the class).""" self.label = "WATTE" self.description = "Mapping Wave Attenuation with a GIS Toolbox" self.canRunInBackground = False #Declare all parameters: def getParameterInfo(self): """Define parameter definitions""" param0 = arcpy.Parameter() param0.name = "Parameter0" param0.displayName = "Input Raster Dataset:" param0.parameterType = "Required" param0.direction = "Input" param0.datatype = "DERasterDataset" param1 = arcpy.Parameter() param1.name = "Parameter1" param1.displayName = "Classification(s), Water:" param1.parameterType = "Required" param1.direction = "Input" param1.datatype = "String" param2= arcpy.Parameter() param2.name = "Parameter2" param2.displayName = "Classification(s), Other:" param2.parameterType = "Required" param2.direction = "Input" param2.datatype = "String" param3 = arcpy.Parameter() param3.name = "Parameter3" param3.displayName = "Classification(s), Marsh:" param3.parameterType = "Required" param3.direction = "Input" param3.datatype = "String" param4 = arcpy.Parameter() param4.name = "Parameter4" param4.displayName = "Decay constant for each classification of marsh:" param4.parameterType = "Required" param4.direction = "Input" param4.datatype = "String" param5 = arcpy.Parameter() param5.name = "Parameter5" param5.displayName = "Transect length:" param5.parameterType = "Required" param5.direction = "Input" param5.datatype = "Long" param6 = arcpy.Parameter() param6.name = "Parameter6" param6.displayName = "Distance between transects:" param6.parameterType = "Required" param6.direction = "Input" param6.datatype = "GPDouble" param7 = arcpy.Parameter() param7.name = "Parameter7" param7.displayName = "Wave decay stopping point:" param7.parameterType = "Required" param7.direction = "Input" param7.datatype = "GPDouble" param8 = arcpy.Parameter() param8.name = "Parameter8" param8.displayName = "Spacing of points along transect:" param8.parameterType = "Required" param8.direction = "Input" param8.datatype = "GPDouble" param9 = arcpy.Parameter() param9.name = "Parameter9" param9.displayName = "Interpolation method:" param9.parameterType = "Required" param9.direction = "Input" param9.datatype = "String" return [param0,param1,param2,param3, param4, param5, param6, param7, param8, param9] #Set defaults for parameter values def updateParameters(self, parameters): """Modify the values and properties of parameters before internal validation is performed. This method is called whenever a parameter has been changed.""" if parameters[0].value: if not parameters[7].altered: parameters[7].values = "0.02" if not parameters[6].altered: parameters[6].values = "10" if not parameters[5].altered: parameters[5].values = "200" if not parameters[8].altered: parameters[8].values = "5" if not parameters[9].altered: parameters[9].values = "IDW" return def execute(self, parameters, messages): """The source code of the tool.""" input_RD = parameters[0].valueAsText # RD = Raster Dataset #Parameters is a list, so this takes the first item in that list messages.addMessage("INPUT RASTER DATASET = " + input_RD) water_class = parameters[1].valueAsText messages.addMessage("Water classification = " + water_class) water_class = water_class.split(",") #the classifications for "other" are referred to as "upland" in the code upland_class = parameters[2].valueAsText upland_class = upland_class.split(",") marsh_class = parameters[3].valueAsText marsh_class = marsh_class.split(",") marsh_k = parameters[4].valueAsText messages.addMessage("Marsh decay constant(s) = " + marsh_k) marsh_k = marsh_k.split(",") TransecLength = parameters[5].valueAsText TransecLength_Unit = "METERS" DistanceSplit = parameters[6].valueAsText stop_pt = parameters[7].valueAsText point_space = parameters[8].valueAsText messages.addMessage("Stopping point = " + stop_pt) interp_method = parameters[9].valueAsText messages.addMessage("WATTE Copyright (C) 2020 M. Foster-Martinez; This program comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under certain conditions. See Licence file for further details.") import os #Creating the files that will be called/filled later path = arcpy.Describe(input_RD).path #OutDir = env.workspace # Create "WATTE" file geodatabase WorkFolder = path WATTE_GDB=WorkFolder+"\WATTE_GDB.gdb" arcpy.CreateFileGDB_management(WorkFolder, "WATTE_GDB", "CURRENT") env.workspace=WATTE_GDB OutDir = env.workspace InterpArea = OutDir+"/"+"interp_IA" name = arcpy.Describe(input_RD).basename + "_rlayer" output_RL = os.path.join(path,name) # RL = Raster Layer messages.addMessage("output raster layer is " +output_RL) PW = OutDir+"/"+"pw" # PW = Polygon of Water PM = OutDir+"/"+ "pm" # PM = Polygon of Marsh PU = OutDir+"/"+"pu" # PU = Polygon of Upland WaterEdge = OutDir+"/"+"wel" # wel = water edge line GTP = OutDir+"/"+"GTP" #global transect polygon GWT = OutDir+"/"+"GWT" # GWT = global wave transmission GWT_NT = OutDir+"/"+"GWT_NearTable" # GWT = global wave transmission name = "GWT_FL" GWT_FL = os.path.join(path,name) # GWT = global wave transmission, FL = feature layer FWT = OutDir+"/"+ "WaveTranPts" # FWT = final wave transmission FWT_FL = OutDir+"/"+"FWT_FL" # FWT = final wave transmission feature layer name = "Complete_Attenuation_Line" # CAL = complete attenuation line CAL = os.path.join(path,name) name = "GTP_FL" GTP_FL = os.path.join(path,name) name = arcpy.Describe(input_RD).basename + "_OT" # OT = out transects OutputTransect = os.path.join(path,name) name = arcpy.Describe(input_RD).basename + "_OT_FL" # OT = out transects feature layer OutputTransectFL = os.path.join(path,name) TempTransect = OutDir+"/"+"TempT" # TempT = temporary transect DPM = OutDir+"/"+"DPM" #dissolved polygon of marsh MarshEdge = OutDir+"/"+"ME" # ME = Marsh Edge WMEdge = OutDir+"/"+"WME" # WME = WAter-Marsh Edge InterpAreaUW = OutDir+"/"+"IA_UW" # _IA_UW = Interpolation Area with upland and water InterpAreaU = OutDir+"/"+ "IA_U"# _IA_U = Interpolation Area with upland InterpArea = OutDir+"/"+"interp_IA" # IA = Interpolation Area name = "WvHt_Trans" Interp_WT = os.path.join(path,name) Interp_WT_0 = OutDir+"/"+ "Interp_WT_0" # interpolated raster multipled by 0 Interp_WT_0_int = OutDir+"/"+ "Itp_WT0int" # that 0 raster converted to integers Interp_WT_0_poly = OutDir+"/"+ "Inp_WT0p"# that integer raster converted to a polygon Interp_WT_0_line = OutDir+"/"+ "Inp_WT0l" # that polygon converted to a line messages.addMessage("Making Raster Layer...") arcpy.MakeRasterLayer_management(input_RD, output_RL) messages.addMessage("Selecting the water portion of the raster") arcpy.SelectLayerByAttribute_management(output_RL, "NEW_SELECTION",'"VALUE" = ' + water_class[0]) if len(water_class) > 1: for water in water_class: arcpy.SelectLayerByAttribute_management(output_RL, "ADD_TO_SELECTION",'"VALUE" = ' + water) messages.addMessage("Converting to a polygon of only water..." ) arcpy.RasterToPolygon_conversion(output_RL, PW,"SIMPLIFY","VALUE","SINGLE_OUTER_PART") messages.addMessage("Converting the water polygon to a polyline of the water's edge..." ) arcpy.PolygonToLine_management(PW, WaterEdge, "IDENTIFY_NEIGHBORS") messages.addMessage("Selecting the 'other' portion of the raster") arcpy.SelectLayerByAttribute_management(output_RL, "NEW_SELECTION",'"VALUE" = ' + upland_class[0]) if len(upland_class) > 1: for upland in upland_class: arcpy.SelectLayerByAttribute_management(output_RL, "ADD_TO_SELECTION",'"VALUE" = ' + upland) messages.addMessage("Converting to a polygon of only 'other'..." ) arcpy.RasterToPolygon_conversion(output_RL, PU,"SIMPLIFY","VALUE","SINGLE_OUTER_PART") messages.addMessage("Creating a polygon of the marsh areas....") arcpy.SelectLayerByAttribute_management(output_RL, "NEW_SELECTION",'"VALUE" = ' + marsh_class[0]) if len(marsh_class) > 1: for marsh in marsh_class: arcpy.SelectLayerByAttribute_management(output_RL, "ADD_TO_SELECTION",'"VALUE" = ' + marsh) arcpy.RasterToPolygon_conversion(output_RL, PM,"SIMPLIFY","VALUE","SINGLE_OUTER_PART") arcpy.Dissolve_management(PM, DPM,"", "", "MULTI_PART", "DISSOLVE_LINES") arcpy.PolygonToLine_management( DPM, MarshEdge, "IDENTIFY_NEIGHBORS") messages.addMessage("Creating a line to denote the edge of the marsh...") arcpy.Intersect_analysis([MarshEdge,WaterEdge], WMEdge, "ALL", "8 Meters", "INPUT") messages.addMessage("Creating the mask for interpolating the wave decay values...") arcpy.Buffer_analysis(WMEdge, InterpAreaUW, TransecLength, "FULL", "","ALL", "", "PLANAR") arcpy.Clip_analysis(InterpAreaUW, PM,InterpArea,"") Lines = (WMEdge) SplitType="Split at approximate distance" env.workspace = path messages.addMessage("Drawing the transects " + (TransecLength) + "m long every " + (DistanceSplit) + "m along the water's edge..." ) ####################################### Transect Tool Code Source:http://gis4geomorphology.com/stream-transects-partial/ ######################################### # Def splitline module ###START SPLIT LINE CODE IN A SAME DISTANCE### Source: http://nodedangles.wordpress.com/2011/05/01/quick-dirty-arcpy-batch-splitting-polylines-to-a-specific-length/ def splitline (inFC,FCName,alongDist): OutDir = env.workspace outFCName = FCName outFC = OutDir+"/"+outFCName def distPoint(p1, p2): #returns the distance between the two points p1 and p2, which are given in xy coordinates calc1 = p1.X - p2.X calc2 = p1.Y - p2.Y return math.sqrt((calc1**2)+(calc2**2)) def midpoint(prevpoint,nextpoint,targetDist,totalDist): #returns a single point object that is some distance from the previous point (set by the ratio target distance/total distance) newX = prevpoint.X + ((nextpoint.X - prevpoint.X) * (targetDist/totalDist)) newY = prevpoint.Y + ((nextpoint.Y - prevpoint.Y) * (targetDist/totalDist)) return arcpy.Point(newX, newY) def splitShape(feat,splitDist): # Count the number of points in the current multipart feature # partcount = feat.partCount partnum = 0 # Enter while loop for each part in the feature (if a singlepart feature # this will occur only once) # lineArray = arcpy.Array() while partnum < partcount: # Print the part number # #print "Part " + str(partnum) + ":" part = feat.getPart(partnum) #print part.count totalDist = 0 pnt = part.next() pntcount = 0 prevpoint = None shapelist = [] # Enter while loop for each vertex # while pnt: if not (prevpoint is None): thisDist = distPoint(prevpoint,pnt) maxAdditionalDist = splitDist - totalDist print thisDist, totalDist, maxAdditionalDist if (totalDist+thisDist)> splitDist: while(totalDist+thisDist) > splitDist: maxAdditionalDist = splitDist - totalDist #print thisDist, totalDist, maxAdditionalDist newpoint = midpoint(prevpoint,pnt,maxAdditionalDist,thisDist) lineArray.add(newpoint) shapelist.append(lineArray) lineArray = arcpy.Array() lineArray.add(newpoint) prevpoint = newpoint thisDist = distPoint(prevpoint,pnt) totalDist = 0 lineArray.add(pnt) totalDist+=thisDist else: totalDist+=thisDist lineArray.add(pnt) #shapelist.append(lineArray) else: lineArray.add(pnt) totalDist = 0 prevpoint = pnt pntcount += 1 pnt = part.next() # If pnt is null, either the part is finished or there is an # interior ring # if not pnt: pnt = part.next() if pnt: print "Interior Ring:" partnum += 1 if (lineArray.count > 1): shapelist.append(lineArray) return shapelist if arcpy.Exists(outFC): arcpy.Delete_management(outFC) arcpy.Copy_management(inFC,outFC) #origDesc = arcpy.Describe(inFC) #sR = origDesc.spatialReference #revDesc = arcpy.Describe(outFC) #revDesc.ShapeFieldName deleterows = arcpy.UpdateCursor(outFC) for iDRow in deleterows: deleterows.deleteRow(iDRow) try: del iDRow del deleterows except: pass inputRows = arcpy.SearchCursor(inFC) outputRows = arcpy.InsertCursor(outFC) fields = arcpy.ListFields(inFC) numRecords = int(arcpy.GetCount_management(inFC).getOutput(0)) OnePercentThreshold = numRecords // 100 #printit(numRecords) iCounter = 0 iCounter2 = 0 for iInRow in inputRows: inGeom = iInRow.shape iCounter+=1 iCounter2+=1 if (iCounter2 > (OnePercentThreshold+0)): #printit("Processing Record "+str(iCounter) + " of "+ str(numRecords)) iCounter2=0 if (inGeom.length > alongDist): shapeList = splitShape(iInRow.shape,alongDist) for itmp in shapeList: newRow = outputRows.newRow() for ifield in fields: if (ifield.editable): newRow.setValue(ifield.name,iInRow.getValue(ifield.name)) newRow.shape = itmp outputRows.insertRow(newRow) else: outputRows.insertRow(iInRow) del inputRows del outputRows #printit("Done!") ###END SPLIT LINE CODE IN A SAME DISTANCE### # Create "General" file geodatabase WorkFolder=env.workspace General_GDB=WorkFolder+"\General.gdb" arcpy.CreateFileGDB_management(WorkFolder, "General", "CURRENT") env.workspace=General_GDB #Unsplit Line LineDissolve="LineDissolve" arcpy.Dissolve_management (Lines, LineDissolve,"", "", "SINGLE_PART") LineSplit="LineSplit" #Split Line if SplitType=="Split at approximate distance": splitline(LineDissolve, LineSplit, float(DistanceSplit)) else: arcpy.SplitLine_management (LineDissolve, LineSplit) #Add fields to LineSplit FieldsNames=["LineID", "Direction", "Azimuth", "X_mid", "Y_mid", "AziLine_1", "AziLine_2", "Distance"] for fn in FieldsNames: arcpy.AddField_management (LineSplit, fn, "DOUBLE") #Calculate Fields CodeBlock_Direction="""def GetAzimuthPolyline(shape): radian = math.atan((shape.lastpoint.x - shape.firstpoint.x)/(shape.lastpoint.y - shape.firstpoint.y)) degrees = radian * 180 / math.pi return degrees""" CodeBlock_Azimuth="""def Azimuth(direction): if direction < 0: azimuth = direction + 360 return azimuth else: return direction""" CodeBlock_NULLS="""def findNulls(fieldValue): if fieldValue is None: return 0 elif fieldValue is not None: return fieldValue""" arcpy.CalculateField_management (LineSplit, "LineID", "!OBJECTID!", "PYTHON_9.3") arcpy.CalculateField_management (LineSplit, "Direction", "GetAzimuthPolyline(!Shape!)", "PYTHON_9.3", CodeBlock_Direction) arcpy.CalculateField_management (LineSplit, "Direction", "findNulls(!Direction!)", "PYTHON_9.3", CodeBlock_NULLS) arcpy.CalculateField_management (LineSplit, "Azimuth", "Azimuth(!Direction!)", "PYTHON_9.3", CodeBlock_Azimuth) arcpy.CalculateField_management (LineSplit, "X_mid", "!Shape!.positionAlongLine(0.5,True).firstPoint.X", "PYTHON_9.3") arcpy.CalculateField_management (LineSplit, "Y_mid", "!Shape!.positionAlongLine(0.5,True).firstPoint.Y", "PYTHON_9.3") CodeBlock_AziLine1="""def Azline1(azimuth): az1 = azimuth + 90 if az1 > 360: az1-=360 return az1 else: return az1""" CodeBlock_AziLine2="""def Azline2(azimuth): az2 = azimuth - 90 if az2 < 0: az2+=360 return az2 else: return az2""" arcpy.CalculateField_management (LineSplit, "AziLine_1", "Azline1(!Azimuth!)", "PYTHON_9.3", CodeBlock_AziLine1) arcpy.CalculateField_management (LineSplit, "AziLine_2", "Azline2(!Azimuth!)", "PYTHON_9.3", CodeBlock_AziLine2) arcpy.CalculateField_management (LineSplit, "Distance", TransecLength, "PYTHON_9.3") #Generate Azline1 and Azline2 spatial_reference=arcpy.Describe(Lines).spatialReference Azline1="Azline1" Azline2="Azline2" arcpy.BearingDistanceToLine_management (LineSplit, Azline1, "X_mid", "Y_mid", "Distance", TransecLength_Unit, "AziLine_1", "DEGREES", "GEODESIC", "LineID", spatial_reference) arcpy.BearingDistanceToLine_management (LineSplit, Azline2, "X_mid", "Y_mid", "Distance", TransecLength_Unit, "AziLine_2", "DEGREES", "GEODESIC", "LineID", spatial_reference) #Create Azline and append Azline1 and Azline2 Azline="Azline" arcpy.CreateFeatureclass_management(env.workspace, "Azline", "POLYLINE", "", "", "", spatial_reference) arcpy.AddField_management (Azline, "LineID", "DOUBLE") arcpy.Append_management ([Azline1, Azline2], Azline, "NO_TEST") #Dissolve Azline Azline_Dissolve="Azline_Dissolve" arcpy.Dissolve_management (Azline, Azline_Dissolve,"LineID", "", "SINGLE_PART") #Add Fields to Azline_Dissolve FieldsNames2=["x_start", "y_start", "x_end", "y_end"] for fn2 in FieldsNames2: arcpy.AddField_management (Azline_Dissolve, fn2, "DOUBLE") #Calculate Azline_Dissolve fields arcpy.CalculateField_management (Azline_Dissolve, "x_start", "!Shape!.positionAlongLine(0,True).firstPoint.X", "PYTHON_9.3") arcpy.CalculateField_management (Azline_Dissolve, "y_start", "!Shape!.positionAlongLine(0,True).firstPoint.Y", "PYTHON_9.3") arcpy.CalculateField_management (Azline_Dissolve, "x_end", "!Shape!.positionAlongLine(1,True).firstPoint.X", "PYTHON_9.3") arcpy.CalculateField_management (Azline_Dissolve, "y_end", "!Shape!.positionAlongLine(1,True).firstPoint.Y", "PYTHON_9.3") #Generate output file arcpy.XYToLine_management (Azline_Dissolve, OutputTransect,"x_start", "y_start", "x_end","y_end", "", "", spatial_reference) #Delete General.gdb arcpy.Delete_management(General_GDB) ################end transect tool code############################################################################################################################ total_n = arcpy.GetCount_management(OutputTransect+".shp") arcpy.MakeFeatureLayer_management(OutputTransect+".shp", OutputTransectFL) mid = (int(float(TransecLength)))/int(float(point_space)) n = 0 while n < int(float(total_n.getOutput(0))): messages.addMessage("Currently on transect number " + str(n+1) + " of " + str(total_n)) arcpy.SelectLayerByAttribute_management(OutputTransectFL, "NEW_SELECTION", "FID = " + str(n)) arcpy.GeneratePointsAlongLines_management(OutputTransectFL, TempTransect, "DISTANCE", point_space +" Meters") arcpy.AddXY_management(TempTransect) #Make a list of x and y coordinates x = "POINT_X" y = "POINT_Y" X = [] Y = [] cursor = arcpy.SearchCursor(TempTransect) for row in cursor: X.append(row.getValue(x)) Y.append(row.getValue(y)) #Make a list of the raster values for each of those coordinates RV = [] i = 0 while i < len(X): #Note: ExtractMultiValuesToPoints was attempted to increase speed but did not work well (created complications with point orders) V = arcpy.GetCellValue_management(input_RD, str(X[i]) + " " + str(Y[i])) RV.append(V) i = i + 1 #Now to the good stuff...wave transmission WT = [9999] * len(RV) LO = [9999] * len(RV) if water_class.count(str(RV[mid+2])) == 1 and water_class.count(str(RV[mid-2])) != 1 :#figuring out how the transect is oriented type = "countdown" else: type = "countup" WT_previous = 1 c = mid-1 lo = 1 while c < ((mid*2)-2) and c >= 0 and WT_previous>= float(stop_pt): #for pt in RV: if water_class.count(str(RV[c])) == 1 and 'm' in locals(): WT_previous = 1 c = mid*4 #if it runs into water after going over marsh, it should stop. Increasing the value of c will kick it out of the loop. elif water_class.count(str(RV[c])) == 1: WT_previous = 1 elif upland_class.count(str(RV[c])) == 1: WT_previous = 1 c = mid*4 #if it has run into upland, it should just be kicked out of the loop, which is what this line essentially does elif str(RV[c]) == "NoData" or str(RV[c]) == "None": WT_previous = 1 else: if marsh_class.count(str(RV[c])) == 1: d = marsh_class.index(str(RV[c])) k = float(marsh_k[d]) m = 10 #trip so that if m does not exist then we know it's a transect we don't care about else: messages.addMessage("This class has not been assigned. The value of this point is " + str(RV[c])) if WT[c-1] == 9999 and WT[c+1] == 9999: wave_transmission = 1 else: wave_transmission = numpy.exp(-k*int(float(point_space))) total_wave_transmission = WT_previous*wave_transmission WT.pop(c) total_wave_transmission_percent = total_wave_transmission*100 WT.insert(c,total_wave_transmission_percent) WT_previous = total_wave_transmission LO.pop(c) LO.insert(c,lo) if type == "countdown": c = c -1 lo = lo+1 else: c = c+1 lo = lo+1 if 'm' in locals(): arcpy.AddField_management(TempTransect,"Wave_Trans","FLOAT") with arcpy.da.UpdateCursor(TempTransect,"Wave_Trans") as cursor: j = 0 for row in cursor: row[0]=WT[j] cursor.updateRow(row) j = j+1 arcpy.AddField_management(TempTransect,"LineOrder","SHORT") with arcpy.da.UpdateCursor(TempTransect,"LineOrder") as cursor: j = 0 for row in cursor: row[0]=LO[j] cursor.updateRow(row) j = j+1 #if it's the first transect, make the global point file. if it's not, add the transect points to the global point file. if 'mr' in locals(): arcpy.Append_management(TempTransect,GTP, "TEST","","") else: arcpy.Merge_management(TempTransect, GTP) mr = [] n = n+1 del(m) else: n = n+1 arcpy.Delete_management(OutputTransect+'.shp',"") messages.Addmessage("Getting rid of points with no value...") arcpy.MakeFeatureLayer_management(GTP,GTP_FL) arcpy.SelectLayerByAttribute_management(GTP_FL+".shp" ,"NEW_SELECTION", '"Wave_Trans" < 9999') arcpy.CopyFeatures_management(GTP_FL, GWT) ### messages.Addmessage("Checking License type") if arcpy.ProductInfo() == "ArcInfo": messages.Addmessage("Advanced License found. Getting rid of points that are too close...") spacing = int(float(point_space))/2 arcpy.GenerateNearTable_analysis(GWT, GWT, GWT_NT, str(spacing)+" Meters", "NO_LOCATION", "NO_ANGLE", "ALL", "0", "PLANAR") LO_Values = [r[0] for r in arcpy.da.SearchCursor(GWT, ["LineOrder"])] Overlap = [1] * len(LO_Values) cursor = arcpy.SearchCursor(GWT_NT) for row in cursor: point1 = (row.getValue("IN_FID")) point2 = (row.getValue("NEAR_FID")) if LO_Values[point1] <= LO_Values[point2]: Overlap[point2] = 0 else: Overlap[point1] = 0 arcpy.AddField_management(GWT + ".shp","Overlap","SHORT") with arcpy.da.UpdateCursor(GWT,"Overlap") as cursor: j = 0 for row in cursor: row[0]=Overlap[j] cursor.updateRow(row) j = j+1 arcpy.MakeFeatureLayer_management(GWT,GWT_FL) arcpy.SelectLayerByAttribute_management(GWT_FL,"NEW_SELECTION", '"Overlap" > 0') arcpy.CopyFeatures_management(GWT_FL, FWT) else: messages.Addmessage("Advanced License not found, moving to next step") arcpy.MakeFeatureLayer_management(GWT,GWT_FL) arcpy.CopyFeatures_management(GWT_FL, FWT) messages.Addmessage("Interpolating points to create a raster...") arcpy.env.mask = InterpArea + ".shp" #set the mask to be only where we want to interpolate if interp_method == "IDW": arcpy.gp.Idw_sa(FWT + ".shp", "Wave_Trans", Interp_WT , "10", "2", "VARIABLE 4 20", "") elif interp_method == "Kriging": arcpy.gp.Kriging_sa(FWT + ".shp", "Wave_Trans", Interp_WT, "Spherical 1.541455", "1.54145497959852", "VARIABLE 4 20", "") if arcpy.ProductInfo() == "ArcInfo": messages.Addmessage("Advanced license found. Creating the complete attenuation line...") #multiply that raster by 0 so it's all one thing outRas = Raster(Interp_WT)*0 #turn that into an integer (seems like that would be unecessary, but it doesn't work otherwise) arcpy.gp.Int_sa(outRas, Interp_WT_0_int) #convert that raster into a Polygon arcpy.RasterToPolygon_conversion(Interp_WT_0_int, Interp_WT_0_poly, "NO_SIMPLIFY", "Value", "SINGLE_OUTER_PART", "") #convert that polygon into a line arcpy.PolygonToLine_management(Interp_WT_0_poly, Interp_WT_0_line, "IDENTIFY_NEIGHBORS") #erase all the places where the line intersects with the marsh edge line arcpy.Erase_analysis(Interp_WT_0_line,MarshEdge, CAL, "20 Meters") else: messages.Addmessage("Advanced license not found. Complete attenuation line will need to be created manually. Please see WATTE instructions.") #multiply that raster by 0 so it's all one thing outRas = Raster(Interp_WT)*0 #turn that into an integer (seems like that would be unecessary, but it doesn't work otherwise) arcpy.gp.Int_sa(outRas, Interp_WT_0_int) #convert that raster into a Polygon arcpy.RasterToPolygon_conversion(Interp_WT_0_int, Interp_WT_0_poly, "NO_SIMPLIFY", "Value", "SINGLE_OUTER_PART", "") #convert that polygon into a line arcpy.PolygonToLine_management(Interp_WT_0_poly, Interp_WT_0_line, "IDENTIFY_NEIGHBORS") #erase all the places where the line intersects with the marsh edge line #arcpy.Delete_management(WATTE_GDB) #add this line back if you do not want to output the intermediate steps