Exercise 7: Risk Model for Landslide Susceptibility in Oregon

The Objectives

The goal of this assignment was to take raster files from previous exercises as well as a Major Roads and Landslide points shape files to create a landslide susceptibility map depicting risk levels of major roads. This exercise used fishnet and zonal statistics functions to obtain the resulting map.

The Process

In order to accomplish this task, rasters from exercises 5 and 6 were used as well as roads and landslide points shape files. First, the standard modules and extensions were imported and variables were created. A square Kilometer fishnet was then created to be combined with rasters and other data later on to form a grid pattern for landslide susceptibility of roads. Next, a road buffer was created and intersected with the fishnet. Then, values were reclassified into five breaks, each with a different risk factor (slope and landcover rasters were used for this). Lastly, zonal statistics were used to find median values for the unit of analysis and the resulting table was joined with the analysis layer.

The Results

#-------------------------------------------------------------------------------
# Name:        Ex7
# Purpose:     Generates Risk Model for Landslides in Oregon
#
# Author:      Zach Miller
#
# Date:        7/26/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\Ex6.gdb"
arcpy.env.overwriteOutput = True
arcpy.env.outputCoordinateSystem = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex5\Ex5.gdb\OregonBoundary"

#NAD_1983_HARN_Oregon_Statewide_Lambert
arcpy.CheckOutExtension('3D')
arcpy.CheckOutExtension('spatial')
#-------------------------------------------------------------------------------

#create variables
geodatabasePath = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex6.gdb"
featureDatasetPath = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex6.gdb\OregonFCs"

majorRoads = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex6.gdb\OregonFCs\MajorRoads_NW"
landslidePts = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex6.gdb\OregonFCs\LandslidePoints_NW"
precipData = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex6.gdb\NW_Oregon_PrecipitationData_800m"
slopeData = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex5\Ex5.gdb\NW_Oregon_Slope_30m"
landUseData = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex6.gdb\NLCD2011_Oregon_NW"

boundaryName = "ProcessingBoundary"
fishnetName = "Oregon_NW_1kmFishnet"
roadsBufferName = "MajorRoads_NW_Buffered"
roadsBuffIntName = "MajorRoads_NW_Buffer_1kmGrid"
riskRasterName = "NW_Oregon_LandslideRisk"
pointName = os.path.basename(landslidePts).rstrip(os.path.splitext(landslidePts)[1])

processingBoundary = os.path.join(featureDatasetPath,boundaryName)
outputFishnet = os.path.join(featureDatasetPath,fishnetName)
roadBufferOutput = os.path.join(featureDatasetPath,roadsBufferName)
roadBuffGridOutput = os.path.join(featureDatasetPath,roadsBuffIntName)
slopeReclassRaster = slopeData + "_Reclass"
lulcReclassRaster = landUseData + "_Reclass"

landslideRiskRaster = os.path.join(geodatabasePath,riskRasterName)

landslideRiskStatistics = landslideRiskRaster + "_Stats"

#field names
bufferedRoadsUniqueIDField = "FID_Oregon_NW_1kmFishnet"
#-------------------------------------------------------------------------------

#Create processing boundary by extracting the footprint from the slope raster
print "Creating Processing boundary!" ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.RasterDomain_3d(slopeData,processingBoundary,"POLYGON")

#Create fishnet (Grid) with 1km spacing
#variables for fishnet
desc = arcpy.Describe(processingBoundary)
gridSize = 1000.0
print "Creating Fishnet!" ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.CreateFishnet_management(outputFishnet,str(desc.extent.lowerLeft),str(desc.extent.XMin) + " " + str(desc.extent.YMax+10),gridSize,gridSize,"0","0",str(desc.extent.upperRight),"NO_LABELS",processingBoundary,"POLYGON")
#-------------------------------------------------------------------------------

# Process: Buffer
arcpy.Buffer_analysis(majorRoads, roadBufferOutput, "100 Meters", "FULL", "ROUND", "ALL", "", "PLANAR")
#-------------------------------------------------------------------------------

# Process: Intersect
print "Intersecting!" ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Intersect_analysis([roadBufferOutput,outputFishnet], roadBuffGridOutput, "ALL", "", "INPUT")
#-------------------------------------------------------------------------------

#This string specifies different value ranges and the values these ranges will be assigned.
remapStringSlope = "0 1.484953 1;1.484954 6.184316 2;6.184317 10.883696 3;10.883697 15.583077 4;15.583078 29.681219 5;29.681220 34.380600 4; 34.380601 39.079981 3;39.079982 43.779362 2;43.779363 90 1"

#This string specifies risk values to the different nlcd classes based on their values.
remapStringNLCD = "11 1;21 4;22 4;23 3;24 1;31 2;41 4;42 5;43 5;52 5;71 4;81 3;82 3;90 3;95 1"
#-------------------------------------------------------------------------------

print "Reclassifying slope into individual run type classes " ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Reclassify_3d(slopeData,"VALUE",remapStringSlope,slopeReclassRaster,"NODATA")

print "Reclassifying landuse into individual run type classes " ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Reclassify_3d(landUseData,"VALUE",remapStringNLCD,lulcReclassRaster,"NODATA")
#-------------------------------------------------------------------------------

print "Multiplying the reclassed rasters to obtain a risk value"
arcpy.Times_3d(slopeReclassRaster,lulcReclassRaster,landslideRiskRaster)
#-------------------------------------------------------------------------------

# Process: Zonal Statistics as Table
arcpy.gp.ZonalStatisticsAsTable_sa(roadBuffGridOutput, bufferedRoadsUniqueIDField, landslideRiskRaster, landslideRiskStatistics, "DATA", "MEDIAN")
#-------------------------------------------------------------------------------

# Process: Join Field
arcpy.JoinField_management(roadBuffGridOutput, bufferedRoadsUniqueIDField, landslideRiskStatistics, bufferedRoadsUniqueIDField, "Median")

print "Script is complete"

Sources

All sources provided by Dr. Christina Hupy

Oregon Spatial Data Library
USGS National Map