#!/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 os, time, sys import arcpy from arcpy.sa import * # Print system start time and date tstart = time.strftime('%X %x %Z') print (tstart) # Check out the ArcGIS Spatial Analyst extension license arcpy.CheckOutExtension("Spatial") ########################SET CODE PARAMETERS################################# # Specify Lake lake = "" # Lake abbreviation lakefull = "" # Full name of lake method = 'REI' # Method: 'EF' or 'REI' # List of years to process years = ['2006', '2007'] # Raster depicting areas of land or water # also used to set the reference raster settings for # the Environments.... Raster Snap & Raster Cell Size & Raster Extent # water pixel value = 1; land pixel value = 2 refRas = r"" # Point feature class containg buoy locations buoypath = r"" # Name of field containing the buoy name (e.g., "45002") buoyName = "" # File geodatabase containing effective fetch or REI lakewide rasters lakeDataGDB = r"" # Output path outputGDB = r"" ############################################################ # Environment settings arcpy.env.workspace = lakeDataGDB arcpy.env.cellSize = refRas arcpy.env.snapRaster = refRas # Begin cost distance allocation processing # Run cost allocation tool costAllocate = lake + "_cost_allocate" ### Cost allocation tool ##arcpy.env.extent = refRas ##print ('Running cost allocation tool... may take 5-10 minutes...') ##print ('') ##arcpy.gp.CostAllocation_sa(buoypath, refRas, costAllocate, ## "1000000", "", "OBJECTID", "","") ##print ('...Cost Allocation complete') ##print ('') # Create feature layer from buoy points arcpy.MakeFeatureLayer_management (buoypath, "buoy_lyr", "") # Create buoy list and use Search Cursor to # populate with buoy names from buoy_lyr buoy_List = list() rows = arcpy.SearchCursor("buoy_lyr") fields = arcpy.ListFields("buoy_lyr", buoyName) for row in rows: for field in fields: if field.name == buoyName: buoy_List.append(row.getValue(field.name)) print ('Buoy list is: ') print (buoy_List) print ('') # Use the Search Cursor to move through the cost allocate # attribute table, and populate a list using the 'Value' column costVal_list = list() rows = arcpy.SearchCursor(costAllocate) fields = arcpy.ListFields(costAllocate,"Value") # Extract out each buoy region from the Cost Allocate raster # and use that as a mask to extract the appropriate region from # each single-buoy, lake-wide weighted fetch raster for row in rows : for field in fields : if field.name == "Value" : costVal = row.getValue(field.name) costVal_list.append(costVal) arcpy.env.workspace = lakeDataGDB where_clause = "Value=" + str(costVal) # Match the 'Value' with the buoy in the buoy list currentBuoy = buoy_List[costVal-1] print ('Current buoy is ' + currentBuoy) costExtName = "extracted_by_" + currentBuoy costExtractPath = os.path.join(lakeDataGDB, costExtName) print ('') # Creates an extraction mask for each buoy region from the cost Allocation surface if costExtName in arcpy.ListRasters("extracted*", "" ) : print (costExtName + ' costAllocate extracted raster already exists. Skipping...') else : print (where_clause + 'and corresponding buoy is: ' + currentBuoy + ' ...extracting region from cost allocation surface...') print ('') costExtract = ExtractByAttributes(costAllocate, where_clause) costExtract.save(costExtractPath) print ('') print ('') for year in years : year = str(year) print ('Lake is ' + lakefull + ' and wind/fetch method is ' + method + ' and year is '+ year + '.') # Check if final composite raster exists lakewide_compositeRas = lakefull + '_' + method + '_' + year arcpy.env.workspace = outputGDB outRasterList = arcpy.ListRasters("","") if lakewide_compositeRas in outRasterList : print ('Final lakewide composite for ' + lakefull + ' ' + year + ' ' + method + ' already exsits.') else : print ('Begin geoprocessing steps to generate final composite') # Next, loop through the fetch raster list (which specifically contains only the fetch rasters # for the specified method and year) matching and extracting the correct fetch layer # Generates list of the appropriate wind weighted fetch rasters # sorted by lake, method and year range arcpy.env.workspace = lakeDataGDB for Ras in arcpy.ListRasters("","") : if method in Ras : if year in Ras: ExtractOutName = method + "_" + year + "_" + Ras.split("_")[-2] + "_ex_Ras" ExtractOutPath = os.path.join(outputGDB, ExtractOutName) if ExtractOutName in outRasterList: print ("Raster " + ExtractOutName + " already exists, moving to next raster.") else : print ("Extracting " + Ras + " with mask from " + Ras.split("_")[-2]) costExtract = arcpy.Raster("extracted_by_" + Ras.split("_")[-2]) ExtractOut = ExtractByMask(Ras, costExtract) ExtractOut.save(ExtractOutPath) # Switch workspace arcpy.env.workspace = outputGDB arcpy.env.extent = refRas # Combine each extracted raster into final composite print ("Combining each raster into final composite: ") i = 0 for raster in arcpy.ListRasters("","") : if year not in raster : print ("Skipping " + raster) continue if method not in raster : print ("Skipping " + raster) continue if "*composite" in raster : print ("Skipping " + raster) if 'ex_Ras' in raster : if i == 0 : print ("Using " + raster + " as part of composite. First raster.") conRas = arcpy.Raster(raster) i = i + 1 print (i) else : print ("Using " + raster + " as part of composite.") conRas = Con(IsNull(arcpy.Raster(raster)), conRas, arcpy.Raster(raster)) i = i + 1 print (i) conRas.save(lakewide_compositeRas) # Delete input partial rasters arcpy.env.workspace = outputGDB for raster in arcpy.ListRasters("","") : if 'ex_Ras' in raster : print ('Deleting raster inputs of composite... ' + raster) arcpy.Delete_management(raster) arcpy.Delete_management("buoy_lyr") # Print system end time and date print ("End of processing.") tend = time.strftime('%X %x %Z') print (tend)