Exercise 5: Preparing Raster Data For Analysis Using Raster Lists and Loops

The Objectives

The goal of this assignment was to set up a loop for formatting, projecting, clipping, hill shading, and calculating a list of raster files.

The Process

In order to complete this task, a series of commands were put into Python as script. First, the arcpy, date and time, spatial analysis extension, and 3-D analysis extension modules were imported to Python. Then, a work space was set up and overwrite was enabled to allow for data manipulation and saving: "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex5\oregon". After, the coordinate system was established and variables were created. Next, raster files from the Ex5 geodatabase were listed, pulling only .img files. A "FOR IN" loop was implemented to allow the following list of rasters to undergo all of the manipulations (ie. clip, project, hill shade, etc) established in the variable creation step. This has a similar function as the iterator tool in model builder used in exercise 2. Lastly, the files were merged by mosaicing the manipulated raster files.

The Results

The script generated four shape files and four raster files to be used together.

#-------------------------------------------------------------------------------
# Name:        Exercise 5
# Purpose:     Prepares Raster data for analysis
#
# Author:      Zach Miller
#
# Date:        7/20/2016
#-------------------------------------------------------------------------------

#ensure the program is working
print"script initialized"

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

# Allow the script to overwrite existing files
arcpy.env.workspace = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex5\oregon"
arcpy.env.overwriteOutput = True

#The next line needs to be set to your EX5 geodatabase's Oregon Boundary feature class
arcpy.env.outputCoordinateSystem = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex5\Ex5.gdb\OregonBoundary"

#"NAD_1983_HARN_Oregon_Statewide_Lambert" projection
arcpy.CheckOutExtension('3D')
arcpy.CheckOutExtension('spatial')
#-------------------------------------------------------------------------------

#create variables
print "Obtaining a List of rasters" ,datetime.datetime.now().strftime("%H:%M:%S")

#The workspace is the folder holding the rasters
workspace = r"Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex5\oregon"
geodatabasePath = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex5\Ex5.gdb"

# Set the boundary used to clip all of the rasters to the Oregon State boundary
clipBoundary = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex5\Ex5.gdb\OregonBoundary"

#This list will hold all of the clipped rasters for a future mosaicing operation
clipList = []

#This list is similar to the clip list, but is for hillshades
hsList = []

#This list is for slope rasters
slopeList = []

#The following names will be used for the three mosaic operations: Elevation, Hillshade, and Slope
mosaicedDEM = "NW_Oregon_Elevation_30m"
mosaicedHS = "NW_Oregon_Hillshade_30m"
mosaicedSlope = "NW_Oregon_Slope_30m"
#-------------------------------------------------------------------------------

#create a list of rasters with .img extension
print "Obtaining a list of rasters"
rasters = arcpy.ListRasters("","IMG")
#-------------------------------------------------------------------------------

#Set up "for in" loop
for raster in rasters:
    print "Formatting the raster's name"
#-------------------------------------------------------------------------------
#using os.path.basename... we obtain just the base name of the raster, without the file extension or source path.
outName = os.path.basename(raster).rstrip(os.path.splitext(raster)[1])
print "Unformatted name: " + outName
outName = str(outName)

#We use a string 'slice' to remove the first three characters from the filename...essentially removing the img, and leaving only the coordinates
fileName=str(outName[3:])
print "Formatted name: " + fileName

#Next, using the newly reformatted filename, we set all of the names for the new rasters we are to create.
rasterProj = os.path.join(geodatabasePath,fileName) + "_Proj"
rasterClip = os.path.join(geodatabasePath,fileName) + "_Clip"
rasterHS = os.path.join(geodatabasePath,fileName)+"_HS"
rasterSlope = os.path.join(geodatabasePath,fileName) + "_Slope"

#print raster projection to make sure the name formatting is working.
print rasterProj
#-------------------------------------------------------------------------------

#We already set the coordinate system using the 'output coordinate system' environment setting so we don't need to worry about it here.
# Be sure to use "BILINEAR" or "CUBIC" resampling when projecting elevation and other continuous datasets, as they will otherwise have incorrect values.
print "Projecting " + fileName , datetime.datetime.now().strftime("%H:%M:%S")
arcpy.ProjectRaster_management(raster,rasterProj,"","BILINEAR")
#-------------------------------------------------------------------------------

#The projected raster will be clipped to the state boundary, and "NO_MAINTAIN_EXTENT" should be set so the alignment of the pixels isn't changed by the clip function
print "Clipping " + fileName
arcpy.Clip_management(rasterProj,"",rasterClip,clipBoundary,"","ClippingGeometry","NO_MAINTAIN_EXTENT")
clipList.append(rasterClip)
#-------------------------------------------------------------------------------

#run a hillshade on the clipped rasters
print "Hillshading " + fileName
arcpy.HillShade_3d(rasterClip,rasterHS)
hsList.append(rasterHS)
#-------------------------------------------------------------------------------

#calculate slope on the clipped rasters
print "Calculating slope of " + fileName
arcpy.Slope_3d(rasterClip,rasterSlope)
slopeList.append(rasterSlope)
print "loop is complete"
#-------------------------------------------------------------------------------

#merging the tiles
print "Mosaicing the clipped DEMS"
arcpy.MosaicToNewRaster_management(clipList,geodatabasePath,mosaicedDEM,"","32_BIT_FLOAT","",1)
print "Mosaicing the hillshaded DEMS"
arcpy.MosaicToNewRaster_management(hsList,geodatabasePath,mosaicedHS,"","8_BIT_UNSIGNED","",1)
print "Mosaicing the slope rasters"
arcpy.MosaicToNewRaster_management(slopeList,geodatabasePath,mosaicedSlope,"","32_BIT_FLOAT","",1)

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

Sources

All sources provided by Dr. Christina Hupy