Included here is the code for the python script used in ArcGIS. As previously mentioned this is a site specific code and requires the presense of certain GIS grid files to be present within the working geodatabase to execute properly. But in the spirit of open source, here it is anyway. The author of the code in no way suggests that this is an accurate representation of the true avalanche hazard, it is simply a tool that can be used in helping to assess the current conditions.
”’
Created on Nov 16, 2012
@author: wesleyashwood
”’
import urllib2
import arcpy
from arcpy.sa import *
#function to extract string from a given line of text
def extract(raw_string, start_marker, end_marker):
start = raw_string.index(start_marker) + len(start_marker)
end = raw_string.index(end_marker, start)
return raw_string[start:end]
#get text from url
url = urllib2.urlopen(“http://images.drivebc.ca/bchighwaycam/pub/weather/26224.html#”)
text = url.readlines()
url.close()
#specific lines to search
lines = [3,9,11,13,15,17,19]
#what to extract from each line
for i in lines:
if i == 3:
line = text[i]
start_marker = ‘right;background-color:#B9E5FB;padding:2px;”>’
end_marker = ‘</td></tr>’
time = extract(line,start_marker,end_marker)
if i == 9:
line = text[i]
start_marker = ‘0px;”>’
end_marker = ‘<sup>o</sup>’
temp = float(extract(line,start_marker,end_marker))
if i == 11:
line = text[i]
start_marker = ‘0px;”>’
end_marker = ‘mm (<a href=”#”‘
precip = int(extract(line,start_marker,end_marker))
if i == 13:
line = text[i]
start_marker = ‘0px;”>’
end_marker = ‘</td></tr>’
loading = extract(line,start_marker,end_marker)
if i == 15:
line = text[i]
start_marker = ‘0px;”>’
end_marker = ‘ km/h</td></tr>’
wind = float(extract(line,start_marker,end_marker))
if i == 17:
line = text[i]
start_marker = ‘0px;”>’
end_marker = ‘ km/h</td></tr>’
wind_max = float(extract(line,start_marker,end_marker))
if i == 19:
line = text[i]
start_marker = ‘0px;”>’
end_marker = ‘</td></tr>’
wind_dir = extract(line,start_marker,end_marker)
print time
print “Temperature: ” + str(temp) + ” C”
print “Precipitation since 0600 or 1800: ” + str(precip) + ” mm”
print “Precipitation in last hour: ” + loading
print “Wind speed: ” + str(wind) + ” km/h”
print “Max. wind: ” + str(wind_max) + ” km/h”
print “Wind direction: ” + wind_dir + “\n”
#Get average precipitation
start_marker = ‘ ‘
end_marker = ‘:00:00’
T = extract(time,start_marker,end_marker)
T = int(T)
if T > 6:
T = T – 6
else:
T = T + 6
P = 100*float(precip/T) #mm/hr water equivalent
#Precipitation weight value
if P == 0:
wt_precip = 0
elif 0 < P <= 20:
wt_precip = 6
elif 20 < P <= 40:
wt_precip = 12
elif 40 < P <= 60:
wt_precip = 16
elif 60 < P <= 80:
wt_precip = 24
elif P > 80:
wt_precip = 30
#Assign aspect values to wind direction
if wind_dir == “N”:
aspect = 180
elif wind_dir == “NE”:
aspect = 225
elif wind_dir == “E”:
aspect = 270
elif wind_dir == “SE”:
aspect = 315
elif wind_dir == “S”:
aspect = 0
elif wind_dir == “SW”:
aspect = 45
elif wind_dir == “W”:
aspect = 90
elif wind_dir == “NW”:
aspect = 135
#Aspect weight value
if wind_max == 0:
wt_wind = 2
elif 0 < wind_max <= 5:
wt_wind = 6
elif 5 < wind_max <= 10:
wt_wind = 12
elif 10 < wind_max <= 15:
wt_wind = 16
elif 15 < wind_max <= 20:
wt_wind = 24
elif wind_max > 20:
wt_wind = 30
#Weight for slope
wt_slope = 30
#Elevation bonus
wt_elev = 5
#Immediate loading bonus
if loading == “Yes”:
wt_loading = 5
else:
wt_loading = 0
#ArcGIS processing
#Set environment
arcpy.env.workspace = “C:\data\project\Duffy_backcountry.gdb”
arcpy.env.overwriteOutput = True
arcpy.env.addOutputsToMap = False
#Check out extension
arcpy.CheckOutExtension(“Spatial”)
#Set clipping mask
mask = “bndy_mask”
#Fuzzy Membership for slope, Gaussian Distribution centered at 38 degrees with spread of 0.01
temp_FM_Slope = FuzzyMembership(“Slope_30m”, FuzzyGaussian(38, 0.01))
#Extract by Mask for slope
FM_bndy_Slope = ExtractByMask(temp_FM_Slope, mask)
#Fuzzy Membership for aspect, Gaussian Distribution based on wind direction
if aspect == 0:
#Rotate values 180 degrees and perform FM
temp_Aspect1 = Raster(“Aspect_30m”) + 180
temp_Aspect2 = Con(temp_Aspect1 > 360, temp_Aspect1 – 360, temp_Aspect1)
temp_FM_Aspect = FuzzyMembership(temp_Aspect2, FuzzyGaussian(180, 0.001))
#Extract by Mask for Aspect
FM_bndy_Aspect = ExtractByMask(temp_FM_Aspect, mask)
arcpy.Delete_management(temp_Aspect1)
arcpy.Delete_management(temp_Aspect2)
else:
temp_FM_Aspect = FuzzyMembership(“Aspect_30m”, FuzzyGaussian(aspect, 0.001))
FM_bndy_Aspect = ExtractByMask(temp_FM_Aspect, mask)
#Fuzzy Membership for elevation, need to get local max and min
radius = 100
temp_stat_min = FocalStatistics(“DEM_30m”, NbrCircle(radius, “MAP”), “MINIMUM”, “”)
temp_stat_max = FocalStatistics(“DEM_30m”, NbrCircle(radius, “MAP”), “MAXIMUM”, “”)
temp_FM_Elev = 1.0 * (Raster(“DEM_30m”) – temp_stat_min) / (temp_stat_max – temp_stat_min)
FM_bndy_Elev = ExtractByMask(temp_FM_Elev, mask)
print “Precipitation weight = ” + str(wt_precip) + “%”
print “Slope weight = ” + str(wt_slope) + “%”
print “Aspect weight based on max. wind speed = ” + str(wt_wind) + “%”
print “Elevation bonus = 5%”
print “Immediate loading bonus = ” + str(wt_loading)
#Weighted Sum
WSumTableObj = WSTable([[FM_bndy_Slope, “VALUE”, wt_slope], [FM_bndy_Aspect, “VALUE”, wt_wind],
[FM_bndy_Elev, “VALUE”, wt_elev]])
temp_WSum = WeightedSum(WSumTableObj)
Avalanche_Hazard = temp_WSum + wt_precip + wt_loading
arcpy.env.addOutputsToMap = True
Avalanche_Hazard.save(“Avalanche_Hazard”)
#Delete temp layers
arcpy.Delete_management(temp_FM_Slope)
arcpy.Delete_management(FM_bndy_Slope)
arcpy.Delete_management(temp_FM_Aspect)
arcpy.Delete_management(FM_bndy_Aspect)
arcpy.Delete_management(temp_stat_min)
arcpy.Delete_management(temp_stat_max)
arcpy.Delete_management(temp_FM_Elev)
arcpy.Delete_management(FM_bndy_Elev)
arcpy.Delete_management(temp_WSum)
arcpy.RefreshTOC()
arcpy.RefreshActiveView()
print “Avalanche hazard updated”