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")
_______________________________________________
Tutor maillist - Tutor@python.org
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor
