## MGLPhase2p.py

 1 ```# MGLPhase2.py ``` ```# Created on: Wed Oct 7 2009 09:35:35 AM ``` ```# by William Perry & Robert Lugo ``` ```# Usage: MGLPhase2 ``` ```# Description: Master python script for second phase of the MGL model ``` ```# This program does; ``` ```# 1. Creates a surface with genetic values and values of 0 for areas outside of species range but within total analysis extent for each species ``` ```# 2. Creates a surface with values of 1 for each cell with a genetic value, 0 for other cells in analysis extent ``` ```# 3. Sum of all genetic values created in step 1. ``` ```# 4. Sum of counts created in step 2. ``` ```# 5. Divide step 3 layer by step 4 layer to get the average of all species ``` ```# 6. Calculate the variance using this equation: s^2 = (sum of x^2 - ((sum of x)^2/n)/n-1, where ``` ```# s^2 = variance ``` ```# x = input genetic values (from the individual species genetic grids) ``` ```# n= number of input values (from the result of step 4) ``` ```# 7. Clean up files ``` ```# ``` ```# Parameter list ``` ```# Parameter Properties ``` ```# Display Name Data type Type Direction MultiValue ``` ```# argv Input folder Folder Required Input No ``` ```# argv Output folder Folder Required Input No ``` ```# argv Extent grid Grid Required Input No ``` ```# Import system modules ``` ```import sys, string, os, arcgisscripting ``` ```# Create the Geoprocessor object ``` ```gp = arcgisscripting.create(9.3) ``` ```# set variables ``` ```inFolder = sys.argv ``` ```outFolder = sys.argv ``` ```outExtent = sys.argv ``` ```# Check out any necessary licenses ``` ```gp.CheckOutExtension("spatial") ``` ```############################### future code to implement ``` ```#try: ``` ```## gp.CheckOutExtension("spatial") ``` ```# else: ``` ```# raise "LicenseError" ``` ``` #gp.Workspace = "D:/GrosMorne" ``` ``` #gp.HillShade_sa("WesternBrook", "westbrook_hill", 300) ``` ``` #gp.Aspect_sa("WesternBrook", "westbrook_aspect") ``` ``` #gp.CheckInExtension("spatial") ``` ```#except "LicenseError": ``` ```# print "Spatial Analyst license is unavailable" ``` ```#except: ``` ```# print gp.GetMessages(2) ``` ```############################### ``` ```# Load required toolboxes... ``` ```ArcGISHome = os.environ["ARCGISHOME"] ``` ```Toolpath = ArcGISHome + "ArcToolbox\\Toolboxes\\Spatial Analyst Tools.tbx" ``` ```print Toolpath ``` ```gp.AddToolbox(Toolpath) ``` ```# Set the Geoprocessing environment... ``` ```# set the analysis extent ``` ```gp.extent = outExtent ``` ```#gp.extent = 'c:/files/gis/mglmodel/tooldata/ecogridarea' ``` ```# select directory to look in ``` ```gp.Workspace = inFolder ``` ```Outpath = outFolder # select directory to write to ``` ```### Step 1. Reads in the genetic surfaces and changes the null values to 0 for areas ``` ```### outside of species range but within total analysis extent ``` ```gp.Addmessage("Changing Null values to zeros...") ``` ```# Get a list of grids in the workspace. ``` ```rasters = gp.ListRasters("", "GRID") ``` ```# create a variable to count the number of grids ``` ```x = 0 ``` ```for raster in rasters: ``` ``` # loop through rows and process data.... ``` ``` print "selected grid for calculations is " + raster ``` ``` ``` ``` # Local variables... ``` ``` #out_grid = fc[:13] #only use the first 13 characters ``` ``` temp05 = Outpath + "/" + raster + "_0" ``` ``` #temp05 = Outpath + "/" + raster + "_0" ``` ``` gp.SingleOutputMapAlgebra_sa("con(isnull(" + raster + " ), 0 , " + raster + ")", temp05, "") ``` ``` x = x + 1 ``` ``` ``` ```### Step 2. Creates a surface with values of 1 for each cell with a genetic value, ``` ```### 0 for other cells in analysis extent ``` ```gp.Addmessage("Changing Genetic values to ones...") ``` ```# Get a list of grids in the workspace. ``` ```rasters = gp.ListRasters("", "GRID") ``` ```for raster in rasters: ``` ``` # loop through rows and process data.... ``` ``` print "selected grid for calculations is - " + raster ``` ``` # Local variables... ``` ``` temp05 = Outpath + "/" + raster + "_1" ``` ``` ``` ``` #gp.SingleOutputMapAlgebra_sa("con(isnull(" + raster + " ), 0 , " + raster + ")", temp05, "") ``` ``` gp.SingleOutputMapAlgebra_sa("con(isnull(" + raster + " ), 0 , 1)", temp05, "") ``` ``` ``` ```### Step 3. Sum of all genetic values created in step 1 ``` ```# now need to switch the workspace to the outfolder ``` ```gp.Workspace = outFolder ``` ```# Get a list of grids in the workspace. ``` ```rasters = gp.ListRasters("*_0", "GRID") ``` ```# Local variables... ``` ```theString = "" ``` ```#out_grid = fc[:13] #only use the first 13 characters ``` ```gp.Addmessage("Adding all the genetic grids...") ``` ```for raster in rasters: ``` ``` # loop through rows and process data.... ``` ``` theString = theString + str(raster) + ", " ``` ``` ``` ```theString = theString.rstrip(", ") ``` ```print theString ``` ``` ``` ```try: ``` ``` Output = "allsumgv" ``` ``` InExpression = "Sum(" + theString + ")" ``` ``` gp.SingleOutputMapAlgebra_sa(InExpression, Output, theString) ``` ``` ``` ```except: ``` ``` print "exception:" ``` ``` print gp.GetMessages(2) ``` ```# now we need to square the sum of all grids from SumAllGeneticValues.py, AllSumGV. ``` ```try: ``` ``` Inraster = 'allsumgv' ``` ``` temp06 = "allsumgvsq" ``` ``` #temp06 = Outpath + "/" + "AllSumGVSq" ``` ``` # Process: Square ``` ``` gp.Square_sa(Inraster, temp06) ``` ```except: ``` ``` print "exception:" ``` ``` print gp.GetMessages(2) ``` ``` ``` ```### Step 4. Sum of counts created in step 2. ``` ```# Get a list of grids in the workspace. ``` ```rasters = gp.ListRasters("*_1", "GRID") ``` ```# Local variables... ``` ```theString = "" ``` ```#out_grid = fc[:13] #only use the first 13 characters ``` ```gp.Addmessage("Adding all the count grids...") ``` ```for raster in rasters: ``` ``` # loop through rows and process data.... ``` ``` print "selected grid for calculations is " + raster ``` ``` theString = theString + str(raster) + ", " ``` ``` ``` ```theString = theString.rstrip(", ") ``` ```print theString ``` ``` ``` ```try: ``` ``` Output = "allcountgv" ``` ``` InExpression = "Sum(" + theString + ")" ``` ``` gp.SingleOutputMapAlgebra_sa(InExpression, Output, theString) ``` ``` ``` ```except: ``` ``` print "exception:" ``` ``` print gp.GetMessages(2) ``` ``` ``` ```### Step 5. Divide step 3 layer by step 4 layer to get the average of all species ``` ```## ``` ```# The first step is to take all the single species genetic divergence surfaces and square them, then sum them. ``` ```# Get a list of grids in the workspace. The '_0' grids contain 0 instead of null values. ``` ```rasters = gp.ListRasters("*_0", "GRID") ``` ```theString = "" ``` ```gp.Addmessage("Squaring all the genetic values grids, then adding them...") ``` ```# loop through rows and process data.... ``` ```for raster in rasters: ``` ``` ``` ``` # Local variables... ``` ``` temp05 = raster[:-2] + "_2" # here we are stripping off the '_0' and replacing with '_2' to indicated a squared grid. ``` ``` # Process: Square ``` ``` gp.Square_sa(raster, temp05) ``` ``` ``` ```# Now get a list of squared grids in the workspace and add them together. ``` ```rasters2 = gp.ListRasters("*_2", "GRID") ``` ```gp.Addmessage("Adding all the squared genetic values grids...") ``` ```for raster in rasters2: ``` ``` # loop through rows and process data.... ``` ``` print "selected grid for calculations is " + raster ``` ``` theString = theString + str(raster) + ", " # builds expression of raster names separated by commas. ``` ``` ``` ```theString = theString.rstrip(", ") ``` ```print theString ``` ```try: ``` ``` Output = "allsqsumgv" ``` ``` InExpression = "Sum(" + theString + ")" ``` ``` gp.SingleOutputMapAlgebra_sa(InExpression, Output, theString) ``` ``` ``` ```except: ``` ``` print "exception:" ``` ``` print gp.GetMessages(2) ``` ``` ``` ```print ``` ```print "Finished" ``` ```print "Number of Grids Processed: " + str(x) ``` ```# Next we need to change all the zero values in allcountgv back to null and remove all the grid cells with count < 2. ``` ```try: ``` ``` # Set local variables ``` ``` InRaster = "allcountgv" ``` ``` OutRaster = "allcountgvn" ``` ``` # Process: SetNull ``` ``` gp.SetNull_sa(InRaster, InRaster, OutRaster, "Value < 2") ``` ```except: ``` ``` # If an error occurred while running a tool, then print the messages. ``` ``` print gp.GetMessages() ``` ```gp.Addmessage("Creating average genetic value grid...") ``` ``` ``` ```try: ``` ``` # Set the input raster dataset ``` ``` inputRaster1 = "allsumgv" ``` ``` inputRaster2 = "allcountgvn" ``` ``` # Set the output raster name ``` ``` outputRaster = "average" ``` ``` ``` ``` # Process: Divide ``` ``` gp.Divide_sa(inputRaster1, inputRaster2, outputRaster) ``` ```except: ``` ``` # If an error occurred while running a tool, then print the messages. ``` ``` gp.AddMessage("Failed to create raster: " + outputRaster) ``` ``` ``` ```### Step 6. Calculate the variance using this equation: s^2 = (sum of x^2 - ((sum of x)^2/n)/n-1, where ``` ```## s^2 = variance ``` ```## x = input genetic values (from the individual species genetic grids) ``` ```## n= number of input values (from the result of step 4) ``` ``` ``` ```# Now we need to build a single output Map Algebra expression for the variance equation below; ``` ```#Gvariance = ([allsqsumgv] - ([allsumgvsq] / [allcountgvn])) / ([allcountgvn] - 1) ``` ``` ``` ```gp.Addmessage("Creating genetic variance grid...") ``` ```try: ``` ```# Set the input raster dataset ``` ``` inputRaster1 = "AllSumGVsq" ``` ``` inputRaster2 = "AllCountGVn" ``` ``` # Set the output raster name ``` ``` outputRaster = "SumCntDivGV" ``` ``` # Process: Divide ``` ``` gp.Divide_sa(inputRaster1, inputRaster2, outputRaster) ``` ``` ``` ``` inputRaster3 = "Allsqsumgv" ``` ``` inputRaster4 = outputRaster ``` ``` ``` ``` #set the output raster ``` ``` outputRaster2 = "AminusB" ``` ``` ``` ``` # Process: Minus ``` ``` gp.Minus_sa(inputRaster3, inputRaster4, outputRaster2) ``` ``` inputRaster5 = "allcountgvn" ``` ``` ``` ``` outputRaster3 = "MinusCount" ``` ``` ``` ``` # Process: Minus ``` ``` gp.Minus_sa(inputRaster5, 1, outputRaster3) ``` ``` ``` ``` FinalOutput = "variance" ``` ``` ``` ``` # Process: Divide ``` ``` gp.Divide_sa(outputRaster2, outputRaster3, FinalOutput) ``` ``` ``` ``` ``` ```except: ``` ``` print "exception:" ``` ``` print gp.GetMessages(2) ``` ```### Step 7. Clean up and delete all the temporary grids from the Output folder ``` ```# first we need to copy the average and variance grids to another directory ``` ```gp.workspace = Outpath ``` ```gp.Addmessage("Cleaning up and removing temporary files....") ``` ```average_n = Outpath + "\\FinalResults\\average" ``` ```if gp.Exists(average_n): ``` ``` gp.delete_management(Outpath + "\\FinalResults\\average") ``` ```variance_n = Outpath + "\\FinalResults\\variance" ``` ```if gp.Exists(variance_n): ``` ``` gp.delete_management(Outpath + "\\FinalResults\\variance") ``` ```Final_n = Outpath + "\\FinalResults" ``` ```if gp.Exists(Final_n): ``` ``` gp.delete_management(Outpath + "\\FinalResults") ``` ```output_folder = "FinalResults" ``` ```gp.CreateFolder_management(Outpath, output_folder) ``` ```gp.copy_management("average", Outpath + "\\FinalResults\\average") ``` ```gp.copy_management("variance", Outpath + "\\FinalResults\\variance") ``` ```gp.copy_management("allcountgv", Outpath + "\\FinalResults\\sumcount") ``` ```rasters = gp.ListRasters("", "GRID") ``` ```for raster in rasters: ``` ``` # loop through and delete data.... ``` ``` print "deleting raster... " + raster ``` ``` ``` ``` gp.delete_management(raster) ``` ``` ``` ``` ```