#!/usr/bin/env python # -*- coding: utf-8 -*- # # Written By: A. Layman # ArcGIS for Desktop + Spatial Analyst Extension required # More information: DOI LINK HERE # # # Import modules import arcpy from arcpy import env from arcpy.sa import * import string, os, sys, time, array import re from math import * import csv import fnmatch import numpy # Check out any necessary licenses arcpy.CheckOutExtension("Spatial") # Print system start time and date tstart = time.strftime('%X %x %Z') print (tstart) ################################################################################### # Set Global Variables # Fetch Raster Weighting Methods # EF = multiplies directional fetch rasters by % frequency # that wind blows in that direction, then sums all rasters # # REI = Relative Exposure Index, after Keddy (1982). Multiplies directional fetch rasters # by both % frequency that wind blows along the direction and average wind velocity # in that direction, then sums all rasters # Input text file ontaining direction wind-weighting # formatted according to manuscript wind_text = r"" # Base path for all fetch data and workspace env setting BaseFetchPath = r"" lake = "" # Lake abbreviation lakefull = "" # Full name of lake buoy = "" # Buoy name method ="REI" # Method: 'EF' or 'REI' # List of years to process years = ['2006', '2007'] # File geodatabase (.gdb) containing the 36 directional rasters # for each lake from Rohweder et al. 2012 tool output inputDirRasFolder = r"" # Output file geodatabase path (.gdb) outputGDB = r"" # Final lake-wide composite geodatabse lakewide_compositeGDB = r"" ################################################################################### # Begin data processing # Loop through each year in the years list for year in years : print ('\n', lake, 'buoy is', buoy, 'year is', year, 'and method is', method, '\n') # Generate (temporary) fetch calc GDB fetchCalcGDB = "fetch_calc" + ".gdb" fetchCalcPath = os.path.join(BaseFetchPath, fetchCalcGDB) # If temporary fetch calc GDB exists, then do not create arcpy.env.workspace = BaseFetchPath if fetchCalcPath in arcpy.ListWorkspaces("", "FileGDB") : pass else: arcpy.CreateFileGDB_management (BaseFetchPath, fetchCalcGDB) # Read in the text file containing wind direction, directional frequency #(as %), and average velocity with open(wind_text) as F : reader = csv.reader(F, quoting=csv.QUOTE_NONNUMERIC) wndDir, dirFreq, wndSpd = zip(*reader) wndDir = list(wndDir) dirFreq = list(dirFreq) wndSpd = list(wndSpd) # Naturally sort directional rasters def sorted_nicely( l ): """ Sort the given iterable in the way that humans expect.""" convert = lambda text: int(text) if text.isdigit() else text alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ] return sorted(l, key = alphanum_key) # Read in each directional raster arcpy.env.workspace = inputDirRasFolder DirectionalRasterList = list() for x in sorted_nicely(arcpy.ListRasters("","")): DirectionalRasterList.append(x) # Check to see if final composites exist before running # Generate list of final composite rasters finalCompositeList = list() arcpy.env.workspace = lakewide_compositeGDB lakewideGDB = os.path.join(lakewide_compositeGDB, lakefull + '_lakewide_fetch_costAlloc_composite.gdb') lakewideGDBsList = arcpy.ListWorkspaces("", "FileGDB") finalCompositeName = lakefull + '_' + method + '_' + year + '_lakewide_costAlloc_composite' if lakewideGDB in lakewideGDBsList : arcpy.env.workspace = lakewideGDB finalCompositeList = arcpy.ListRasters("","") if finalCompositeName in finalCompositeList : print ('Final lakewide compiled raster,', finalCompositeName, 'exists. Will not perform weighting of dir rasters') CompositeExists = True else : CompositeExists = False else : CompositeExists = False if CompositeExists == False : print ('Final lakewide compiled raster does not exist. Proceeding with operations.') # Read in names of all final summed rasters in output GDB arcpy.env.workspace = outputGDB outRasSumList = arcpy.ListRasters("","") # Before weighting directional rasters, check to see if final summed fetch output exists. outRasSumName = lake + '_' + 'fetch' + '_' + method + '_' + buoy + '_' + year arcpy.env.workspace = fetchCalcPath if outRasSumName in outRasSumList : print ('Summed raster', outRasSumName, 'exists, will not perform map algebra.') else : print ('Existing summed rasters are:') for each in outRasSumList : print (each) # Read in the directional rasters and perform map algebra on each (weighting by wind data) i = 0 for raster in DirectionalRasterList : print (i+1, 'of 36') outputWeightRasFilename = raster + 'w_' + lake + '_' + buoy + '_' + method + '_' + year inRas = os.path.join(inputDirRasFolder, raster) outPath = os.path.join(fetchCalcPath, outputWeightRasFilename) # if weighted output already exists, do not run if outputWeightRasFilename in arcpy.ListRasters("","") : print ('weighted raster', outputWeightRasFilename, 'already exists') i = i + 1 continue if method == 'EF' : # Weighting by frequency of direction only print ('Method is EF. Direction is', wndDir[i], 'and %frequency is', dirFreq[i]) tproc = time.strftime('%X %x %Z') print ('start of map algebra operation: ' + tproc) print (outputWeightRasFilename + ' = ' + str(DirectionalRasterList[i]) + '*' + str(dirFreq[i]) + ' for direction ' + str(wndDir[i])) outRaster = Raster(inRas) * (dirFreq[i]/100) outRaster.save(outPath) tproc_end = time.strftime('%X %x %Z') print ('End of map algebra operation: ' + tproc_end) print ('\n') i = i + 1 if method == 'REI' : # For weighting by frequency of direction AND velocity : print ('Method is REI. Direction is', wndDir[i], '%frequency is', dirFreq[i], 'and avg. wind speed is', wndSpd[i]) tproc = time.strftime('%X %x %Z') print ('start of map algebra operation: ' + tproc) print (outputWeightRasFilename + ' = ' + str(DirectionalRasterList[i]) + '*' + str(wndSpd[i]) + '*' + str(dirFreq[i]) + ' for direction ' + str(wndDir[i])) outRaster = Raster(inRas) * (dirFreq[i]/100) * (wndSpd[i]) outRaster.save(outPath) tproc_end = time.strftime('%X %x %Z') print ('End of map algebra operation: ' + tproc_end) print ('\n') i = i + 1 # Set output extent to "MINOF", which is the smallest # extent that all input rasters overlap (necessary for LM and LH) arcpy.env.extent = "MINOF" # Assign output name and directory outRasSumName = lake + '_' + method + '_' + buoy + '_' + year outRasSumPath = os.path.join(outputGDB, outRasSumName) # Weighted raster List W_rasterList = arcpy.ListRasters("fet*", "") key = lake + '_' + buoy + '_' + method + '_' + year keyList = [Ras for Ras in W_rasterList if key == Ras.split('_', 2)[2]] # Check to see if output raster already exists, if so continue # else complete cell statitics tool run. if outRasSumName in outRasSumList : print ('file ' + outRasSumName + ' already exists, will not calculate sum.') else : print ("Weighted rasters to be summed are:") for item in keyList : print (item) tproc = time.strftime('%X %x %Z') print ('Calculating Sum: ' + outRasSumName) outRasSum = CellStatistics(keyList, "SUM", "NODATA") outRasSum.save(outRasSumPath) # Delete all the directional weighted rasters in the keyList for tempRas in keyList : arcpy.Delete_management(tempRas) print ('\n', lake, ': Buoy was', buoy, '; year was', year, '; method was', method, '\n') # Print system end time and date print ("End of processing.") tend = time.strftime('%X %x %Z') print (tend)