The goal of this exercise was to clip scanned raster files and merge them together using Booleans and loops.
The Process
- Import modules.
- Create variables, set up the work space and create empty lists.
- Create a spatial reference object (one tile index's coordinate system was used).
- Set up a loop.
- Put files of the specified criteria to a raster list.
- Set up processing loops for each raster.
- Format the raster name.
- Project all of the rasters.
- Create a polygon of the raster's topographic footprint.
- Add the name of the raster to a new field.
- Create a list of tiles for a merge function.
- Merge the tiles.
- Create "IF" statement to set up a "for, in" loop.
- Clip the rasters.
- Merge the rasters together.
The Results
#-------------------------------------------------------------------------------
# Name: Ex8
# Purpose: Processes rasters with booleans and loops
#
# Author: Zach Miller
#
# Created: 7/28/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 = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Exercise_8\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"Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Exercise_8\Topos"
geodatabasePath = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Exercise_8\Ex8.gdb"
#make a feature dataset in your geodatabase using the WGS84 UTM Zone 16N coordinate system
featureDatasetPath = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Exercise_8\Ex8.gdb\WisconsinTiles"
#Make a scratch geodatabase for the clipped topos
clipGeodatabasePath = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Exercise_8\Ex8.gdb"
rasters = []
tileList = []
mergedFolderLocation = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Exercise_8\Ex8.gdb"
mergedName = "KettleMoraine_Merged"
field = "FILENAME"
clipTileIndex = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Exercise_8\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):
#The walk function traverses a workspace, and finds files with the specified filetype
#in this instance, it finds TIF files. It contains the directory path, directory name, and file names within the 'walk' object
walk = arcpy.da.Walk(workspace, type="TIF")
#-------------------------------------------------------------------------------
#i, j, and k are variables that will be used used to know the progress of the loop
i = 0
j = 0
k = 0
#-------------------------------------------------------------------------------
for dirpath, dirnames, filenames in walk:
for filename in filenames:
#The script will iterate through all of the subfolders, and add all valid tifs to the list 'rasters'
rasters.append(os.path.join(dirpath, filename))
#For each valid raster, 1 value will be added to i...which will ultimately tell how many valid files there are
i = i+1
#-------------------------------------------------------------------------------
#The following operations will be performed on all of the rasters in 'rasters'
for raster in rasters:
#-------------------------------------------------------------------------------
#for the first iteration, j = 0 value 'k' will be # of valid rasters- 0.
#for subsequent iterations, it will count down until no rasters remain.
k = i-j
print str(k) + "rasters remaining"
#-------------------------------------------------------------------------------
#Formatting the name of the topo, in an initial attempt to replace the spaces with underscores
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:])
#Then we use the string 'index' function to identify the index of the next underscore in the string (it directly follows the end of the filename!)
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])
#use the 'replace' string function to completely remove the spaces from the resulting filenames
name2 = name1.replace(" ","_")
rasterProj = os.path.join(clipGeodatabasePath,name2) + "_Proj"
rasterTile = os.path.join(featureDatasetPath,name2)+"_Footprint"
#-------------------------------------------------------------------------------
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)
#-------------------------------------------------------------------------------
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])
#obtain only the base filename from the previous string
inputRasterName = os.path.basename(inputRaster).rstrip(os.path.splitext(inputRaster)[1])
clipBoundary="ClipBoundary"
arcpy.MakeFeatureLayer_management(clipTileIndex,clipBoundary)
rasterBoundary="RasterBoundary"
rasterNameField=arcpy.AddFieldDelimiters(mergedTiles,field)
#Within python, when one wants to put a string within certain quotes, they can specify this by using '%s' as a placeholder, then using a single % before the string they wish to substitute
#This will allow for the input raster to be recorded as 'inputRaster', within the string.
value = "'%s'" % inputRaster
#The next line allows for the creation of a sql query by using string substitution to replace both of the placeholders with the
#string name field and the newly formatted 'inputRaster' value
layerCreationSQLExpression = "%s = %s" % (rasterNameField, value)
#This line uses the newly created sql statement to select only the tile that matches the cursor value, which is incredibly
#important for future steps.
arcpy.MakeFeatureLayer_management(mergedTiles,rasterBoundary,layerCreationSQLExpression)
#This line selects the feature from the state tile index that has its center in the previously selected polygon
arcpy.SelectLayerByLocation_management(clipBoundary,"HAVE_THEIR_CENTER_IN",rasterBoundary,"","NEW_SELECTION")
outputRaster = os.path.join(clipGeodatabasePath,inputRasterName)+"_Clip"
# This logic exists to identify if the raster has already been generated for a specific tile, and skips the operation if it has already been performed
if arcpy.Exists(outputRaster):
print inputRasterName + " has already been clipped"
rastersToMosaic.append(outputRaster)
else:
# The raster will be clipped by the boundary identified using the select layer by location function.
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)
#This del function is essential, otherwise a schema lock may remain on your geodatabase, preventing you from accessing its contents until ArcMap is restarted.
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")
Sources
All sources provided by Dr. Christina Hupy
USGS National Map - http://nationalmap.gov/
No comments:
Post a Comment