The goals of this exercise were to use the raster files created in exercise 5 in combination with precipitation and land use rasters provided in this exercise to calculate susceptibility of landslides in the study area. Ultimately, this produced an excel table with fields including: slope, precipitation, slide area, and land cover for landslides within the study area.
The Process
In order to produce this table raster files (and their data tables) were used from exercise 5 as well as new raster files "". First, the standard modules and extensions were imported and the work space was set up. Next, variables were created, this is where the rasters were brought into the script to be manipulated. Output files were then named and given a list of feature classes to belong to. Field names for the newly-created feature classes' attribute tables were then created. Then, the select tool was used to narrow down the type of landslide that is desired and those selected values were put into a point feature class. After, a buffer distance and a slide length field were created. The buffer was then executed and the resulting values added to the field. An UpdateCursor function was used to manage the values in the field. Next, summary stats were calculated for all of the landslide variables. Lastly, these values were tabulated to correlate land cover type and landslide area. The resulting table was then exported to Microsoft Excel for further analysis.
The Results
#-------------------------------------------------------------------------------
# Name: Ex6
# Purpose: Analyzes Raster Data
#
# Author: Zach Miller
#
# Date: 7/22/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"
otherDataFolder = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex6.gdb"
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"
landUseData = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex6.gdb\NLCD2011_Oregon_NW"
slopeData = "Q:\StudentCoursework\CHupy\GEOG.491.801.2167\MILLERZM\Ex5\Ex5.gdb\NW_Oregon_Slope_30m"
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"
buffered_TabulateLULC = pointsBuffered + "_TabulateLULC"
tabulateLULC_ExcelTable = os.path.join(otherDataFolder,os.path.basename(buffered_TabulateLULC).rstrip(os.path.splitext(buffered_TabulateLULC)[1]))+".xls"
createdFCs = []
createdFCs.append(pointSelectOutput)
createdFCs.append(pointsWithLandUse)
createdFCs.append(pointsWithSlope)
createdFCs.append(pointsWithPrecip)
# field names
uniqueIDField = "UNIQUE_ID"
bufferDistanceField = "BufferDistance"
slideLengthmFieldName = 'Slide_Length_m'
rasterValueFieldName = "RASTERVALU"
lulcFieldName = "NLCD_2011"
moveClassFieldName = "MOVE_CLASS"
pointSlopeFieldName = 'Slope_Pt'
pointPrecipFieldName = 'PRECIP_mm'
bufferedSlopeFieldName = 'Slope_Buffered_Mean'
calculatedSlopeFieldName = 'Slope_Mean_Calculated'
#-------------------------------------------------------------------------------
# Process: Select
print "Selecting!" ,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 information to the landslide points only for their specific locations
ExtractValuesToPoints(pointSelectOutput,landUseData,pointsWithLandUse,"NONE","ALL")
arcpy.AlterField_management(pointsWithLandUse,rasterValueFieldName,"LULC_Value","LULC Value")
# Add slope information to the landslide points using interpolation
ExtractMultiValuesToPoints(pointsWithLandUse,[[slopeData,pointSlopeFieldName],[precipData, pointPrecipFieldName]],"BILINEAR")
#-------------------------------------------------------------------------------
# Process: Add Buffer Field
arcpy.AddField_management(pointsWithLandUse, bufferDistanceField, "FLOAT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
# Process: Calculate Field
arcpy.CalculateField_management(pointsWithLandUse, bufferDistanceField, "(!LENGTH_ft!+ !WIDTH_ft!)/2 *0.3048", "PYTHON_9.3", "")
# Process: Add Slide Length Field
arcpy.AddField_management(pointsWithLandUse, slideLengthmFieldName, "FLOAT", "","", "", "", "NULLABLE", "NON_REQUIRED", "")
# Process: Calculate Field
arcpy.CalculateField_management(pointsWithLandUse,slideLengthmFieldName,"!LENGTH_ft!*0.3048", "PYTHON_9.3", "")
# Process: Buffer (2)
print "Buffering!", datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Buffer_analysis(pointsWithLandUse, pointsBuffered, bufferDistanceField, "FULL","ROUND", "NONE", "", "PLANAR")
#-------------------------------------------------------------------------------
# Process: Zonal Statistics as Table
arcpy.gp.ZonalStatisticsAsTable_sa(pointsBuffered, uniqueIDField, slopeData, bufferedSlopeStatistics, "DATA", "ALL")
# Process: Join Field
arcpy.JoinField_management(pointsBuffered, uniqueIDField, bufferedSlopeStatistics, uniqueIDField, "Min;Max;Range;Mean;Std")
#alter the joined field name so we can easily identify it later
arcpy.AlterField_management(pointsBuffered,"MEAN",bufferedSlopeFieldName,bufferedSlopeFieldName)
# Process: Add Field for the search and update cursor
arcpy.AddField_management(pointsBuffered, calculatedSlopeFieldName, "FLOAT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
# Process: Add Field for slide length
arcpy.AddField_management(pointsBuffered, slideLengthmFieldName, "FLOAT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")
# calculate field for slide length
arcpy.CalculateField_management(pointsBuffered, slideLengthmFieldName, "!Length_ft!*0.3048", "PYTHON_9.3", "")
#-------------------------------------------------------------------------------
# Create and use an update cursor to set the slope value to the buffered 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]
# The update cursor takes a list of fields and allows for the values to be changed by specific functions
with arcpy.da.UpdateCursor(pointsBuffered,listOfFields) as cursor:
# it allows for iteration by row (so each value will be calculated independently
for row in cursor:
# if the buffered slope value [1] is null (None), set the value for the calculated slope [2] to the point value [0]
# (FYI: lists of values use indicies of 0,1,2,etc)
if(row[1]== None):
#print "Setting value to the point's slope value"
row[2] = row[0]
cursor.updateRow(row)
# if the buffered slope value [1] is valid, set the value for the calculated slope [2] to the buffered value [1]
else:
#print "Setting value to the buffer's slope value"
row[2] = row[1]
cursor.updateRow(row)
del cursor
#-------------------------------------------------------------------------------
# Calculate summary statistics for precip, slope, and slide area, based on movement class
print "Calculating summary statistics for precip, slope, and slide area, based on movement class" ,datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Statistics_analysis(pointsBuffered,bufferedSlideTypeSummary,[[pointPrecipFieldName,"MEAN"],[pointSlopeFieldName,"MEAN"],[bufferDistanceField,"MEAN"],[slideLengthmFieldName,"MEAN"]],moveClassFieldName)
# Calculate summary statistics for precip, slope, and slide area, based on LULC class
print "Calculating summary statistics for precip, slope, and slide area, based on LULC",datetime.datetime.now().strftime("%H:%M:%S")
arcpy.Statistics_analysis(pointsBuffered,bufferedLandUseSummary,[[pointPrecipFieldName,"MEAN"],[pointSlopeFieldName,"MEAN"],[bufferDistanceField,"MEAN"],[slideLengthmFieldName,"MEAN"]],lulcFieldName)
# Calculate summary statistics for precip, slope, and slide area, based on movement class and LULC class
print "Calculating summary statistics 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"],[bufferDistanceField,"MEAN"],[slideLengthmFieldName,"MEAN"]],[lulcFieldName,moveClassFieldName])
#-------------------------------------------------------------------------------
#Tabulate area is used to identify how much of each buffer falls within different LULC classes
print "Tabulating the area of each LULC class for each buffered point"
TabulateArea(pointsBuffered,uniqueIDField,landUseData,lulcFieldName,buffered_TabulateLULC,15)
#Export the table to excel for additional analysis!
arcpy.TableToExcel_conversion(buffered_TabulateLULC,tabulateLULC_ExcelTable)
for fc in createdFCs:
if arcpy.Exists(fc):
print fc + "Exists, and will now be deleted"
arcpy.Delete_management(fc)
print "Script is complete"
Sources
All sources provided by Dr. Christina Hupy
No comments:
Post a Comment