Geog 491: Exercise 7- Risk Model for Landslide Susceptibility in Oregon

Goal and Background

The main purpose of this exercise is to generate a risk model for landslide in Northwestern Oregon.  The data used was prepared in previous exercises by projecting and clipping slope, land cover, precipitation, and hillshade.  The objectives are stated below:

• Set up a script, import modules, set environments and check out extensions.
• Set up smart variables and create a list.
• Create a fishnet to be used as the unit of analysis in combination with the roadway buffer. 
• Buffer the roads.
• Intersect the road buffer and the fishnet to create the unit of analysis feature class.
• Create a variable for the reclass values for slope and landcover. 
• Apply the reclass values using the reclassify function. 
• Multiply the rasters together to create the risk raster. 
• Use zonal statistics as a table to tally the risk values for the unit of analysis. 
• Join the results back to the unit of analysis. 
• Create a map of the results

Methods

Variables were created in order to give data to work with.  Field names were set as well in order to make things neat in the beginning.  A processing boundary was set in order to only get data in a certain area.  Geoprocessing tools buffer, intersect, and more were used.  

Results

The map and script for this exercise are pictured below:

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

# Name:        Exercise 7

# Purpose:     Generate a risk model for landslides in Oregon

#

# Author:      Laura Hartley, hartlelj

#

# Date:        7/26/2016

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

#Objective One

#Ensure the program is working

print "Script initialized"

#import system modules

import arcpy

from arcpy import env

import os

import shutil

import time

import datetime

from arcpy.sa import*

from arcpy.ddd import*

#Allow the script to overwrite existing files

arcpy.env.overwriteOutput = True

arcpy.env.workspace = "C:\Users\lhartley13\Documents\ArcGIS\EX6\Ex6.gdb"

arcpy.env.outputCoordinateSystem = "C:\Users\lhartley13\Documents\ArcGIS\EX5\Ex5\Ex5.gdb\OregonBoundary"

#NAD_1983_HARN_Oregon_Statewide_Lambert

arcpy.CheckOutExtension('3D')

arcpy.CheckOutExtension('spatial')

#Objective Two

#Create Variables

geodatabasePath = "C:\Users\lhartley13\Documents\ArcGIS\EX6\Ex6.gdb"

featureDatasetPath = "C:\Users\lhartley13\Documents\ArcGIS\EX6\Ex6.gdb\OregonFCs"

majorRoads = "C:\Users\lhartley13\Documents\ArcGIS\EX6\Ex6.gdb\OregonFCs\MajorRoads_NW"

landslidePts = "C:\Users\lhartley13\Documents\ArcGIS\EX6\Ex6.gdb\OregonFCs\LandslidePoints_NW"

precipData = "C:\Users\lhartley13\Documents\ArcGIS\EX6\Ex6.gdb\NW_Oregon_PrecipitationData_800m"

slopeData = "C:\Users\lhartley13\Documents\ArcGIS\EX5\Ex5\Ex5.gdb\NW_Oregon_Slope_30m"

landUseData = "C:\Users\lhartley13\Documents\ArcGIS\EX6\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"

#Set up field names

bufferedRoadsUniqueIDField = "FID_Oregon_NW_1kmFishnet"

#Objective Three

#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")

#Objective Four

#Process- Buffer

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

#Objective Five

#Process- Intersect

print "Intersecting!" ,datetime.datetime.now().strftime("%H:%M:%S")

arcpy.Intersect_analysis([roadBufferOutput,outputFishnet],roadBuffGridOutput,"ALL","","INPUT")

#Objective Six

#This string specifies different value ranges and the values these ranges will be assigned.

#The range 0-1.484953 will be assigned a value of 1.

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 '11' and is 1 on the risk scale. 'Developed, open space' has a raster values of '21' and a 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"

#Objective Seven

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

#Objective Eight

print "Multiplying the reclasses rasters to obtain a risk value"

arcpy.Times_3d(slopeReclassRaster,lulcReclassRaster,landslideRiskRaster)

#Objective Nine

#Process- Zonal statistics as table

arcpy.gp.ZonalStatisticsAsTable_sa(roadBuffGridOutput,bufferedRoadsUniqueIDField,landslideRiskRaster,landslideRiskStatistics,"DATA","MEDIAN")

#Objective Ten

#Process- Join Field

arcpy.JoinField_management(roadBuffGridOutput,bufferedRoadsUniqueIDField,landslideRiskStatistics,bufferedRoadsUniqueIDField,"MEDIAN")

print "script is complete" ,datetime.datetime.now().strftime("%H:%M:%S")

Sources

http://spatialdata.oregonexplorer.info/geoportal/catalog/main/home.page

http://nationalmap.gov/