Inspired and based on the USGS Earthquakes example from the Mapbox GL Jupyter library, this notebook visualizes data on accidents by snow avalanches causing at least one fatality in Switzerland. See Why the Swiss are experts at predicting avalanches (swissinfo) for some interesting background on the subject.
The data used here is sourced from the following open dataset provided by the Swiss Federal Institute for Forest, Snow and Landscape Research WSL, downloadable at envidat.ch. Citation:
SLF (2017): Fatal avalanche accidents in Switzerland since 1995-1996; WSL Institute for Snow and Avalanche Research SLF; doi:10.16904/13.
The CSV file needs to be downloaded and converted to UTF-8, then placed in the same folder as this script under the filename avalancheaccidentsswitzerlandsince1995_utf8.csv
Additionally, State of natural hazard mapping (Cantons): Avalanches (opendata.swiss) gives an overview of the StorMe dataset, which could be explored on the national geoportal in a bit more detail, but is only available on request at the Federal Office for the Environment BAFU, the Institute for Snow and Avalanche Research SLF, or cantonal offices e.g. for the canton of Bern via the Geodata portal - Ereigniskataster. To give an impression of this statistical data, we include a map of it below.
Make sure that aiohttp
, pandas
and pysal
are installed in your environment in addition to the new mapboxgl
library.
import asyncio
from aiohttp import ClientSession
import json, geojson, os, time
import pandas as pd
from datetime import datetime, timedelta
from mapboxgl.viz import *
from mapboxgl.utils import *
import pysal.esda.mapclassify as mapclassify
Set a MAPBOX_ACCESS_TOKEN
environment variable, or copy your token to use this notebook.
token = os.getenv('MAPBOX_ACCESS_TOKEN')
The following conversion function for Swiss coordinates (CH1903) to World geodetic system (WGS1984) was adapted from a script by Valentin Minder.
def CHtoWGSlat(y, x):
# Axiliary values (% Bern)
y_aux = (y - 600000) / 1000000
x_aux = (x - 200000) / 1000000
lat = (16.9023892 + (3.238272 * x_aux)) + \
- (0.270978 * pow(y_aux, 2)) + \
- (0.002528 * pow(x_aux, 2)) + \
- (0.0447 * pow(y_aux, 2) * x_aux) + \
- (0.0140 * pow(x_aux, 3))
# Unit 10000" to 1" and convert seconds to degrees (dec)
lat = (lat * 100) / 36
lng = (2.6779094 + (4.728982 * y_aux) + \
+ (0.791484 * y_aux * x_aux) + \
+ (0.1306 * y_aux * pow(x_aux, 2))) + \
- (0.0436 * pow(y_aux, 3))
# Unit 10000" to 1" and convert seconds to degrees (dec)
lng = (lng * 100) / 36
return (lat, lng)
# Quick test: result should approx. (46.42, 7.12)
CHtoWGSlat(575300,141010)
(46.419983053690665, 7.117350584670833)
Here we access the CSV file obtained from envidat.ch, and converted to UTF-8 using LibreOffice:
df = pd.read_csv('avalancheaccidentsswitzerlandsince1995_utf8.csv')
df.head(5)
avalanche.id | date | date.quality | hydrological.year | canton | local.name | start.zone.coordinates.x | start.zone.coordinates.y | coordinates.quality | start.zone.elevation | start.zone.slope.aspect | start.zone.inclination | forecasted.dangerlevel | number.dead | number.caught | number.fully.buried | activity | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 13007 | 1995-12-25 | 0 | 1995/96 | VS | Chetseron / Vallon de l`Ertentse | 602380 | 131230 | 25 | 2060 | NW | NaN | NaN | 2 | 2 | 0 | offpiste |
1 | 13014 | 1995-12-28 | 0 | 1995/96 | VS | Verbier / Les Ruinettes | 586070 | 104160 | 25 | 2340 | NNE | 39.0 | NaN | 1 | 1 | 1 | offpiste |
2 | 13028 | 1996-01-14 | 0 | 1995/96 | GR | Chilchalphorn | 731720 | 155110 | 25 | 2960 | E | 35.0 | NaN | 1 | 4 | 2 | tour |
3 | 13038 | 1996-02-14 | 0 | 1995/96 | VD | La Lécherette | 575300 | 141010 | 25 | 1580 | NNW | 42.0 | NaN | 1 | 2 | 1 | offpiste |
4 | 13040 | 1996-02-15 | 0 | 1995/96 | OW | Pilatus / Matthorn / Ruessiflue | 661640 | 202180 | 25 | 1860 | NW | 40.0 | NaN | 1 | 4 | 3 | tour |
We create a toll
column by adding the counts of people affected by the avalance, as well as an elevation
and geoposition columns using our conversion function.
df['toll'] = df['number.dead'] + df['number.caught'] + df['number.fully.buried']
df['elevation'] = df['start.zone.elevation']
for i in df.index:
x = df.at[i, 'start.zone.coordinates.x']
y = df.at[i, 'start.zone.coordinates.y']
lat, lon = CHtoWGSlat(x, y)
df.at[i, 'latitude'] = lat
df.at[i, 'longitude'] = lon
dfc = df[['date', 'toll', 'elevation', 'latitude', 'longitude']]
print(dfc.shape)
dfc.head(3)
(381, 5)
date | toll | elevation | latitude | longitude | |
---|---|---|---|---|---|
0 | 1995-12-25 | 4 | 2060 | 46.332456 | 7.469545 |
1 | 1995-12-28 | 3 | 2340 | 46.088807 | 7.258541 |
2 | 1996-01-14 | 7 | 2960 | 46.534313 | 9.155739 |
min(df['elevation']), max(df['elevation'])
(725, 4000)
min(df['toll']), max(df['toll'])
(2, 31)
Now a GeoJSON file is generated from the DataFrame using a Mapbox utility.
df_to_geojson(df, filename='points.geojson', precision=4,
lon='longitude', lat='latitude',
properties=['toll','elevation','date'])
{'feature_count': 381, 'filename': 'points.geojson', 'type': 'file'}
# Calculate gradients and sizes
color_breaks = mapclassify.Natural_Breaks(df['toll'], k=6, initial=0).bins
color_stops = create_color_stops(color_breaks, colors='YlOrRd')
radius_breaks = mapclassify.Natural_Breaks(df['elevation'], k=10, initial=0).bins
#radius_breaks = flip(radius_breaks, 0) # causes a MapBox error (strictly ascending order)
radius_stops = create_radius_stops(radius_breaks, 1, 10)
# Create the visualization
viz = GraduatedCircleViz('points.geojson',
color_property = 'toll',
color_stops = color_stops,
radius_property = 'elevation',
radius_stops = radius_stops,
opacity=0.8,
zoom=6.5,
center= [8.4, 46.8],
access_token=token)
viz.show()
For comparison, here is the abovementioned national dataset of avalanche risks, displayed on the Swiss geoportal (direct link):
from IPython.display import IFrame
IFrame("https://map.geo.admin.ch/embed.html?lang=en&topic=ech&bgLayer=ch.swisstopo.pixelkarte-farbe&layers=ch.swisstopo.zeitreihen,ch.bfs.gebaeude_wohnungs_register,ch.bav.haltestellen-oev,ch.swisstopo.swisstlm3d-wanderwege,ch.bafu.showme-gemeinden_sturzprozesse&layers_visibility=false,false,false,false,true&layers_timestamp=18641231,,,,&E=2663226.61&N=1193507.44&zoom=1&layers_opacity=1,1,1,1,0.75",
width='100%', height=500)
© Oleg Lavrovsky, February 2018. This notebook is licensed under a Creative Commons Attribution 4.0 International License.