Monday, July 25, 2016

Exercise 6: Analyzing Raster Data for Landslide Susceptibility in Oregon

Goals and Objectives

The purpose of this exercise is to use a suite of raster and vector tools to identify common characteristics of landslides in Oregon.  Characteristics examined included precipitation, slope and land use/land cover.


Methods

In this exercise, raster data concerning the state of Oregon was used to examine landslide characteristics with the help of Python.  While developing this 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
  • adding raster values from precipitation, slope and land use to point feature class
  • adding and calculating a field to determine buffer radius
  • buffering landslides
  • using zonal statistics to calculate median slope of buffered landslides
  • joining results table to buffered landslides
  • creating and using a search and update cursor to update null values
  • calculating summary statistics for raster values of precipitation, slope and slide area based on movement class, land cover type and the combination of these two factors
  • using the tabulate area tool to calculate how much of each buffer falls within different analysis classes
  • exporting table to MS Excel
  • using a loop to delete intermediate feature classes
  • debugging and running script

Results

The final script is displayed below. 

 #-------------------------------------------------------------------
# Name:        Ex 6
# Purpose: To model landslide susceptibility 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"
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"
otherDataFolder = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\WILLIAAN\Ex_6"

 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"
#naming conventions
print "Setting up naming conventions" ,datetime.datetime.now().strftime("%H:%M:%S")
pointName = os.path.basename(landslidePts).rstrip(os.path.splitext(landslidePts)[1])
pointSelectOutput = os.path.join(geodatabasepath,pointName) + "_Selected"
pointsWithLandUse = os.path.join(geodatabasepath,pointName) + "_LU"
pointsWithSlope = os.path.join(geodatabasepath,pointName) + "_Slope"
pointsWithPrecip = os.path.join(geodatabasepath,pointName) + "_Precip"

 pointsBuffered = pointSelectOutput + "_Buffer"

 bufferedSlopeStatistics = pointsBuffered + "_Slope"
bufferedSlideTypeSummary = pointsBuffered + "_Summary_SlideType"
bufferedLandUseSummary = pointsBuffered + "_Summary_LULC"
bufferedBothSummary = pointsBuffered + "_Summary_Both"
bufferedTabulateLULC = pointsBuffered + "_TabulateLULC"

 tabulateLULC_ExcelTable = os.path.join(otherDataFolder,os.path.basename(bufferedTabulateLULC).rstrip(os.path.splitext(bufferedTabulateLULC)[1])) + ".xls"

 #creating list of feature classes to be created
print "Creating list of feature classes" ,datetime.datetime.now().strftime("%H:%M:%S")
createdFCs = []
createdFCs.append(pointSelectOutput)
createdFCs.append(pointsWithLandUse)
createdFCs.append(pointsWithSlope)
createdFCs.append(pointsWithPrecip)

 #field names
uniqueIDField = "UNIQUE_ID"
bufferDistanceField = "BufferDistance"
slideLengthmFieldName = 'Slide_Lenth_m'
rasterValueFieldName = "RASTERVALU"
lulcFieldName = "NLCD_2011"
moveClassFieldName = "MOVE_CLASS"
pointSlopeFieldName = 'Slope_Pt'
pointPrecipFieldName = 'PRECIP_mm'
bufferedSlopeFieldName = 'Slope_Buffered_Mean'
calculatedSlopeFieldName = 'Slope_Mean_Calculated'

 #processing of data will begin
#selecting landslides that have width and length greater than 0 and are in the debris flow, earth flow, or flow mvmt classes
print "Selecting landslides" ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Select_analysis(landslidePts,pointSelectOutput, "LENGTH_ft <>0 AND WIDTH_ft <>0 AND MOVE_CLASS IN ('Debris Flow', 'Earth Flow', 'Flow')")

 #add land use class info to the landslide points only for their specific location
ExtractValuesToPoints(pointSelectOutput,landUseData,pointsWithLandUse,"NONE","ALL")
arcpy.AlterField_management(pointsWithLandUse,rasterValueFieldName,"LULC_Value", "LULC Value")
#add slope info to the landslide points using interpolation
ExtractMultiValuesToPoints(pointsWithLandUse,[[slopeData,pointSlopeFieldName],[precipData,pointPrecipFieldName]],"BILINEAR")

 #add field
arcpy.AddField_management(pointsWithLandUse,bufferDistanceField, "FLOAT", "","","","","NULLABLE", "NON_REQUIRED","")
#calculate field
arcpy.CalculateField_management(pointsWithLandUse,bufferDistanceField, "(!LENGTH_ft!+!WIDTH_ft!)/2 *0.3048","PYTHON_9.3","")
#add field
arcpy.AddField_management(pointsWithLandUse,slideLengthmFieldName, "FLOAT", "","","","","NULLABLE", "NON_REQUIRED","")
#calculate field
arcpy.CalculateField_management(pointsWithLandUse,slideLengthmFieldName, "(!LENGTH_ft!+!WIDTH_ft!)*0.3048","PYTHON_9.3","")

 #buffer landslides
print "Buffering landslides" ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Buffer_analysis(pointsWithLandUse,pointsBuffered,bufferDistanceField,"FULL", "ROUND", "NONE", "", "PLANAR")

 #zonal stats as table
arcpy.gp.ZonalStatisticsAsTable_sa(pointsBuffered,uniqueIDField,slopeData,bufferedSlopeStatistics, "DATA", "ALL")

 #join field
arcpy.JoinField_management(pointsBuffered,uniqueIDField,bufferedSlopeStatistics,uniqueIDField,"Min;Max;Range;Mean;Std")

 #alter the joined field name for easier identification
arcpy.AlterField_management(pointsBuffered,"MEAN",bufferedSlopeFieldName,bufferedSlopeFieldName)

 #add field for the search and update cursor
arcpy.AddField_management(pointsBuffered,calculatedSlopeFieldName, "FLOAT", "","","","", "NULLABLE", "NON_REQUIRED","")
#create and use an update cursor to set the slope value for all that are not null, and the point value for features that have a null buffered value
print "Setting values based on buffered and point slope" ,datetime.datetime.now().strftime("%H:%M:%S")
listOfFields = [pointSlopeFieldName,bufferedSlopeFieldName,calculatedSlopeFieldName]

 with arcpy.da.UpdateCursor(pointsBuffered,listOfFields) as cursor:
    for row in cursor:
        if(row[1]==None):
            #if the buffered slope value [1] is null (None), set the value for the calculated slope [2] to the point value [0]
            print "Setting value to the point's slope value" ,datetime.datetime.now().strftime("%H:%M:%S")
            row[2]=row[0]
            cursor.updateRow(row)
        else:
                #if the buffered slope value [1] is valid, set the value for the calculated slope [2] to the buffered value [1]
                print "Setting value to the buffer's slope value" ,datetime.datetime.now().strftime("%H:%M:%S")
                row[2]=row[1]
                cursor.updateRow(row)
del cursor
#calculate summary stats for precip, slope, and slide area based on mvmt class
print "Calculating summary stats for precip, slope and slide area based on mvmt class" ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Statistics_analysis(pointsBuffered,bufferedSlideTypeSummary,[[pointPrecipFieldName,"MEAN"],[pointSlopeFieldName,"MEAN"],[bufferDistanceField,"MEAN"],[slideLengthmFieldName,"MEAN"]],moveClassFieldName)

 #calculate summary stats for precip, slope and slide area based on LULC class
print "Calculating summary stats for precip, slope and slide area based on LULC" ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Statistics_analysis(pointsBuffered,bufferedBothSummary,[[pointPrecipFieldName,"MEAN"],[pointSlopeFieldName,"MEAN"],[slideLengthmFieldName,"MEAN"]],[lulcFieldName,moveClassFieldName])

 #tabulate area
print "Tabulating the area of each LULC class for each buffered point" ,datetime.datetime.now().strftime("%H:%M:%S")
TabulateArea(pointsBuffered,uniqueIDField,landUseData,lulcFieldName,bufferedTabulateLULC,15)

 #export table to excel for additional analysis
arcpy.TableToExcel_conversion(bufferedTabulateLULC,tabulateLULC_ExcelTable)

 #delete intermediate feature classes
for fc in createdFCs:
    if arcpy.Exists(fc):
        print fc + "Exists, and will now be deleted"
        arcpy.Delete_management

 print "The script is complete."
Posted by at

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)