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

[P McCombs]

On Mar 21, 2015 1:22 AM, "Michael Omohundro"
<momohund1972 at yahoo.com.dmarc.invalid> wrote:
>
> Does anyone know how to create a python script to generate an aspect
raster from the input elevation in a digital elevation model?
>
This topic is out of the normal scope of the mailing list. It's hard to say
if you are having a python issue or an ArcGIS issue.

I could possibly help with the latter,  but please provide the trace back
you get when it errors out,  or how the output differs from what you
expect.

Without the data you are using I'm not going to attempt to run your code.

>
> 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 *
In my experience this makes debugging harder, I can't tell if you are using
local functions or module functions without a detailed knowledge of the
module you imported.

> 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")
>
I'm guessing that you've resolved this problem in the past week,  but for
future problems I encourage you to check what data types are expected for
the functions you are calling. They can be unexpected sometimes in arcpy.

Paul McCombs