Table of contents

• Environment set up
• OSM OverPass API Turbo for Alpine huts
  • Load the alpine huts data
  • Alpine huts provinces
  • Alpine huts elevation
  • SAT huts prices
• Dolomities huts
  • Dolomities huts elevation

• OSM OverPass API Turbo for Peas

  • Load the peaks data
  • Dolomities Peaks
  • Peaks' 3000m
  • SAT huts prices

• Dolomities Huts and Peaks

  • Huts close to peaks
  • Absolute distances

@Gabriele Ghisleni, Data Science.

---

  0. Environment set up

Install the required dependecies on google colab and mount the drive.

%%capture
!pip install geopandas==0.10.2
!pip install mapclassify==2.4.3
!pip install pygeos==0.10.2
!pip install rasterio==1.2.3
!pip install owslib==0.25.0
!pip install -U folium

import warnings, sys, os, folium, owslib, rasterio, geopy, shapely
from folium import FeatureGroup, GeoJson, Marker, Icon
import matplotlib.pyplot as plt
from google.colab import drive
from tqdm.auto import tqdm
import geopandas as gpd
import pandas as pd
import numpy as np

warnings.simplefilter("ignore")
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)

%%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/geo_utils/utils.py"
  url_style = "https://raw.githubusercontent.com/GabrieleGhisleni/Geo-Dolomities/master/geo_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'

  1. OSM OverPass API Turbo for Alpine huts

First we will extract all the nodes marked as 'alpine_hut' from OSM for having a available more data about Veneto, Trentino-Alto Adige and Friuli-Venezia Giulia.

Thanks to Overpass api we could query from geocodeArea, this is the reason why we manually download them as GeoJson.

@https://overpass-turbo.eu/

[out:json][timeout:25];
{{geocodeArea:Veneto}}->.searchArea;
(
  node["tourism"="alpine_hut"](area.searchArea);
  way["tourism"="alpine_hut"](area.searchArea);
);

out body;

Doing this we will be able to retrieve more alpine huts rather than rely only on the CAI resources.

  1.1 Load the alpine huts data

osm_rifugi = pd.concat([
  gpd.read_file('data/alpine_huts/bolzano_rifugi_osm.geojson'),
  gpd.read_file('data/alpine_huts/trento_rifugi_osm.geojson'),
  gpd.read_file('data/alpine_huts/friuli_rifugi_osm.geojson'),
  gpd.read_file('data/alpine_huts/veneto_rifugi_osm.geojson')                       
])

# load previusly prepared gdf
dolomiti_df = gpd.read_file('data/geo_dolomiti/dolomiti_geo.shp')
dolomiti_df.loc[2, ["name"]] = "Sistema 3 - Pale di San Martino"

# create the folium geojson for next rendering
dolomiti_geo = GeoJson(
    data = dolomiti_df,
    style_function = styles.style_dolomiti,
    tooltip = folium.GeoJsonTooltip(fields=['name', 'url_info', 'area'])
  )

# load italian geospatial data
dolomiti_regions, dolomiti_provincies, dolomiti_municipalities = utils.load_italian_north_east_area()

# Some preprocessing
osm_rifugi.replace({np.nan: None}, inplace=True)
osm_rifugi.dropna(subset=['name'], inplace=True)

selection = ['addr:city', 
             'name',
             'ele', 
             'operator', 
             'opening_hours',
             'capacity', 
             'contact:email', 
             'contact:mobile', 
             'website',
             'geometry',
            ]

osm_rifugi = osm_rifugi.loc[:, selection]

# removing 'malga'
osm_rifugi = osm_rifugi.loc[osm_rifugi.name.apply(
    lambda x: ('hütte' in x.lower() or 'rifugio' in x.lower()) and 'malga' not in x.lower()
  )]

# be sure that all the hut are marked as points and not multipolygon
osm_rifugi.geometry = osm_rifugi.geometry.apply(utils.from_poly_to_point)
osm_rifugi.reset_index(drop=True, inplace=True)

# set proper type of 'ele' which stand for the elevation on the sea level
osm_rifugi.ele = osm_rifugi.ele.str.replace('m', '')
osm_rifugi.ele.fillna(0, inplace=True)
osm_rifugi.ele = osm_rifugi.ele.astype(float)

# set the province
assert dolomiti_provincies.crs == osm_rifugi.crs
osm_rifugi['DEN_PROV'] = osm_rifugi.geometry.apply(lambda x: utils.element_within_area(x, dolomiti_provincies, 'DEN_PROV'))
osm_rifugi = osm_rifugi[~osm_rifugi.DEN_PROV.isna()]

unique_provice_dolomities = osm_rifugi.DEN_PROV.unique().tolist()

if not 'osm_rifugi_processed.shp' is os.listdir('data/alpine_huts'):
  osm_rifugi.to_file('data//alpine_huts/osm_rifugi_processed.shp')

  1.2 Alpine huts provinces

First explorations to understand where the majority of the huts are located.

fig, axes = plt.subplots(1,1, figsize=(15,15))

dolomiti_regions.to_crs(epsg=4326).plot(ax=axes, cmap="Pastel1_r", column='DEN_REG')
osm_rifugi.to_crs(epsg=4326).plot(
    ax=axes, 
    cmap='tab10', 
    markersize=70,
    marker='H',
    column='DEN_PROV', 
    legend=True
  )



plt.title('Alpine huts divided by province', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()

# explore number of huts for each provinces
hut_per_province = {k : len(osm_rifugi.loc[osm_rifugi.DEN_PROV == k]) for k in osm_rifugi.DEN_PROV.dropna().unique()}
for name, value in reversed(sorted(hut_per_province.items(), key=lambda item: item[1])):
  print(f'Number of Hut in {name}: {value}')

Number of Hut in Trento: 136
Number of Hut in Bolzano: 117
Number of Hut in Belluno: 103
Number of Hut in Udine: 35
Number of Hut in Vicenza: 29
Number of Hut in Verona: 19
Number of Hut in Pordenone: 11
Number of Hut in Treviso: 4
Number of Hut in Trieste: 1

Since they are many we will try to create an interactive map with MarkerCluster for better visualize them.

#  https://stackoverflow.com/questions/37466683/create-a-legend-on-a-folium-map
from folium.plugins import MarkerCluster

osm_hut_map = utils.get_new_map(title='Hut divided by region of provenience', tiles=tiles)
colors = ['gray', 'lightgreen', 'green', 'orange', 'darkred','black']

feature_layers = {
    group: FeatureGroup(name=f"<span style='color:{styles.c_colors[idx]}'>Hut in {group}</span>") 
    for idx, group in enumerate(unique_provice_dolomities)
  }

marker_clusters_ops = {
    'zoomToBoundsOnClick': False, 'maxClusterRadius': 110, 'spiderfyOnMaxZoom': 2, 'disableClusteringAtZoom': 10
  }

marker_clusters = {
    group: MarkerCluster(
        options = marker_clusters_ops,
        name = f"<span style='color:{styles.c_colors[idx]}'>Hut in {group}</span>"
      ) 
    for idx, group in enumerate(unique_provice_dolomities)
  }

feature_layers['dolomiti_area'] =  FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])

# add huts to the map
for idx, row in osm_rifugi.iterrows():
  marker_group = row.DEN_PROV
  marker = Marker(
      location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row), 
      tooltip = utils.get_info_from_row(row),
      icon = Icon(color = styles.c_colors[unique_provice_dolomities.index(row.DEN_PROV)])
  )
  marker.add_to(marker_clusters[marker_group])

for k,v in marker_clusters.items():
  marker_clusters[k].add_to(feature_layers[k])

utils.add_map_infos(feature_layers, osm_hut_map, collapsed_legend=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

  1.3 Alpine huts elevation

Explore the elevation of the huts.

osm_hut_dolomities_map_ele = utils.get_new_map(title='Hut colored according to the elevation', tiles=tiles)

altitudes = [500, 1000, 1500, 2000, 2500, 2800]
colors = ['gray', 'lightgreen', 'green', 'orange', 'darkred','black']

feature_layers = {
    group:FeatureGroup(name=f"<span style='color:{colors[idx]}'>Hut higher than {group}m</span>") 
    for idx, group in enumerate(altitudes)
  }

feature_layers['dolomiti_area'] =  FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])

# huts according to the elevation
for idx, row in osm_rifugi.iterrows():
  if row.ele != 0:
    color, marker_group = styles.hut_style_ele(row.ele)
    marker = Marker(
      location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row), 
      tooltip = utils.get_info_from_row(row),
      icon = Icon(color = color)
    )
    marker.add_to(feature_layers[marker_group])

utils.add_map_infos(feature_layers, osm_hut_dolomities_map_ele, collapsed_legend=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

fig, axes = plt.subplots(1,1, figsize=(15,15))

dolomiti_regions.to_crs(epsg=4326).plot(ax=axes, cmap="Pastel1_r", column='DEN_REG')

osm_rifugi.to_crs(epsg=4326).plot(
    ax=axes, 
    cmap='magma', 
    markersize=70,
    column='ele', 
    legend=True,
    scheme='JenksCaspall',
    marker='X'
  )
    
plt.title('Alpine huts according to the elevation', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()

  1.4 SAT huts prices

# prices manually collected from the SAT information on prices
huts_prices = pd.read_csv('data/huts_prices.csv', encoding='latin-1',  on_bad_lines='skip', header=[0], sep=";")
sat_rifugi = gpd.read_file("data/alpine_huts/rifugi_cai.shp") 


# create common columns to merge the datasets
idxs = []
for _, row in huts_prices.iterrows():
  for idx, sat in sat_rifugi.iterrows():
    if row["name"] in sat.display_na:
      idxs.append(idx)

# here we miss many alpine huts compared to osm data.
huts_prices.index = idxs
huts_prices_merged = huts_prices.merge(sat_rifugi, left_index=True, right_index=True).reset_index(drop=True)
huts_prices_merged = gpd.GeoDataFrame(huts_prices_merged)
huts_prices_merged = huts_prices_merged.to_crs(epsg=4326)
huts_prices_merged.tail(2)

	name	elevation	categoria	posto_riposo_emergenze	posto_letto_soci	posto_letto_non_soci	mezza_pensione_soci	mezza_pensione_non_soci	acqua	municipali	...	extended_n	opening_ho	contact_mo	winter_roo	capacity	website	email	display_na	type_descr	geometry
26	Rifugio XII Apostoli	2489	d	5.0	17	34	48	58	NaN	San Lorenzo Dorsino	...	Rifugio XII Apostoli	Jun 20-Sep 20	+39 339 8405054	yes	38	http://www.dodiciapostoli.it/	rifugio@dodiciapostoli.it	Rifugio XII Apostoli, Via ferrata sentiero Bre...	alpine_hut	POINT (10.84743 46.15064)
27	Rifugio Vajolet	2243	c	5.0	16	32	47	58	NaN	San Giovanni di Fassa	...	Rifugio Vajolet	Jun 20-Sep 20	+39 335 7073258	yes	130	https://www.rifugiovajolet.com/	info@rifugiovajolet.com	Rifugio Vajolet, 546, San Giovanni di Fassa, C...	alpine_hut	POINT (11.63276 46.45860)

2 rows × 23 columns

# summarinzing the prices information
cost_indicator_field = [
                        "posto_riposo_emergenze",	
                        "posto_letto_soci",	
                        "posto_letto_non_soci",	
                        "mezza_pensione_soci",	
                        "mezza_pensione_non_soci"
                      ]
cost_indicator = []

for _, row in huts_prices.iterrows():
  cost_indicator.append(
      (row[cost_indicator_field].sum())
  )

huts_prices_merged['cost_indicator'] = cost_indicator

osm_hut_dolomities_map_ele = utils.get_new_map(title='Hut colored according to the elevation with prices information', tiles=tiles)

altitudes = [500, 1000, 1500, 2000, 2500, 2800]
colors = ['gray', 'lightgreen', 'green', 'orange', 'darkred','black']

feature_layers = {
    group:FeatureGroup(name=f"<span style='color:{colors[idx]}'>Hut higher than {group}m</span>") 
    for idx, group in enumerate(altitudes)
  }

feature_layers['dolomiti_area'] =  FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])

# hut in dolomities according to the elevation
for idx, row in huts_prices_merged.loc[:, ['name', 'elevation', 'categoria', 'extended_n', 'geometry', 'cost_indicator'] + cost_indicator_field].iterrows():
  color, marker_group = styles.hut_style_ele(row.elevation)
  marker = Marker(
      location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row), 
      tooltip = utils.get_info_from_row(row),
      icon=folium.Icon(color = color)
  )
  
  marker.add_to(feature_layers[marker_group])

utils.add_map_infos(feature_layers, osm_hut_dolomities_map_ele, collapsed_legend=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

  2. Dolomities huts

Select only the huts that are nearby the dolomites

buffered_dolomiti = dolomiti_df.copy()

# buffering on fly on epsg 32632
buffered_dolomiti.geometry = dolomiti_df.to_crs(epsg=32632).geometry.buffer(500).to_crs(epsg=4326).values
hut_in_dolomities_buffer = osm_rifugi.geometry.apply(lambda x: utils.element_within_area(x, buffered_dolomiti, 'name'))
assert buffered_dolomiti.crs == osm_rifugi.crs

# selecting only the huts that are closeby the dolomities
osm_rifugi['dolomities'] = hut_in_dolomities_buffer
osm_dolomites = osm_rifugi.loc[~osm_rifugi.dolomities.isna()].reset_index(drop=True)

m_km = lambda x: round(int(x.split(' ')[0]) * 1e-6)
# exploring the huts for each group of dolomities
hut_per_dolomities = {k : len(osm_dolomites.loc[osm_dolomites.dolomities == k]) for k in osm_dolomites.dolomities.unique()}
for name, value in reversed(sorted(hut_per_dolomities.items(), key=lambda item: item[1])):
  print(f"Number of Hut in{name.split('-')[1]}: {value}, total territory of respectively dolomies group: {m_km(dolomiti_df[dolomiti_df['name'] == name]['area'].values[0])} km^2")

Number of Hut in Dolomiti settentrionali: 36, total territory of respectively dolomies group: 536 km^2
Number of Hut in Pale di San Martino: 19, total territory of respectively dolomies group: 317 km^2
Number of Hut in Sciliar: 18, total territory of respectively dolomies group: 93 km^2
Number of Hut in Dolomiti di Brenta: 9, total territory of respectively dolomies group: 111 km^2
Number of Hut in Dolomiti friulane e d'Oltre Piave: 7, total territory of respectively dolomies group: 215 km^2
Number of Hut in Pelmo & Croda da Lago: 5, total territory of respectively dolomies group: 43 km^2
Number of Hut in Marmolada: 2, total territory of respectively dolomies group: 22 km^2
Number of Hut in Puez: 2, total territory of respectively dolomies group: 79 km^2

print(f'Number of hut in the dolomities: {len(osm_dolomites)}, total hut: {len(osm_rifugi)}')
osm_hut_dolomities_map = utils.get_new_map(title='Hut divided by Dolomities group', tiles=tiles)
unique_dolomities_group = osm_dolomites.dolomities.unique().tolist()

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'])

# hut in the dolomities
for idx, row in osm_dolomites.iterrows():
  marker_group = row.dolomities
  marker = Marker(
      location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row), 
      tooltip = 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, osm_hut_dolomities_map, collapsed_legend=True)

Number of hut in the dolomities: 98, total hut: 455

Make this Notebook Trusted to load map: File -> Trust Notebook

# explore with static maps
fig, axes = plt.subplots(1,1, figsize=(15,15))

dolomiti_df.to_crs(epsg=4326).plot(
    ax=axes,
    column='name',
    cmap = "Set2",
    legend=True,  
    categorical=True,

    legend_kwds= {
        'bbox_to_anchor':(.3, 1.05),
        'fontsize':10,
        'frameon':False
        }
)

osm_dolomites.to_crs(epsg=4326).plot(ax=axes, color = "red", markersize=70, marker='X')

plt.title('Alpine Huts in the Dolomities Area', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()

  2.1 Dolomities huts elevation

osm_hut_dolomities_map_ele = utils.get_new_map(title='Hut colored according to the elevation', tiles=tiles)

altitudes = [500, 1000, 1500, 2000, 2500, 2800]
colors = ['gray', 'lightgreen', 'green', 'orange', 'darkred','black']

feature_layers = {
    group:FeatureGroup(name=f"<span style='color:{colors[idx]}'>Hut higher than {group}m</span>") 
    for idx, group in enumerate(altitudes)
  }

feature_layers['dolomiti_area'] =  FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])

# hut in dolomities according to the elevation
for idx, row in osm_dolomites.iterrows():
  color, marker_group = styles.hut_style_ele(row.ele)
  marker = Marker(
      location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row), 
      tooltip = utils.get_info_from_row(row),
      icon=folium.Icon(color = color)
  )
  
  marker.add_to(feature_layers[marker_group])

utils.add_map_infos(feature_layers, osm_hut_dolomities_map_ele, collapsed_legend=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

  3. OSM OverPass API for Peaks

Similar to before we are now able to retrieve the peaks in the selected regions.

query:

[out:json][timeout:25];
{{geocodeArea:Trento}}->.searchArea;
(
  // query part for: “peak”
  node["natural"="peak"](area.searchArea);
);
out body;
>;
out skel qt;

  3.1 Load the peaks data

osm_peak = pd.concat([
  gpd.read_file('data/peaks/osm_peak_bolzano.geojson'),
  gpd.read_file('data/peaks/osm_peak_trento.geojson'),
  gpd.read_file('data/peaks/osm_peak_veneto.geojson'),                      
  gpd.read_file('data/peaks/osm_peak_friuli.geojson')                        
])

# some preprocess
osm_peak.drop(columns=['id', '@id'], inplace=True)
osm_peak.replace({np.nan: None}, inplace=True)
osm_peak.dropna(subset=['ele'], inplace=True)
osm_peak.ele = osm_peak.ele.astype(float)
osm_peak.head(2)

	alt_name	alt_name:de	alt_name:it	alt_name:lld	alt_name:vi	alt_name_2	artist_name	comment	description	direction	...	name:es	name:sl	name:slo	notes	old_name:it	old_name:sl	old_source:name	population	short_name:fur	summit:bell
0	None	None	None	None	None	None	None	None	None	None	...	None	None	None	None	None	None	None	None	None	None
1	None	None	None	None	None	None	None	None	None	None	...	None	None	None	None	None	None	None	None	None	None

2 rows × 102 columns

  3.2 Dolomities Peaks

#check peaks in dolomities
assert osm_peak.crs == buffered_dolomiti.crs
peak_in_dolomities_buffer = osm_peak.geometry.apply(lambda x: utils.element_within_area(x, buffered_dolomiti, 'name'))

osm_peak['dolomities'] = peak_in_dolomities_buffer
osm_dolomites_peak = osm_peak.loc[~osm_peak.dolomities.isna()].reset_index(drop=True)

# explore with static maps
fig, axes = plt.subplots(1,1, figsize=(15,15))


dolomiti_df.to_crs(epsg=4326).plot(
    ax=axes,
    column='name',
    cmap = "Set2",
    legend=True,  
    categorical=True,
    legend_kwds= {
        'bbox_to_anchor':(.3, 1.05),
        'fontsize':10,
        'frameon':False
        }
  )

osm_dolomites_peak.to_crs(epsg=4326).plot(ax=axes, color='red', markersize=40, marker=10)

plt.title('Peaks in Dolomities', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()

# number of peaks directly related to the extensions.
for _, row in osm_dolomites_peak.groupby('dolomities').size().sort_values(0, ascending=False).reset_index().iterrows():
  print(f"'{row.dolomities}' has {row[0]} peaks; total territory of respectively dolomies group: {m_km(dolomiti_df[dolomiti_df['name'] == row.dolomities]['area'].values[0])} km^2")

'Sistema 5 - Dolomiti settentrionali' has 342 peaks; total territory of respectively dolomies group: 536 km^2
'Sistema 3 - Pale di San Martino' has 271 peaks; total territory of respectively dolomies group: 317 km^2
'Sistema 4 - Dolomiti friulane e d'Oltre Piave' has 262 peaks; total territory of respectively dolomies group: 215 km^2
'Sistema 9 - Dolomiti di Brenta' has 113 peaks; total territory of respectively dolomies group: 111 km^2
'Sistema 7 - Sciliar-Catinaccio & Latemar' has 104 peaks; total territory of respectively dolomies group: 93 km^2
'Sistema 6 - Puez-Odle' has 72 peaks; total territory of respectively dolomies group: 79 km^2
'Sistema 1 - Pelmo & Croda da Lago' has 38 peaks; total territory of respectively dolomies group: 43 km^2
'Sistema 2 - Marmolada' has 27 peaks; total territory of respectively dolomies group: 22 km^2
'Sistema 8 - Bletterbach' has 1 peaks; total territory of respectively dolomies group: 3 km^2

osm_peak_dolomities_map = utils.get_new_map(title='Peak colored according to the elevation', tiles=tiles)

altitudes = [500, 1000, 1500, 2000, 2500, 2800]
colors = ['gray', 'lightgreen', 'green', 'orange', 'darkred','black']

feature_layers = {
    group: FeatureGroup(name=f"<span style='color:{colors[idx]}'>Hut higher than {group}m</span>") 
    for idx, group in enumerate(altitudes)
  }

marker_clusters = {
    group: MarkerCluster(name=f"<span style='color:{colors[idx]}'>Hut higher than {group}m</span>") 
    for idx, group in enumerate(altitudes)
  }

feature_layers['dolomiti_area'] =  FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])

# all peaks according to the elevation
for idx, row in osm_dolomites_peak.iterrows():
  color, marker_group = styles.hut_style_ele(row.ele)
  marker = Marker(
      location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row), 
      tooltip = utils.get_info_from_row(row),
      icon=folium.Icon(color=color)
  )

  marker.add_to(marker_clusters[marker_group])

for k,v in marker_clusters.items():
  marker_clusters[k].add_to(feature_layers[k])

utils.add_map_infos(feature_layers, osm_peak_dolomities_map, collapsed_legend=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

  3.3 Peaks' 3000m

Let's inspect the most intersting peaks, the over 3000 meters peaks.

# very high peaks
osm_peak_dolomities_map = utils.get_new_map(title='Dolomities very high peak > 3000', tiles=tiles)

altitudes = [3000, 3100, 3200, 3300]
colors = ['orange', 'red', 'darkred','black']

feature_layers = {
    group:FeatureGroup(name=f"<span style='color:{colors[idx]}'>Hut higher than {group}m</span>") 
    for idx, group in enumerate(altitudes)
  }

feature_layers['dolomiti_area'] =  FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])

for idx, row in osm_dolomites_peak.iterrows():
  if row.ele > 3000:
    color, marker_group = styles.hut_style_ele_high_peak(row.ele)
    marker = Marker(
        location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row), 
      tooltip = utils.get_info_from_row(row),
        icon=folium.Icon(color=color)
    )
    
    marker.add_to(feature_layers[marker_group])

utils.add_map_infos(feature_layers, osm_peak_dolomities_map, collapsed_legend=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

# municipality for each region that touches the dolomities
to_iter = osm_dolomites_peak.loc[osm_dolomites_peak.ele > 3000].groupby('dolomities').size().sort_values(0, ascending=False)
for _, row in to_iter.reset_index().iterrows():
  print(f"'{row.dolomities}' has {row[0]} high peaks")

'Sistema 5 - Dolomiti settentrionali' has 24 high peaks
'Sistema 3 - Pale di San Martino' has 10 high peaks
'Sistema 2 - Marmolada' has 8 high peaks
'Sistema 9 - Dolomiti di Brenta' has 8 high peaks
'Sistema 1 - Pelmo & Croda da Lago' has 3 high peaks
'Sistema 6 - Puez-Odle' has 2 high peaks
'Sistema 7 - Sciliar-Catinaccio & Latemar' has 1 high peaks

  4. Dolomities Huts and Peaks

Let's try to understand the Hut that are very close to some peaks

  4.1 Huts close to Peaks

# https://github.com/napo/geospatial_course_unitn/blob/master/code/solutions/02_exercise_spatial_relationship_and_operations.ipynb
from shapely.ops import nearest_points

def get_closest_point(peak, hut_unary, hut_dataset):
  _, hut = nearest_points(peak, hut_unary)
  return hut_dataset.loc[hut_dataset.geometry == hut].name.values[0]

hut_unary = osm_dolomites.geometry.unary_union
assert osm_dolomites_peak.crs == osm_dolomites.crs

osm_dolomites_high_peak = osm_dolomites_peak.copy()
osm_dolomites_high_peak = osm_dolomites_high_peak.loc[osm_dolomites_high_peak.ele > 3000]
osm_dolomites_high_peak['closest_hut'] = osm_dolomites_high_peak.geometry.apply(lambda x: get_closest_point(x, hut_unary, osm_dolomites))

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 = "Set2",
    legend=True,  
    categorical=True,
    legend_kwds= {
        'bbox_to_anchor':(.3, 1.05),
        'fontsize':10,
        'frameon':False
        }
)

osm_dolomites_high_peak.to_crs(epsg=4326).plot(ax=axes, color = "red", markersize=40, marker=10)
osm_dolomites.to_crs(epsg=4326).plot(ax=axes, color = "blue", markersize=70, marker="X")

plt.title('Peaks-Huts in Dolomities', fontdict=dict(size=25))
plt.tight_layout()
plt.axis('off')
plt.show()

for idx, row in osm_dolomites_high_peak.iterrows():
  ele_hut = osm_dolomites.loc[osm_dolomites['name'] == row['closest_hut']].ele.values[0]
  group_1 = osm_dolomites.loc[osm_dolomites['name'] == row['closest_hut']].dolomities.values[0] 
  group_2 = row.dolomities
  same_group = group_1 and group_2
  if ele_hut and same_group:
    diff = row.ele - ele_hut
    s = f"Peak '{row['name']}'"
    s += " "*(50 - len(s))
    s += f"-> '{row['closest_hut']}'"
    s += " "*(110 - len(s))
    s += f"-> elevation from hut: {int(diff)}m"
    print(s)

osm_peak_hut_dolomities_map = utils.get_new_map(title='Hut close to the highest Peaks', tiles=tiles)

groups = ['peak', 'hut']
colors = ['darkred', 'blue']

feature_layers = {
    group:FeatureGroup(name=f"<span style='color:{colors[idx]}'>{group.capitalize()}</span>") 
    for idx, group in enumerate(groups)
  }

feature_layers['dolomiti_area'] =  FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])

# high peak
for idx, row in osm_dolomites_high_peak.iterrows():
  marker = Marker(
      location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row), 
      tooltip = utils.get_info_from_row(row),
      icon=folium.Icon(color='darkred')
  )

  marker.add_to(feature_layers['peak'])

# hut closesest to the highest peaks
set_uniques = osm_dolomites_high_peak['closest_hut'].unique()
for idx, row in osm_dolomites.loc[osm_dolomites.name.isin(set_uniques)].iterrows():
  marker = Marker(
      location = [i[0] for i in reversed(row.geometry.xy)],
      popup = utils.get_info_from_row(row), 
      tooltip = utils.get_info_from_row(row),
      icon=folium.Icon(color='blue')
  )
  marker.add_to(feature_layers['hut'])

utils.add_map_infos(feature_layers, osm_peak_hut_dolomities_map, collapsed_legend=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

  4.2 Absolute distances

compute the absolute distance between the peaks and its closests hut. Moreover we define the shortest absolute path between the peak and the hut (straight line from one point to the other).

We will perform this operation better in the next notebook where we will take in count all trails as a network.

res = []
osm_dolomites_32632 = osm_dolomites.to_crs(epsg=32632)
osm_dolomites_high_peak.reset_index(inplace=True, drop=True)

# find the distance in absolute terms
assert osm_dolomites_32632.crs == osm_dolomites_high_peak.to_crs(epsg=32632).crs
for idx, row in osm_dolomites_high_peak.to_crs(epsg=32632).iterrows():
  peak = row["geometry"]
  hut = osm_dolomites_32632.loc[osm_dolomites_32632['name'] == row["closest_hut"]].geometry.values[0]
  res.append({'distance' : round(peak.distance(hut))})

# find the absolute path from a -> b (for very proficient trekkers :P).
assert osm_dolomites_high_peak.crs == osm_dolomites.crs
for idx, row in osm_dolomites_high_peak.iterrows():
  peak = row["geometry"]
  hut = osm_dolomites.loc[osm_dolomites['name'] == row["closest_hut"]].geometry.values[0]
  ele_hut = osm_dolomites.loc[osm_dolomites['name'] == row['closest_hut']].ele.values[0]
  res[idx]['geometry'] = shapely.geometry.LineString([peak,hut])
  
  if ele_hut:
    res[idx]['difference'] = - (ele_hut - row.ele)
  
        
linestring_distance = gpd.GeoDataFrame(res, crs="EPSG:4326")

osm_peak_hut_dolomities_map = utils.get_new_map(title='Shortest Path between peaks and huts', tiles=tiles)
groups = ['peak', 'hut', 'peak_hut_path']
colors = ['darkred', 'blue', 'black']

feature_layers = {
    group:FeatureGroup(name=f"<span style='color:{colors[idx]}'>{group.capitalize()}</span>") 
    for idx, group in enumerate(groups)
  }

feature_layers['dolomiti_area'] =  FeatureGroup(name='Area Dolomiti')
dolomiti_geo.add_to(feature_layers['dolomiti_area'])

# peaks
for idx, row in osm_dolomites_high_peak.iterrows():
  marker = Marker(
    location = [i[0] for i in reversed(row.geometry.xy)],
    popup = utils.get_info_from_row(row), 
    tooltip = utils.get_info_from_row(row),
    icon=folium.Icon(color='darkred')
  )
  marker.add_to(feature_layers['peak'])


# huts close to the peaks
set_uniques = osm_dolomites_high_peak['closest_hut'].unique()
for idx, row in osm_dolomites.loc[osm_dolomites.name.isin(set_uniques)].iterrows():
  marker = Marker(
    location = [i[0] for i in reversed(row.geometry.xy)],
    popup = utils.get_info_from_row(row), 
    tooltip = utils.get_info_from_row(row),
    icon=folium.Icon(color='blue')
  )
  marker.add_to(feature_layers['hut'])


# Linestring huts-peaks
for idx, row in linestring_distance.iterrows():
  line = folium.PolyLine(
    locations = [list(reversed(i)) for i in list(row.geometry.coords)],
    popup = utils.get_info_from_row(row), 
    tooltip = utils.get_info_from_row(row),
    color='black',
    weight=10,
  )

  line.add_to(feature_layers['peak_hut_path'])

feature_layers['dolomiti_area'].add_to(osm_peak_hut_dolomities_map)
utils.add_map_infos(feature_layers, osm_peak_hut_dolomities_map, collapsed_legend=True)

Make this Notebook Trusted to load map: File -> Trust Notebook

in the next notebook we will consider the SAT trail of the Trentino area.