QGIS » QGIS Application

# MGLPhase2.py

# Created on: Wed Oct 7 2009 09:35:35 AM

#   by William Perry & Robert Lugo

# Usage: MGLPhase2  <Input_folder> <Output_folder> <Extent>

# 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[1]  Input folder         Folder           Required  Input      No

#  argv[2]  Output folder        Folder           Required  Input      No

#  argv[3]  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[1]

outFolder = sys.argv[2]

outExtent = sys.argv[3]

# 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)