Geog 491: Exercise 5- Preparing Raster Data for Analysis using Raster Lists and Loops

Goals and Background

This exercise focused on using standard raster geoprocessing tools and basic raster analysis tools in a FOR IN loop by creating a list of rasters.  The objectives are stated below,

• Set up a script, import modules, set environments and check out extensions.
• Set up smart variables.
• Create a list of rasters.
• Set up a FOR IN loop.
• Format the raster names in the list of rasters.
• Project a list of raster
• Clip a list of raster.
• Create a hillshade for a list of rasters.
• Calculate slope for a list of raster.
• Merge the resuts together.

Methods

Lists were created in order to list the final rasters when the script was done.  Then the loops were set up in order to run the script multiple times on all of the rasters.  The script can be seen below.

Results

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

# Name:        Exercise 5

# Purpose:     To use standard raster geoprocessing preparation     #              tools, such as

#              project, clip, as well as basic raster analysis tools #              such as,

#              hillshade and slope, in a FOR IN loop by creating a   #              list of rasters.

#

# Author:      Laura Hartley, hartlelj

#

# Date:        7/20/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\EX5\Ex5\oregon"

arcpy.env.overwriteOutput = True

#The next line needs to be set to your EX5 geodatabase's Oregon Boundary FC

arcpy.env.outputCoordinateSystem = "C:\Users\lhartley13\Documents\ArcGIS\EX5\Ex5\Ex5.gdb\OregonBoundary"

arcpy.CheckOutExtension('3D')

arcpy.CheckOutExtension('spatial')

#create variables

print "Creating Variables" ,datetime.datetime.now().strftime("%H:%M:%S")

#The workspace is the folder holding the rasters

workspace = r"C:\Users\lhartley13\Documents\ArcGIS\EX5\Ex5\oregon"

geodatabasePath = "C:\Users\lhartley13\Documents\ArcGIS\EX5\Ex5\Ex5.gdb"

#Set the boundary used to clip all of the raster to the Oregon State boundary

clipBoundary = "C:\Users\lhartley13\Documents\ArcGIS\EX5\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"

print "Obtaining a list of rasters"

rasters = arcpy.ListRasters("", "IMG")

#This command selects just a single raster from the list named 'rasters'

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

    #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 rasterProj, to make sure the name formatting is working

    print rasterProj

    #we already set the CS using the output CS env setting so we         #don't need to worry about it here.

    #Be sure to use 'BILLINEAR' 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_NAINTAIN_EXTENT" should be set so

    #the alignment of the pixels isn't changed

    print "Clipping" +fileName 
    ,datetime.datetime.now().strftime("%H:%M:%S")

    arcpy.Clip_management(rasterProj,"",rasterClip,clipBoundary,
    "","ClippingGeometry","NO_MAINTAIN_EXTENT")

    clipList.append(rasterClip)

    #creating a hillshade from the clipped raster

    print "Hillshading" +fileName
    ,datetime.datetime.now().strftime("%H:%M:%S")

    arcpy.HillShade_3d(rasterClip,rasterHS)

    hsList.append(rasterHS)

    #create a slope raster from the clipped raster

    print "Calculating slope of" + fileName

    arcpy.Slope_3d(rasterClip, rasterSlope)

    slopeList.append(rasterSlope)

    print "the loop has completed once"

#merging the tiles

print "Mosacing the clipped DEMs"

arcpy.MosaicToNewRaster_management(clipList,geodatabasePath,mosaicedDEM,"","32_BIT_FLOAT","",1)

print "Mosacing the hillshaded DEMs"

arcpy.MosaicToNewRaster_management(hsList,geodatabasePath,mosaicedHS,"","8_BIT_UNSIGNED","",1)

print "Mosaicing the slope raster"

arcpy.MosaicToNewRaster_management(slopeList,geodatabasePath,mosaicedSlope,"","32_BIT_FLOAT","",1)

print "Script Complete" ,datetime.datetime.now().strftime("%H:%M:%S")