Exercise 5: Preparing Raster Data for Analysis Using Lists & Loops

Goals and Objectives

The purpose of this exercise was to gain experience with using standard raster geoprocessing and raster analysis tools.  This included project, clip, hillshade, slope, creating a list of rasters and a FOR IN loop.  

Methods

In this exercise, raster data about the state of Oregon was prepared for analysis with the help of Python.  While developing this script, comments were utilized to annotate each section of code.  The following steps were taken: 

• importing system modules
• importing Python modules: os, time, datetime, shutil, ArcGIS extensions
• creating variables
• creating a list of rasters
• setting up a FOR IN loop
• formatting raster names
• projecting and clipping a list of rasters
• creating a hillshade and calculating slope for a list of rasters
• merging the tiles together
• debugging and running script

Results 

The final script is displayed below.  After successfully running the script, the results were viewed in ArcMap. 

#-------------------------------------------------------------------
# Name:         Ex5
# Purpose:      Using standard raster geoprocessing tools in a FOR IN loop
#
# Author:      williaan
#
# Created:     18/07/2016

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

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

arcpy.env.workspace = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\WILLIAAN\Ex_5\Ex5\oregon"

#setting output coordinate system to NAD_1983_HARN_Oregon_Statewide_Lambert
arcpy.env.outputCoordinateSystem = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\WILLIAAN\Ex_5\Ex5\Ex5.gdb\OregonBoundary"

#check out Spatial Analyst and 3D Extensions
arcpy.CheckOutExtension('3D')
arcpy.CheckExtension('spatial')


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

workspace = r"Q:\StudentCoursework\CHupy\GEOG.491.801.2167\WILLIAAN\Ex_5\Ex5\oregon"
geodatabasePath = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\WILLIAAN\Ex_5\Ex5\Ex5.gdb"


#set boundary for clipping all rasters (Oregon State boundary)
clipBoundary = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\WILLIAAN\Ex_5\Ex5\Ex5.gdb\OregonBoundary"

#list for all clipped rasters (for future mosaicing)
clipList = []

#list for hillshades
hsList = []

#list for slopes
slopeList = []

#names for mosaic operations
mosaicedDEM = "NW_Oregon_Elevation_30m"
mosaicedHS = "NW_Oregon_Hillshade_30m"
mosaicedSlope = "NW_Oregon_Slope_30m"

print "Obtaining a list of rasters"
rasters = arcpy.ListRasters("","IMG")

for raster in rasters:
    print "Formatting the raster's name"
    #obtaining just basename of raster (no file extension or path)
    outName = os.path.basename(raster).rstrip(os.path.splitext(raster)[1])
    print "Unformatted name:" + outName
    outName = str(outName)
    #using string 'slice' to remove first three characters from filename
    fileName = str(outName[3:])
    print "Formatted name:" + fileName
    #set all new raster names
    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 rasterProj to make sure name formatting is working
    print rasterProj
    #project raster
    print "Projecting" + fileName ,datetime.datetime.now().strftime("%H:%M:%S")
    arcpy.ProjectRaster_management(raster,rasterProj,"","BILINEAR")
    #clip projected raster to state boundary
    print "Clipping" + fileName
   arcpy.Clip_management(rasterProj,"",rasterClip,clipBoundary,"","ClippingGeometry","NO_MAINTAIN_EXTENT")
    clipList.append(rasterClip)
    #run hillshade on clipped rasters
    print "Hillshading" + fileName
    arcpy.HillShade_3d(rasterClip,rasterHS)
    hsList.append(rasterHS)
    #calculate slope on list of clipped rasters
    print "Calculating slope of" + fileName
    arcpy.Slope_3d(rasterClip,rasterSlope)
    slopeList.append(rasterSlope)
    print "Loop is finished." ,datetime.datetime.now().strftime("%H:%M:%S")

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