%%capture
!pip install geopandas==0.10.2
!pip install mapclassify==2.4.3
!pip install selenium==4.2.0
!pip install pygeos==0.10.2
!pip install git+https://github.com/python-visualization/branca
!pip install -U folium
!apt-get install libspatialindex-dev
!pip install rtree==0.9.7
!pip install -U osmnx==1.1.1
!pip install --upgrade numpy
!pip install matplotlib==3.1.3
from folium import GeoJson, GeoJsonTooltip, Map, Marker, Icon, PolyLine, FeatureGroup
import sys, os, re, mapclassify, datetime, folium, pprint, rtree, requests
from shapely.geometry import LineString, Point
import matplotlib.pyplot as plt
from google.colab import drive
from tqdm.auto import tqdm
import geopandas as gpd
import pandas as pd
import osmnx as ox
import numpy as np
pprint = pprint.PrettyPrinter(indent=4).pprint
drive.mount('/content/drive')
# change this to your folder directory
folder_directory = "/content/drive/MyDrive/Geo-Dolomities"
sys.path.append(folder_directory)
os.chdir(folder_directory)
Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount("/content/drive", force_remount=True).
%%capture
#download the necessary utils files
if 'geo_utils' not in os.listdir():
os.mkdir('geo_utils')
url_utils = "https://raw.githubusercontent.com/GabrieleGhisleni/Geo-Dolomities/master/utils/utils.py"
url_style = "https://raw.githubusercontent.com/GabrieleGhisleni/Geo-Dolomities/master/utils/styles.py"
!wget --no-cache --backups=1 -P'geo_utils' {url_utils}
!wget --no-cache --backups=1 -P'geo_utils' {url_style}
!touch utils/__init__.py
assert 'utils.py' in os.listdir('geo_utils'), 'utils.py not present'
assert 'styles.py' in os.listdir('geo_utils'), 'styles.py not present'
#in case ModuleError restart the kernel
import geo_utils.utils as utils
import geo_utils.styles as styles
tiles = 'Stamen Terrain'
Load the data from https://www.sat.tn.it/sentieri/mappa-sentieri/; this data are available under the license 'Open Data Commons Open Database License (ODbL)' and distributed by "© Società degli Alpinisti Tridentini (SAT)".
Alternatively they were provided (but not updated trough the CatalogueServiceWeb("http://geodati.gov.it/RNDT/csw") under the denomination of 'p_TN_f3547bc8-bf1e-4731-85d2-2084d1f4ba07'.
In this data we have the mountain trail maintaned and handled by SAT.
# load the data and selecting fields
sentieri_SAT = gpd.read_file('data/SAT/sentieri_tratte.shp')
selected_cols = ['numero', 'competenza', 'denominaz', 'difficolta', 'loc_fine', 'quota_iniz', 'quota_fine', 't_andata', 't_ritorno', 'geometry']
sentieri_SAT = sentieri_SAT.loc[:, selected_cols]
sentieri_SAT.to_crs(epsg=4326, inplace=True)
# load dolomities data
dolomiti_df = gpd.read_file('data/geo_dolomiti/dolomiti_geo.shp')
dolomiti_df.loc[2, ["name"]] = "Sistema 3 - Pale di San Martino"
# load italian regions & provinces
dolomiti_regions, dolomiti_provincies, dolomiti_municipalities = utils.load_italian_north_east_area()
# buffering on fly on epsg 32632
dolomiti_df.geometry = dolomiti_df.to_crs(epsg=32632).geometry.buffer(300).to_crs(epsg=4326).values
sentieri_SAT_in_dolomities = sentieri_SAT.geometry.apply(lambda x: utils.element_crosses_or_within_area(x, dolomiti_df, 'name'))
assert dolomiti_df.crs == sentieri_SAT.crs
# selecting only the path that are closeby the dolomities
sentieri_SAT['dolomities'] = sentieri_SAT_in_dolomities
sentieri_SAT_in_dolomities = sentieri_SAT.loc[~sentieri_SAT.dolomities.isna()].reset_index(drop=True)
sentieri_SAT_in_dolomities['DEN_REG'] = sentieri_SAT_in_dolomities.geometry.apply(lambda x: utils.element_within_area(x, dolomiti_regions, 'DEN_REG'))
sentieri_SAT_in_dolomities = sentieri_SAT_in_dolomities.dropna(subset=['DEN_REG'])
sentieri_SAT_in_dolomities = sentieri_SAT_in_dolomities.loc[sentieri_SAT_in_dolomities.geometry.is_valid]
unique_dolomities_group = sentieri_SAT_in_dolomities.dolomities.unique().tolist()
sentieri_SAT_in_dolomities.head(2)
numero | competenza | denominaz | difficolta | loc_fine | quota_iniz | quota_fine | t_andata | t_ritorno | geometry | dolomities | DEN_REG | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | E511 | SEZ. SAT PREDAZZO | VIA FERRATA DEI CAMPANILI DEL LATEMAR | EEA-D | FORCELLA GRANDE - Bivacco "Mario Rigatti" | 2570 | 2634 | 01:30 | 01:30 | LINESTRING (11.56540 46.38048, 11.56547 46.380... | Sistema 7 - Sciliar-Catinaccio & Latemar | Trentino-Alto Adige |
1 | E714 | SAT O.C. | VIA FERRATA "NICO GUSELLA" | EEA-D | FORCELLA DEL PORTON | 2451 | 2420 | 02:00 | 02:00 | LINESTRING (11.84688 46.24668, 11.84683 46.246... | Sistema 3 - Pale di San Martino | Trentino-Alto Adige |
# Explore in which provinces are the SAT paths
path_per_dolomities = {k : len(sentieri_SAT_in_dolomities.loc[sentieri_SAT_in_dolomities.DEN_REG == k]) for k in sentieri_SAT_in_dolomities.DEN_REG.unique()}
for name, value in sorted(path_per_dolomities.items(), key=lambda item: item[1]):
print(f'Number of path in {name}: {value}')
print('-'*50)
# Explore in which group of dolomities are the SAT paths
path_per_dolomities = {k : len(sentieri_SAT_in_dolomities.loc[sentieri_SAT_in_dolomities.dolomities == k]) for k in unique_dolomities_group}
for name, value in sorted(path_per_dolomities.items(), key=lambda item: item[1]):
print(f'Number of path in {name}: {value}')
Number of path in Trentino-Alto Adige: 119 -------------------------------------------------- Number of path in Sistema 2 - Marmolada: 7 Number of path in Sistema 3 - Pale di San Martino: 21 Number of path in Sistema 7 - Sciliar-Catinaccio & Latemar: 29 Number of path in Sistema 9 - Dolomiti di Brenta: 62
# filter all the dolomities only for the intersted provinces:
dolomiti_df = dolomiti_df.loc[dolomiti_df['name'].isin(sentieri_SAT_in_dolomities.dolomities.unique())]
dolomiti_geo = GeoJson(data = dolomiti_df, style_function = styles.style_dolomiti, tooltip = GeoJsonTooltip(fields=['name', 'url_info', 'area']))
sat_dolomities_map = utils.get_new_map(title='SAT Trail on the Dolomities', tiles=tiles)
# create relevant feature layers
feature_layers = {group:FeatureGroup(name=f"Group of dolomities: {group}") for group in unique_dolomities_group}
feature_layers['dolomiti_area'] = FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])
# SAT Paths
for _, row in sentieri_SAT_in_dolomities.iterrows():
line = PolyLine(
locations = [list(reversed(i)) for i in list(row.geometry.coords)],
tooltip = utils.get_info_from_row(row),
popup = utils.get_info_from_row(row),
color='darkred', weight=3,
)
line.add_to(feature_layers[row.dolomities])
utils.add_map_infos(feature_layers, sat_dolomities_map)
# load previusoly processed data about huts
osm_rifugi = gpd.read_file('data/alpine_huts/osm_dolomities.shp')
# filter the huts so that we have the same dolomities group
osm_rifugi = osm_rifugi.loc[osm_rifugi.dolomities.isin(unique_dolomities_group)]
# huts and paths together
sat_hut_dolomities_map = utils.get_new_map(title='Path & Huts according to Dolomities group', tiles=tiles)
feature_layers = {group:FeatureGroup(name=f"Group of dolomities: {group}") for group in unique_dolomities_group}
feature_layers['dolomiti_area'] = FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])
# SAT Paths
for _, row in sentieri_SAT_in_dolomities.iterrows():
line = PolyLine(
locations = [list(reversed(i)) for i in list(row.geometry.coords)],
tooltip = utils.get_info_from_row(row),
popup = utils.get_info_from_row(row),
color='darkred', weight=3,
)
line.add_to(feature_layers[row.dolomities])
# Huts
for idx, row in osm_rifugi.iterrows():
marker = Marker(
location = [i[0] for i in reversed(row.geometry.xy)],
tooltip = utils.get_info_from_row(row),
popup = utils.get_info_from_row(row),
icon = Icon(color = styles.c_colors[unique_dolomities_group.index(row.dolomities)])
)
marker.add_to(feature_layers[row.dolomities])
utils.add_map_infos(feature_layers, sat_hut_dolomities_map, collapsed_legend=True)
osm_rifugi_buffered = osm_rifugi.copy()
osm_rifugi_buffered.geometry = osm_rifugi.to_crs(epsg=32632).buffer(50).to_crs(epsg=4326)
res = []
# filtering only the huts that are reached from a SAT path
for i in range(len(osm_rifugi_buffered)):
intersection_idx = sentieri_SAT_in_dolomities.geometry.crosses(osm_rifugi_buffered.geometry.values[i])
idx_trues = np.nonzero(intersection_idx.values)[0]
if idx_trues.size > 0:
res.append(idx_trues.tolist())
else:
res.append(None)
osm_rifugi_buffered['path'] = res
osm_rifugi_buffered = osm_rifugi_buffered.loc[~osm_rifugi_buffered.path.isna()].reset_index(drop=True)
hut_reachable = osm_rifugi_buffered['name'].unique()
sentieri_SAT_in_dolomities['reachable_huts'] = ""
# add information to the path of the hut that is possible to reach
for idx, row in osm_rifugi_buffered.iterrows():
name = row['name']
for path in row.path:
current = sentieri_SAT_in_dolomities.iloc[path]['reachable_huts']
sentieri_SAT_in_dolomities.loc[path, ['reachable_huts']] = f"{current}, {name}" if current else f"{name}"
# add to huts the trail by which are accessible
trail_to_huts = []
for paths in osm_rifugi_buffered.path:
tmp = []
for path in paths:
tmp.append(sentieri_SAT_in_dolomities.iloc[path].numero)
trail_to_huts.append(', '.join(tmp))
osm_rifugi_buffered["accessible_via_path"] = trail_to_huts
osm_rifugi_buffered['centroid'] = osm_rifugi_buffered.to_crs(epsg=32632).geometry.centroid
sentieri_SAT_in_dolomities = sentieri_SAT_in_dolomities.loc[sentieri_SAT_in_dolomities.geometry.is_valid]
# explore with static maps
fig, axes = plt.subplots(1,1, figsize=(15,15))
dolomiti_df.to_crs(epsg=4326).plot(
ax=axes,
linewidth=2,
edgecolor='black',
column='name',
cmap = "tab10",
legend=True,
categorical=True,
zorder=0,
legend_kwds= {
'bbox_to_anchor':(.3, 1.05),
'fontsize':10,
'frameon':False
}
)
sentieri_SAT_in_dolomities.to_crs(epsg=4326).plot(ax=axes, color = "orange", zorder=1)
osm_rifugi.to_crs(epsg=4326).plot(ax=axes, color='white', markersize=70, marker="X", zorder=2)
plt.title('Huts over the Dolomities area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()
# huts and paths together
sat_hut_path_dolomities_map = utils.get_new_map(title='Path & Huts according to Dolomities group', tiles=tiles)
feature_layers = {group:FeatureGroup(name=f"Hut: {group}") for group in hut_reachable}
feature_layers['dolomiti_area'] = FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])
# Huts - SAT paths
for idx, row in osm_rifugi_buffered.iterrows():
marker_group = row['name']
marker = Marker(
location = [i[0] for i in reversed(row.geometry.centroid.xy)],
tooltip = utils.get_info_from_row(row),
popup = utils.get_info_from_row(row),
icon = Icon(color = styles.c_colors[unique_dolomities_group.index(row.dolomities)])
)
marker.add_to(feature_layers[marker_group])
for path in row.path:
row = sentieri_SAT_in_dolomities.iloc[path]
line = PolyLine(
locations = [list(reversed(i)) for i in list(row.geometry.coords)],
tooltip = utils.get_info_from_row(row),
popup = utils.get_info_from_row(row),
color='darkred', weight=3,
)
line.add_to(feature_layers[marker_group])
utils.add_map_infos(feature_layers, sat_hut_path_dolomities_map, collapsed_legend=True)
Generate the trails network of the main dolomities in Trentino-Alto Adige
import networkx
# for pale di san martino we have some troubles while retreving the trails.
def custom_pale_di_san_martino(dolomiti_df, buffer, network='walk'):
pale = dolomiti_df.loc[dolomiti_df['name'] == 'Sistema 3 - Pale di San Martino']
geom = pale.to_crs(epsg=32632).buffer(buffer).to_crs(epsg=4326).values[0]
e, s, w, n = pale.total_bounds
#generate from bbox rather than polygon
pale_graph = ox.graph_from_bbox(n, s, e, w, network_type=network)
pale_nodes, pale_edges = ox.graph_to_gdfs(pale_graph)
pale_nodes = pale_nodes.loc[pale_nodes.geometry.within(geom)]
pale_edges = pale_edges.loc[pale_edges.geometry.within(geom)]
return ox.graph_from_gdfs(pale_nodes, pale_edges)
# generate trails graph, if already present then just load it.
if 'trail_network.graphml.xml' not in os.listdir('data'):
buffer = 1_500
all_graphs = []
dolomiti_df_buf = dolomiti_df.copy()
dolomiti_df_buf.geometry = dolomiti_df_buf.to_crs(epsg=32632).buffer(buffer).to_crs(epsg=4236)
for idx, row in dolomiti_df_buf.iterrows():
if row['name'] == 'Sistema 3 - Pale di San Martino':
graph = custom_pale_di_san_martino(dolomiti_df, buffer)
else:
graph = ox.graph_from_polygon(
row.geometry,
network_type = 'walk',
)
all_graphs.append(graph)
G_trail = networkx.compose_all(all_graphs)
ox.save_graphml(G_trail, 'data/trail_network.graphml.xml')
else:
G_trail = ox.load_graphml('data/trail_network.graphml.xml')
print(f'total_edges: {len(ox.graph_to_gdfs(G_trail, edges=True, nodes=False))}')
print(f'total_nodes: {len(ox.graph_to_gdfs(G_trail, edges=False, nodes=True))}\n')
#stats
G_trail_stats = ox.basic_stats(G_trail)
pprint(G_trail_stats)
total_edges: 14812 total_nodes: 5804 { 'circuity_avg': 1.309835227220698, 'edge_length_avg': 377.73620969484233, 'edge_length_total': 5595028.738000005, 'intersection_count': 4580, 'k_avg': 5.104066161268091, 'm': 14812, 'n': 5804, 'self_loop_proportion': 0.0005401026194977045, 'street_length_avg': 377.736209694842, 'street_length_total': 2797514.369, 'street_segment_count': 7406, 'streets_per_node_avg': 2.66092350103377, 'streets_per_node_counts': { 0: 0, 1: 1224, 2: 0, 3: 4121, 4: 441, 5: 16, 6: 1, 7: 1}, 'streets_per_node_proportions': { 0: 0.0, 1: 0.21088904203997244, 2: 0.0, 3: 0.7100275671950379, 4: 0.07598208132322536, 5: 0.0027567195037904893, 6: 0.00017229496898690558, 7: 0.00017229496898690558}}
fig, axes = plt.subplots(1,1, figsize=(15,15))
plt.title("Dolomities in Trentino network trails")
dolomiti_df.to_crs(epsg=4326).plot(
ax=axes,
column='name', cmap = "tab10", alpha=0.3,
legend=True, categorical=True,
legend_kwds= {
'bbox_to_anchor':(.3, 1.05),
'fontsize':10,
'frameon':False
}
)
ox.graph_to_gdfs(G_trail, edges=True, nodes=False).plot(ax=axes, color='purple', linewidth=.5)
osm_rifugi.to_crs(epsg=4326).plot(ax=axes, color='orange', markersize=110, marker="X", zorder=3)
plt.tight_layout()
plt.axis('off')
plt.savefig('images/trail.jpg', dpi=300)
plt.show()
trail_maps = utils.get_new_map(title='Trails & Huts according to Dolomities group', tiles='OpenStreetMap')
feature_layers = {group:FeatureGroup(name=f"Group of dolomities: {group}") for group in unique_dolomities_group}
feature_layers['trails'] = FeatureGroup(name='trails')
feature_layers['dolomiti_area'] = FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])
ox.plot_graph_folium(
G_trail,
feature_layers['trails'],
color='darkred',
weight=2
)
for idx, row in osm_rifugi.iterrows():
marker_group = row.dolomities
marker = Marker(
location = [i[0] for i in reversed(row.geometry.centroid.xy)],
tooltip = utils.get_info_from_row(row),
popup = utils.get_info_from_row(row),
icon = Icon(color = styles.c_colors[unique_dolomities_group.index(row.dolomities)])
)
marker.add_to(feature_layers[marker_group])
utils.add_map_infos(feature_layers, trail_maps, collapsed_legend=True)