Geog 491: Exercise 9- Create a Script Tool

Goal and Background

The goal of this exercise is to create a custom script tool.  The tool will generate a swiss hillshade.  A swiss hillshade is a slightly modified hillshade.  It uses several different filters to generate the terrain.  Below are the objectives for this exercise:

• Set the up the script.
• Set up the parameters using the “GetParatermersAsText” option.
• Set local variables.
• Set up a try except loop for the calculation of the hillshade.
• Add the workflow to the loop: divide the DEM, calculate a hillshade, run focal statistics, and add the DEM.
• Set up the except statement to print and error message.
• Create a custom toolbox.
• Use the Add Script Wizzard to add a custom script tool.
• Test your tool.

Methods

The regular parameters were set up to start the script.  Simple parameters were set up using GetParameterAsText.  Then the loops were set up.  A DEM from a previous exercise was created in order to produce the map in the resuts below.

Results

Below is the swiss hillshade created from the script and also the script that was created for this exercise.

#-------------------------------------------------------------------

# Name:        Exercise 9- Create a script tool

# Purpose:      To create a custom script tool.  The tools will generate a swiss hillshade

#

# Author:      Laura Hartley, hartlelj

#

# Date:         8/4/2016

#-------------------------------------------------------------------

#Ensure the program is working

print "Script initialized"

#import system modules

import os

import shutil

import time

import datetime

#Import arcpy module

import arcpy

from arcpy import env

from arcpy.sa import*

from arcpy.ddd import *

arcpy.env.overwriteOutput = True

arcpy.CheckOutExtension('3D')

arcpy.CheckOutExtension('spatial')

#Objective two

input_DEM = arcpy.GetParameterAsText(0)

scratchPath = arcpy.GetParameterAsText(1)

fileName = arcpy.GetParameterAsText(2)

filtered_Hillshade = arcpy.GetParameterAsText(3)

aerial_Perspective = arcpy.GetParameterAsText(4)

#Objective three

Z_factor = "1"

print "Creating local variables!" ,datetime.datetime.now().strftime("%H:%M:%S")

#Local variables:

divide_By = "5"

dEM_div_5 = os.path.join(scratchPath,fileName) +"_divHS"

default_Hillshade = os.path.join(scratchPath,fileName) +"_HS"

#Objective four

try:

    #Objective five, process: divide

    print "dividing input DEM by 5" ,datetime.datetime.now().strftime("%H:%M:%S")

    arcpy.Divide_3d(input_DEM,divide_By,dEM_div_5)

    #Process: Hillshade

    print "Create default hillshade" ,datetime.datetime.now().strftime("%H:%M:%S")

    arcpy.HillShade_3d(input_DEM,default_Hillshade,"315","45","NO_SHADOWS",Z_factor)

    #Process: Focal statistics

    print "Create a median filter for the hillshade effect" ,datetime.datetime.now().strftime("%H:%M:%S")

    outFocalStatistics = FocalStatistics(default_Hillshade,NbrRectangle(4,4,"CELL"),"MEDIAN","NODATA")

    outFocalStatistics.save(filtered_Hillshade)

    #Process: plus

    print "Add the divided dem to the default hilshade, creating an aerial perspective" ,datetime.datetime.now().strftime("%H:%M:%S")

    arcpy.Plus_3d(dEM_div_5,default_Hillshade,aerial_Perspective)

    print "The script is complete"

#Objective six

except:

    #report an error message

    arcpy.AddError("Something didn't work")

    #report any error messages it might have generated

    arcpy.AddMessage(arcpy.GetMessages(2))

Geog 491: Exercise 8- Raster Processing with Booleans and Loops

Goal and Background

The goal of this exercise was to use a suite of loops to clip and project a directory of scanned USGS topographic sheets into a seamless raster.  Below are the objectives for this exercise:

•  Set up the script by importing modules.
• Create variables, set up the workspace and create empty lists.
•  Create a spatial reference object using the coordinate system of a tile index.
• Use an IF statement to set up a loop to read files in a directory by employing a walk statement to traverse a workspace and find files with specific criteria.
• Measure the progress of the loop.
• Add all the files of the specified criteria to a raster list.
• For each of the rasters in the list – set up a processing loop. 
• Add a counter to track our progress.
• Forrmat the name of the raster.
• Project all the rasters.
• Create a polygon of the rasters topo footprint including the collar.
• Add the name of the raster to a new field.
• Add the tile to a list of tiles for a merge function.
• Merge the tiles.
• Create If statement to set up a for in loop.
• Create a search cursor and clip the rasters.
• Merge the rasters together.

Methods

Above is a flow chart of the workflow for this exercise.

Results

Below is the script created for this exercise.

#-------------------------------------------------------------------
# Name:        Exercise 8: Raster processing with Booleans and Loops
# Purpose:     To use a suite of loops to clip and project a directory of
#               scanned USGS topographic sheets into a seamless raster.
#
# Author:      Laura Hartley, hartlelj
#
# Date:         8/3/2016
#-------------------------------------------------------------------
#Ensure the program is working
print "script initialized"

#import system modules
import os
import shutil
import time
import datetime
import arcpy
from arcpy import env
from arcpy.sa import*
from arcpy.ddd import *

#Allow the script to overwrite existing files
arcpy.env.workspace = "C:\Users\lhartley13\Documents\ArcGIS\ex\EXAMPLE8\Topos"
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference("WGS 1984 UTM Zone 16N")
arcpy.CheckOutExtension('3D')
arcpy.CheckOutExtension('spatial')

#Create Variables
print "Obtaining a list of rasters" ,datetime.datetime.now().strftime("%H:%M:%S")

workspace = r"C:\Users\lhartley13\Documents\ArcGIS\ex\EXAMPLE8\Topos"

geodatabasePath = "C:\Users\lhartley13\Documents\ArcGIS\ex\EXAMPLE8\Ex8.gdb"

#Make a feature dataset in your geodatabase using the WGS84 UTM Zone #16N CS
featureDatasetPath = "C:\Users\lhartley13\Documents\ArcGIS\ex\EXAMPLE8\Ex8.gdb\WisconsinTiles"

#Make a scratch geodatabse for the clipped topos
clipGeodatabasePath = "C:\Users\lhartley13\Documents\ArcGIS\ex\EXAMPLE8\Ex8.gdb"

#Set up lists and file locations
rasters = []
tileList = []

mergedFolderLocation = "C:\Users\lhartley13\Documents\ArcGIS\ex\EXAMPLE8\Ex8.gdb"
mergedName = "KettleMoraine_Merged"

field = "FILENAME"

clipTileIndex = "C:\Users\lhartley13\Documents\ArcGIS\ex\EXAMPLE8\Ex8.gdb\WisconsinTiles\Wisconsin_24kTileIndex"

createIndex = True
clipRasters = True
mosaicRasters = True

mergedTiles = os.path.join(featureDatasetPath,mergedName)+"_TileFootprints"

#Create spatial reference object using the coordinate system of the #tile index
sr = arcpy.Describe(clipTileIndex).spatialReference

#If the boolean for create index is true, run the following code.
if(createIndex==True):
    walk = arcpy.da.Walk(workspace, type="TIF")
    i=0
    j=0
    k=0
    for dirpath, dirnames, filenames in walk:
        for filename in filenames:
            rasters.append(os.path.join(dirpath,filename))
            i=i+1
    for raster in rasters:
        k=i-j
        print str(k) + "rasters remaining"

        outName = os.path.basename(raster).rstrip(os.path.splitext (raster)[1])
        print outName
        outName = str(outName)
        outName2 = outName.replace(" ","_")

        print outName2

        #By splitting the string using the index of the 3rd character, we skip the WI_ in the filename
        shortened=str(outName[3:])
        endIndex=shortened.index("_")+3

        #Using these two new indicies, we make a new string of only the file's name
        name1 = str(outName[3:endIndex])
        name2 = name1.replace(" ","_")

        rasterProj = os.path.join(clipGeodatabasePath,name2)+ "_Proj"
        rasterTile = os.path.join(featureDatasetPath,name2)+ "_Footprint"

        #Objective ten
        print "Projecting" + name2 ,datetime.datetime.now().strftime("%H:%M:%S")
        arcpy.ProjectRaster_management(raster,rasterProj,"","BILINEAR")

        #Create a polygon of the topo's footprint, including the #collar
        print "Calculating raster domain"
        arcpy.RasterDomain_3d(rasterProj,rasterTile,"POLYGON")

        #Add field to hold the projected raster's path name in the #polygon's attribute table
        print "Adding a field"
        arcpy.AddField_management(rasterTile,field,"TEXT")

        #The following expression will be used to set the new #field's value to the projected
        #raster's name, using an update cursor
        expression = str(rasterProj)

        print "Calculating field" ,datetime.datetime.now().strftime("%H:%M:%S")
        with arcpy.da.UpdateCursor(rasterTile,field) as cursor:
            for row in cursor:
                row[0] = expression
                #update row must be used for the changes to take #effect
                cursor.updateRow(row)
        del cursor

        #Adding the tile to a list of tiles for a merge function
        tileList.append(rasterTile)

        j = j+1

        print ""

        #Merging the tiles to make a complete tile index
        print "Merging the tiles" ,datetime.datetime.now().strftime("%H:%M:%S")
        arcpy.Merge_management(tileList,mergedTiles)

#Objective fifteen
if(clipRasters==True):
    #creating a cursor to iterate through the rasters, tile by tile
    print "Clipping the collars off of the rasters" ,datetime.datetime.now().strftime("%H:%M:%S")
    cursorForTopos = arcpy.da.SearchCursor(mergedTiles, field)
    rastersToMosaic=[]

    for row in cursorForTopos:

        #Take the name from the filename field
        inputRaster = str(row[0])

        inputRasterName = os.path.basename(inputRaster).rstrip(os.path.splitext(inputRaster)[1])

        clipBoundary = "ClipBoundary"
        arcpy.MakeFeatureLayer_management(clipTileIndex,clipBoundary)

        rasterBoundary = "Raster Boundary"
        rasterNameField=arcpy.AddFieldDelimiters(mergedTiles,field)

        value = "'%s'" % inputRaster

        layerCreationSQLExpression = "%s = %s" % (rasterNameField,value)

        arcpy.MakeFeatureLayer_management(mergedTiles,rasterBoundary,layerCreationSQLExpression)

        arcpy.SelectLayerByLocation_management(clipBoundary,"HAVE_THEIR_CENTER_IN",rasterBoundary,"","NEW_SELECTION")

        outputRaster = os.path.join(clipGeodatabasePath,inputRasterName)+"_Clip"

        if arcpy.Exists(outputRaster):
            print inputRasterName + "has already been clipped"
            rastersToMosaic.append(outputRaster)
        else:
            print "Clipping" + inputRasterName ,datetime.datetime.now().strftime("%H:%M:%S")
            arcpy.Clip_management(inputRaster,"",outputRaster,clipBoundary,"","ClippingGeometry","NO_MAINTAIN_EXTENT")
            print inputRasterName + "Clipped" ,datetime.datetime.now().strftime("%H:%M:%S")
            rastersToMosaic.append(outputRaster)

        del cursorForTopos

        if(mosaicRasters==True):
            print "Mosaicing Topos" ,datetime.datetime.now().strftime("%H:%M:%S")
            arcpy.MosaicToNewRaster_management(rastersToMosaic,mergedFolderLocation,mergedName,sr,"8_BIT_UNSIGNED","",3)
            print "Mosaicing complete" ,datetime.datetime.now().strftime("%H:%M:%S")


print "script complete!" ,datetime.datetime.now().strftime("%H:%M:%S")