Explore province of Bolzano raster data

The Geoportal of the province of Bolzano also provides raster data such as orthophotos (from various years) and Digital Terrain Models (DTM).

By means of these data we are therefore able to visualise what the real morphology of the area actually looks like. Furthermore, thanks to the DTMs we can investigate a third dimension: the altitude of this area.

Set up and data loading

import geopandas as gpd
import pandas as pd
import rasterio as rio
from owslib.wms import WebMapService
from rasterio import MemoryFile
from rasterio.plot import show
from rasterio.mask import mask
import contextily as ctx
from functions import *

province_BZ = load_province('Bolzano', 'data\Limiti01012021_g\ProvCM01012021_g')
geodf_dolomities = load_geodf_dolomities('data\I nove Sistemi delle Dolomiti UNESCO.kml')
dolomities_BZ_clipped = lead_province_clipped_dolomiti('Bolzano', 'data\Limiti01012021_g')
geo_dolomiti_BZ_mun = load_province_dolomiti_mun(21, 'data\Limiti01012021_g')
area_dolomiti_BZ = dissolve_province_dolomiti_area(21)

Explore province's geoservices

Web Map Services

In addition to offering Web Features Services, the Geoportal also has a Web Map Service that makes available data such as ortophotos and DTMs of various types. We can explore the service using owslib.wms' WebMapServices():

import warnings
warnings.filterwarnings('ignore')
wms_BZ = "http://geoservices.buergernetz.bz.it/qgis?SERVICE=WCS"
wms = WebMapService(wms_BZ)

Show all the resources:

list(wms.contents)

['P_BZ_OF_2020',
 'P_BZ_OF_2017_CIR',
 'P_BZ_OF_2014_AGEA',
 'P_BZ_OF_2014_2015_2017',
 'P_BZ_OF_2014_2015',
 'P_BZ_OF_2014',
 'P_BZ_OF_2011_CIR',
 'P_BZ_OF_2011_20cm',
 'P_BZ_OF_2011',
 'P_BZ_OF_2008',
 'P_BZ_OF_2006',
 'P_BZ_OF_2003',
 'P_BZ_OF_2000',
 'Orthophoto_1992_97',
 'Orthophoto_1982_85',
 'P_BZ_OF_1992_1997',
 'P_BZ_OF_1982_1985',
 'BASEMAP_5000_2007',
 'BASEMAP_10000_2007']

P_BZ_OF_2020 represent the 2020's ortophoto, lets' visualize it.

Get bounding box and CRS's options:

wms['P_BZ_OF_2020'].crsOptions

['CRS:84', 'EPSG:3857', 'EPSG:25832', 'EPSG:4326']

wms['P_BZ_OF_2020'].boundingBox

(1153230.0, 5807080.0, 1394440.0, 5964480.0, 'EPSG:3857')

Get the image using owslib.wms' getmap():

bbox = wms['P_BZ_OF_2020'].boundingBox[0:4]
epsg = wms['P_BZ_OF_2020'].boundingBox[4]
request = wms.getmap(
    layers=['P_BZ_OF_2020'],
    srs=epsg,
    format='image/jpeg',
    bbox=bbox,
    size=(833,606)
    )

Get the image with rasterio's MemoryFile() and visualize the ortophoto with show():

import warnings
warnings.filterwarnings('ignore')

image = MemoryFile(request).open()
show(image)

<AxesSubplot:>

The ortophoto appears to be cutted in the area of interest (Dolomiti), so let's try with older data:

wms['P_BZ_OF_2014_AGEA'].boundingBox
bbox = wms['P_BZ_OF_2014_AGEA'].boundingBox[0:4]
epsg = wms['P_BZ_OF_2020'].boundingBox[4]

request = wms.getmap(
    layers=['P_BZ_OF_2014_AGEA'],
    srs=epsg,
    format='image/jpeg',
    bbox=bbox,
    size=(833,606)
    )

raster = MemoryFile(request).open()
show(raster)

<AxesSubplot:>

These ortophoto seems to be more complete. We can therefore get an idea of the landforms in the Dolomites of the province.

Check raster dimension and colors:

raster.width

833

raster.height

606

Visualize RGB colours using rasterio's show() by specifying cmap:

show((raster, 1), cmap='Reds')

<AxesSubplot:>

show((raster, 3), cmap='Greens')

<AxesSubplot:>

show((raster, 3), cmap='Blues')

<AxesSubplot:>

Get the DTM and crop it on the province dolomiti area

We will use the DTM of the whole province of Bolzano's Alpine region to make sure that our area is included.

Open downloaded DTM with rasterio's open():

raster_dtm = rio.open('data\Raster_alpine_region\p_bz-Elevation_DigitalTerrainModel-20m-Hillshade.tif')

show(raster_dtm)

<AxesSubplot:>

Check metadata (width, height, CRS...):

raster_dtm.meta

{'driver': 'GTiff',
 'dtype': 'uint8',
 'nodata': None,
 'width': 7137,
 'height': 4305,
 'count': 1,
 'crs': CRS.from_epsg(25832),
 'transform': Affine(26.267575793, 0.0, 583578.5807889907,
        0.0, -26.267575798, 5225341.55676949)}

We will clip this raster on area_dolomiti_BZ GeoDataFrame:

area_dolomiti_BZ.plot()

<AxesSubplot:>

First we have to to get features from GeoDataFrame so that the rasterio's functions is able to process them. For this procedure we use the function getFeatures(), defined during the lactures.

def getFeatures(gdf):
    import json
    return [json.loads(gdf.to_json())['features'][0]['geometry']]

Convert the GeoDataFrame CRS and extract features:

area_dolomiti_BZ_25832 =area_dolomiti_BZ.to_crs(epsg=25832)
coords = getFeatures(area_dolomiti_BZ_25832)

Crop the raster with raterio mask() function:

raster_area_dolomiti, raster_area_dolomiti_transform = mask(raster_dtm, coords, crop=True)

Visualize Dolomiti area clipped raster:

show(raster_area_dolomiti)

<AxesSubplot:>

Save it as a GeoTiff updating metadata:

raster_area_dolomiti_meta = raster.meta

raster_area_dolomiti_meta.update({"driver": "GTiff",
                 "height": raster_area_dolomiti.shape[1],
                 "width": raster_area_dolomiti.shape[2],
                 "transform": raster_area_dolomiti_transform,
                 "count": 1})

with rio.open("data/raster_area_dolomiti_orthophoto.tif", "w", **raster_area_dolomiti_meta) as dest:
    dest.write(raster_area_dolomiti)

---

Indagate the altitude difference of the longest cableway

In the previous analysis we explored the ski slopes and ski lifts of the area. Since DTM could contain data about the 3rd dimension (altitude), let's indagate it to obtain some interesting information, such as the ski lifts' altitude difference.

In this process we will only consider one cableway, the longest one, however, as the DTM of the whole province is available, it is clearly possible to calculate the height difference of each of them. However, extending the process to all ropeways requires handling and downloading large files which we want to avoid for this notebook.

Get province's skilift data clipped on Dolomiti area:

ski_lifts_BZ_clipped = load_BZ_dolomiti_ski_lift_data("data\Ski_lifts_data")

Find the longest ski lifts/cableway:

longest = ski_lifts_BZ_clipped[ski_lifts_BZ_clipped.LENGTH == ski_lifts_BZ_clipped.LENGTH.max()]
print("The longest ski lift in Bolzano's Dolomiti's area is '", longest.NAME_IT.values[0], "' with lengyt ", round(longest.LENGTH.values[0]/10**3, 3) ," km", sep= "")

The longest ski lift in Bolzano's Dolomiti's area is 'Alpe di Siusi' with lengyt 4.086 km

longest

	OBJECTID	SKIKEY	PLR	SKIZONE	NAME_IT	NAME_DE	LENGTH	SHAPE	geometry
318	69132.0	10.02	10.0	2.0	Alpe di Siusi	Seiseralm	4086.3	None	LINESTRING (11.56396 46.54021, 11.57109 46.540...

Where is this cableway (in which municipality)?

longest_mum = geo_dolomiti_BZ_mun[geo_dolomiti_BZ_mun.contains(longest.geometry.values[0])]
longest_mum

	COD_REG	COD_PROV	PRO_COM_T	COMUNE	Shape_Area	geometry
2773	4	21	021019	Castelrotto	1.177940e+08	POLYGON ((11.56492 46.59659, 11.57231 46.59478...

We can also merge the ski lift data with municipality data using spatial join for possible further investigations: some ski lifts and cableway are contained by more than one municipality.

ski_lifts_with_mun = gpd.sjoin(ski_lifts_BZ_clipped, geo_dolomiti_BZ_mun, how='left', predicate='intersects', lsuffix='dolomities', rsuffix='provinces')

ski_lifts_with_mun[ski_lifts_with_mun.OBJECTID == longest.OBJECTID.values[0]]

	OBJECTID	SKIKEY	PLR	SKIZONE	NAME_IT	NAME_DE	LENGTH	SHAPE	geometry	index_provinces	COD_REG	COD_PROV	PRO_COM_T	COMUNE	Shape_Area
318	69132.0	10.02	10.0	2.0	Alpe di Siusi	Seiseralm	4086.3	None	LINESTRING (11.56396 46.54021, 11.57109 46.540...	2773	4	21	021019	Castelrotto	1.177940e+08

Plot the longest cableway on the dolomiti area and its municipality

base = area_dolomiti_BZ.plot(
    figsize=(15,15),
    facecolor="none",
    edgecolor="k",
    alpha=0.6,
    lw=2,)

ctx.add_basemap(base, crs=area_dolomiti_BZ.crs.to_string(), source=ctx.providers.Stamen.TerrainBackground, alpha=0.7)

geo_dolomiti_BZ_mun.plot(ax=base, 
    cmap="Pastel1", 
    column='COD_PROV',
    alpha=0.5)

dolomities_BZ_clipped.to_crs(epsg=4326).plot(
    ax=base,
    edgecolor="k",
    alpha=0.5,
    lw=0.1,
)

longest_mum.plot(ax=base, 
    color="#F65167",
    alpha=0.7
)

longest.plot(
    ax=base,
    color="red"
)

<AxesSubplot:>

Get longest features with the previous defined function:

longest_4326 = longest.to_crs(epsg=4326)
coords_longest = getFeatures(longest_4326)

Download DTM data with vertical dimension

The province's Geoportal allow to bounding the area of interest to avoid handling the whole province raster (too large a file). So we can download Castelrotto's area DTM.

Open and show Castelrotto's DTM:

raster_castelrotto = rio.open('data\Castelrotto_raster\p_bz-Elevation_DigitalTerrainModel-2.5m.tif')

show(raster_castelrotto, cmap='terrain')

<AxesSubplot:>

Check metadata:

raster_castelrotto.meta

{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': None,
 'width': 6308,
 'height': 4376,
 'count': 1,
 'crs': CRS.from_epsg(25832),
 'transform': Affine(2.4999999997440585, 0.0, 692558.7500021332,
        0.0, -2.50000000033907, 5163678.749876623)}

Extract altitude values vector:

data = raster_castelrotto.read(1)
data

array([[  849.39,   849.4 ,   849.34, ...,  2144.09,  2145.96,  2147.93],
       [  849.34,   849.19,   848.63, ...,  2144.52,  2146.23,  2148.01],
       [  849.2 ,   848.08,   847.47, ...,  2145.48,  2147.23,  2148.62],
       ...,
       [  902.13,   901.15,   900.84, ..., -9999.  , -9999.  , -9999.  ],
       [  901.4 ,   901.03,   900.94, ..., -9999.  , -9999.  , -9999.  ],
       [  902.15,   901.78,   901.4 , ..., -9999.  , -9999.  , -9999.  ]],
      dtype=float32)

Check Maximum height:

data.max()

2956.85

Check Minimum height:

data.min() # due to the shape (the blue hole in the image)

-9999.0

This value is due to the "blue hole" in the DTM image

What is the altitude difference of the longest cableway? 🚠

We need to rasform DTM CRS in order to have meter as unit:

import shapely
import pyproj
from shapely.ops import transform

wgs84 = pyproj.CRS('EPSG:4326')
crs_dtm = pyproj.CRS('EPSG:25832')
projection_transform = pyproj.Transformer.from_crs(wgs84, crs_dtm, always_xy=False).transform

We define a function to convert Points' longitude and latitude values:

def convert(x,y):
    p = shapely.geometry.Point(y,x)
    p = transform(projection_transform,p)
    return(p)

Create two arrays for x and y (longitude and latitude) applying the previous defined function:

pointsx = []
pointsy = []
for coordinate in longest.geometry.values[0].coords:
  x = coordinate[0]
  y = coordinate[1]
  point = convert(x,y)
  pointsx.append(point.x)
  pointsy.append(point.y)

Calculate the distance point-to-point (we basically obtain the length of the linestring updated at each point):

from shapely.geometry import Point, LineString
lengths = []
previousPoint = None
length = 0
for i in range(len(pointsx)):
  point = shapely.geometry.Point(pointsy[i],pointsx[i])
  if previousPoint is None:
    lengths.append(length)
  else:
    length = LineString([previousPoint,point]).length + length
    lengths.append(length)  
  previousPoint = point 

# list of updated lenghts point to point:
lengths

[0,
 548.1321324826865,
 671.6382309110555,
 739.7570003143445,
 845.8157132330692,
 1368.45268512941,
 2003.9914533135588,
 2774.6864939567545,
 2827.774460389724,
 4007.983852219015,
 4021.02225702942,
 4086.297505033714]

Extract the altitude value for each point:

rows, cols = rio.transform.rowcol(raster_castelrotto.transform,(pointsx),(pointsy))

values = []
for i in range(len(rows)):
    values.append(data[rows[i]-1][cols[i]-1]) 

values

[1009.19,
 1062.27,
 1082.27,
 1076.65,
 1097.02,
 1214.6,
 1575.68,
 1693.55,
 1696.76,
 1854.81,
 1854.47,
 1854.85]

Save altitude values for each point in a DataFrame:

siusi_3d = pd.DataFrame()
siusi_3d['value'] = values
siusi_3d['length'] = lengths
siusi_3d

	value	length
0	1009.190002	0.000000
1	1062.270020	548.132132
2	1082.270020	671.638231
3	1076.650024	739.757000
4	1097.020020	845.815713
5	1214.599976	1368.452685
6	1575.680054	2003.991453
7	1693.550049	2774.686494
8	1696.760010	2827.774460
9	1854.810059	4007.983852
10	1854.469971	4021.022257
11	1854.849976	4086.297505

Plot the data:

ax = siusi_3d.plot(y="value",x="length",color='green',figsize=(28,8))
ax.plot()

[]

As mentioned on the Dolomities Seise Alm website, the cableway starts from 1000 metres and goes up to 1800 metres, thus our results seems to be quite accurate:

"L’area escursionistica si raggiunge in tutta comodità e rapidità con moderni impianti di risalita collegati in maniera impeccabile ai sentieri dell’altipiano. I bus-navetta conducono dalle diverse località di Castelrotto, Siusi allo Sciliar, Fiè allo Sciliar e Tires al Catinaccio alla stazione a valle della cabinovia Alpe di Siusi, che consente di passare velocemente da 1.000 a 1.800 m s.l.m."