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