[Tutor] ArcGIS Create a python script to generate north-facing aspect raster from digital elevation model

[Michael Omohundro]

Does anyone know how to create a python script to generate an aspect raster from the input elevation in a digital elevation model?


I need to specify TWO variables as user parameters: input elevation and output north-facing aspect. 
I need to create an aspect raster from the input elevation from a digital elevation model. 
I need to find the north facing aspect is the trick, between 0 and 22.5 and between 337.5 – 360. These are the north facing elevation ranges.
I need to reclass the north facing aspect into a 1/Nodata raster. 
Any cell is north facing should have value 0, other cells should be NoData. 
I need to save the results as a permanent raster named as my second input parameter. 
I can't use ModelBuilder then convert it into Python script.

Here is what I have so far:


# Created on: March 20 2015
# Usage: lab03-5 Aspect from raster surface
# Requirements: ArcGIS Desktop and Spatial Analyst Extension
# ---------------------------------------------------------------------------

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True

# Set environment settings.  dem30 is the digital elevation model that holds the elevation data.
elevation = r"F:\NW_Missouri_State_University\_Python_Class\Week_10\lab10\dem30"

# North facing elevation 1, range 0 to 22.5
inRange_North1 = range (0, 22.5, 0.5)
#North facing elevation 2,range 337.5 - 360
inRange_North2 = range (337.5, 360, 0.5)

# Set local variables
inRaster = "elevation"

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute Aspect
outAspect = Aspect(inRaster)

# Save the output 
outAspect.save("F:\NW_Missouri_State_University\_Python_Class\Week_10\lab10\TestLab10_5")

# Specify the current workspace as the workspace where the input elevation raster is. 
# Extract the base name of the elevation raster.
arcpy.env.workspace = os.path.dirname(elevation)
inputElev = os.path.basename(elevation)

# Specify the output raster name as the input elevation name appending with “NorthFacing”.
saveReclass = arcpy.env.workspace + os.sep + inputElev + "NorthFacing"

# Check out the Spatial Analyst extension license for the script.
arcpy.CheckOutExtension("Spatial")

# Construct a map algebra expression to find the cell with elevation equal to the range values.
minRaster = arcpy.sa.Raster(inputElev) = inRange_North1
maxRaster = arcpy.sa.Raster(inputElev) = inRange_North2

# Construct a map algebra expression to perform a Boolean And 
# operation on the cell values of minRaster and maxRaster. 
# If the cell value is 1 in both raster, then set the output cell value as 1. 
# Otherwise, set the output cell value as 0. Save the output raster as variable outRaster.
outRaster = minRaster & maxRaster

# Create a remap object through RemapValue() function. 
# If the old value is 0, then set the new value as NODATA. 
# If the old value is 1, then set the new value as 1.

remap = arcpy.sa.RemapValue([[0, "NODATA"], [1, 1]])
outReclassify = arcpy.sa.Reclassify(
    outRaster, "Value", remap, "NODATA")

#Call the save method on the raster object to save the raster as permanent dataset.
outReclassify.save(saveReclass)

# Check in the Spatial Analyst extension license.
arcpy.CheckInExtension("Spatial")