Exercise 7: Risk Model for Landslide Susceptibility in Oregon

Goals and Objectives

The purpose of this exercise was to generate a risk model for landslides in Oregon using Python scripting.  

Methods

In this exercise, the data from Exercise 5 and Exercise 6 was used to create a risk model for landslides in Oregon with the help of Python.  While developing the script, comments were utilized to annotate each section of code.  The following steps were taken:

• importing system modules

• importing Python modules: os, time, datetime, shutil, ArcGIS extensions

• creating smart variables and a list

• creating a fishnet

• buffering roads

• intersecting road buffer and fishnet

• creating variables for the reclass values of slope and landcover

• applying the reclass values using the reclassify tool

• multiplying the rasters together to create a risk raster

• using zonal statistics as a table to tally risk values

• joining results 

• debugging and running script

• creating a map in ArcMap

Results

The final script and map displaying my results are shown below.  

Figure 1: Map showing risk model results

#-------------------------------------------------------------------
# Name:        Ex7
# Purpose:    Generate a risk model for landslides in Oregon
#
# Author:      williaan
#
# Created:     25/07/2016
# Copyright:   (c) williaan 2016

#-------------------------------------------------------------------


#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\WILLIAAN\Ex_6\Ex6.gdb"
arcpy.env.overwriteOutput = True
#setting output CS to NAD 1983 HARN Oregon Statewide Lambert
arcpy.env.outputCoordinateSystem = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\WILLIAAN\Ex_5\Ex5\Ex5.gdb\OregonBoundary"
arcpy.CheckOutExtension('3D')
arcpy.CheckOutExtension('spatial')

#create variables
geodatabasepath = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\WILLIAAN\Ex_6\Ex6.gdb"
featureDatasetPath = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\WILLIAAN\Ex_6\Ex6.gdb\OregonFCs"

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

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)
roadBufferGridOutput = 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")

#buffering
arcpy.Buffer_analysis(majorRoads,roadBufferOutput,"100 Meters","FULL","ROUND","ALL","","PLANAR")

#intersecting roadway buffer and fishnet
print "Intersecting" ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Intersect_analysis([roadBufferOutput,outputFishnet],roadBufferGridOutput,"ALL","","INPUT")

#set reclass values for slope and landcover
#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; open water has a raster value of 11, and is 1 on the risk scale; developed open space has a raster value of 11 and risk value of 4
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)

#zonal stats as a table
arcpy.gp.ZonalStatisticsAsTable_sa(roadBufferGridOutput,bufferedRoadsUniqueIDField,landslideRiskRaster,landslideRiskStatistics,"DATA","MEDIAN")

#joining field
arcpy.JoinField_management(roadBufferGridOutput,bufferedRoadsUniqueIDField,landslideRiskStatistics,bufferedRoadsUniqueIDField,"Median")

print "Script is complete."