Analyse des vergangenen Sommers¶

Dieses Notebook kann lokal oder direkt im Browser auf oder ausgeführt werden.

In [119]:

import requests 
import pandas as pd
import geopandas as gpd 
import matplotlib.pyplot as plt # for plotting
import folium
import branca
import math

Einschlaftemperatur Karte¶

In [ ]:

# Auswahl an Zeitraum (z.B. eine Hitzewelle, ein Monat, oder der ganze Sommer)
# Daten (UTC) als String in folgendem Format: "YYYY-MM-DDThh:mm:ssZ"
time_from = "2024-10-01T22:00:00Z"
time_to = "2024-10-31T22:00:00Z"

# Auswahl an Stationen
# Wenn alle Stationen berücksichtigt werden sollen, dann einfach die Liste leer lassen
# Hier ist eine Liste aller Stationen mit Name und ID: https://smart-urban-heat-map.ch/api/v2/latest
# station_ids = ["11001", "11002"] # Beispiel für nur eine Auswahl von zwei Stationen
station_ids = [] # Beispiel für alle Stationen

In [ ]:

def get_stations() -> gpd.GeoDataFrame:
    response = requests.get(
            url=f"https://smart-urban-heat-map.ch/api/v2/latest"
    )
    stations = gpd.read_file(response.text)

    # calculate only at subset of stations 
    # (if a subset was defined, otherwise take all stations)
    if station_ids:
        stations = stations[stations["stationId"].isin(station_ids)]
        
    return stations

def get_station_sleep_temp_analysis(time_from: str, time_to: str) -> pd.DataFrame:
    stations = get_stations()

    # check that time series are at least 83 % complete
    time_from_dt = pd.to_datetime(time_from)
    time_to_dt = pd.to_datetime(time_to)
    time_difference = time_to_dt - time_from_dt
    expected_values = round(time_difference.days * 24 * 2)
    
    stations['einschlaftemperatur'] = None

    for idx, station in stations.iterrows():
        station_id = station.stationId
        response = requests.get(url=f"https://smart-urban-heat-map.ch/api/v2/timeseries?stationId={station_id}&timeFrom={time_from}&timeTo={time_to}")
        
        payload = response.json()

        if payload["values"] is None or not len(payload["values"]): 
            stations = stations.drop(idx)
            continue
            
        if len(payload["values"]) < expected_values:
            stations = stations.drop(idx)
            continue
            
        df = pd.DataFrame(payload["values"])

        df["dateObserved"] = pd.to_datetime(df["dateObserved"])
        df["dateObserved"] = df["dateObserved"].dt.tz_convert("Europe/Zurich")
        
        # hier werden die Hitzetage und Tropennächte berechnet
        einschlaftemperatur = calc_einschlaftemperatur(df)
        
        stations.loc[idx, 'einschlaftemperatur'] = einschlaftemperatur

    return stations

def calc_einschlaftemperatur(df: pd.DataFrame) -> float:
    sleep_time_mask = (df['dateObserved'].dt.hour >= 22) & (df['dateObserved'].dt.hour < 23)
    sleep_time = df.loc[sleep_time_mask]

    nightly_sleep_temperatures = (
        sleep_time.groupby(sleep_time['dateObserved'].dt.date)['temperature'].mean()
    )

    return nightly_sleep_temperatures.mean()

In [122]:

# run 'einschlaftemperatur' analysis
# this will take a few minutes depending on time and nr. of stations
station_analysis = get_station_sleep_temp_analysis(time_from, time_to)
station_analysis

Out[122]:

	name	stationId	geometry	einschlaftemperatur
0	Ausserholligen 2, ewb	11001	POINT (7.40642 46.94542)	21.239814
1	Bundesplatz	11002	POINT (7.44353 46.94692)	22.099029
2	Breitenrain, Waffenweg	11003	POINT (7.45192 46.96173)	21.727543
3	Schosshaldenfriedhof 2	11004	POINT (7.47186 46.95339)	20.17552
4	Monbijou-Park	11005	POINT (7.43462 46.94187)	20.943326
...	...	...	...	...
134	Monopoliplatz Lyss	12006	POINT (7.30588 47.07661)	21.617444
135	Spielplatz Stiglimatt	12007	POINT (7.2992 47.07154)	20.998686
136	Sportanlage Grien	12008	POINT (7.29574 47.07508)	20.725683
138	Reitplatz Grünau	12010	POINT (7.30167 47.07714)	20.846787
139	Sigriswil	13001	POINT (7.71244 46.71345)	20.11358

113 rows × 4 columns

In [123]:

# Einschlaftemperatur Map
m = folium.Map(location=[station_analysis.geometry.y.mean(), station_analysis.geometry.x.mean()], zoom_start=13, tiles="CartoDB positron")

# Add a fixed title to the map
title_html = f'''
     <div style="position: fixed; 
     top: 20px; left: 100px; width: 25%; height: 45px; 
     background-color: #F0F0F0; border: 1px solid black; z-index: 9999; font-size: 14px; font-weight: bold;">
     Einschlaftemperatur zwischen 22:00 und 23:00<br> ({pd.to_datetime(time_from).strftime('%d.%m.%Y')} bis {pd.to_datetime(time_to).strftime('%d.%m.%Y')})
     </div>
     '''
m.get_root().html.add_child(folium.Element(title_html))

colormap = branca.colormap.linear.YlOrRd_09
einschlaftemperatur = station_analysis.einschlaftemperatur.values

# Define colourmap range depending on values of 'einschlaftemperatur
vmin = math.floor(einschlaftemperatur.min())
vmax = math.ceil(einschlaftemperatur.max())

# Define the colormap with the specified range
colormap = branca.colormap.linear.YlOrRd_09.scale(vmin, vmax)

# Convert to step colormap with a specified number of steps
n_steps = int((vmax - vmin) / 0.5)  # Define the number of steps
colormap = colormap.to_step(n_steps)

# colormap = colormap.scale(0, einschlaftemperatur).to_step(einschlaftemperatur) 
colormap.caption = "Mittlere Einschlaftemperatur"
colormap.add_to(m)

# plot each station temperature
for idx, station in station_analysis.iterrows():
    color = colormap(station.einschlaftemperatur)
    # text with temperature value
    folium.Marker(
        location=(station.geometry.y, station.geometry.x),
        icon=folium.DivIcon(
            html=f'<div style="font-size: 10pt; color: {color}; text-shadow: -1px -1px 0 #D3D3D3, 1px -1px 0 #D3D3D3, -1px 1px 0 #D3D3D3, 1px 1px 0 #D3D3D3;">{station.einschlaftemperatur:.1f}°C</div>'
            ),
        tooltip=f"{station['name']}: Einschlaftemperatur {station.einschlaftemperatur:.1f} °C",
    ).add_to(m)

# show map
m

Out[123]:

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

"Urban Heat Island" Effekt Grafik¶

In [ ]:

# Auswahl an Zeitraum
# get data of specific period (i.e. a heatwave)
uhi_time_from = "2024-10-01T22:00:00Z" # startzeit
uhi_time_to = "2024-10-31T23:50:00Z" # endzeit

# Auswahl an Stationen
# Wenn alle Stationen berücksichtigt werden sollen, dann einfach die Liste leer lassen
station_ids = [
                    # # stationen in der stadt
                    "11002", # Bundesplatz
                    "11019", # Breitenrainplatz 127
                    # stationen in quartieren der agglo oder ausserhalb
                    "11090", # Köniz Blinzernplateau Strom-/Telefonmast
                    "11010", # Bümpliz Stöckacker
                    "11098",  # Ostermundigen Schermenweg,
                    "11100", # Umland Bolligen
                    "11121", # Zollikofen Zentrum Laterne,
                    # "11075", # Gurten Kulm, 
                     ] # Beispiel für nur eine Auswahl von zwei Stationen

# hier ist eine Liste alle Stationen mit Name und ID:
# https://smart-urban-heat-map.ch/api/v2/latest

In [ ]:

# code for uhi calculation

# define zollikofen as rural reference station
rural_station =  {
    "name": "Zollikofen 2m",
    "stationId": "11134",
      }

# uhi analysis
def get_station_uhi_analysis(time_from: str, time_to: str) -> pd.DataFrame:
    stations = get_stations()

    # check that time series are at least 83 % complete
    time_from_dt = pd.to_datetime(time_from)
    time_to_dt = pd.to_datetime(time_to)
    time_difference = time_to_dt - time_from_dt
    expected_values = round(time_difference.days * 24 * 2) # minimum two values per hour
    
    stations['uhi'] = None
    # get values of rural reference station
    rural_station_id = rural_station['stationId']
    response = requests.get(url=f"https://smart-urban-heat-map.ch/api/v2/timeseries?stationId={rural_station_id}&timeFrom={time_from}&timeTo={time_to}")

    payload = response.json()
    
    if payload["values"] is None or not len(payload["values"]): 
        raise ValueError(f"No data for the rural reference station {rural_station['name']}. Aborting UHI analysis.")
    
    if len(payload["values"]) < expected_values:
        raise ValueError(f"Insufficient data for the rural reference station {rural_station['name']}. Aborting UHI analysis.")

    df_rural_station = pd.DataFrame(payload["values"])
    df_rural_station["dateObserved"] = pd.to_datetime(df_rural_station["dateObserved"])
    df_rural_station["dateObserved"] = df_rural_station["dateObserved"].dt.tz_convert("Europe/Zurich")
    df_rural_station_mean_hourly_temp = df_rural_station.groupby(df_rural_station['dateObserved'].dt.hour)['temperature'].mean()

    for idx, station in stations.iterrows():
        station_id = station.stationId
        response = requests.get(url=f"https://smart-urban-heat-map.ch/api/v2/timeseries?stationId={station_id}&timeFrom={time_from}&timeTo={time_to}")
        
        payload = response.json()

        if payload["values"] is None or not len(payload["values"]): 
            stations = stations.drop(idx)
            continue
            
        if len(payload["values"]) < expected_values:
            stations = stations.drop(idx)
            continue
            
        df = pd.DataFrame(payload["values"])
        df["dateObserved"] = pd.to_datetime(df["dateObserved"])
        df["dateObserved"] = df["dateObserved"].dt.tz_convert("Europe/Zurich")
        df_urban_station_mean_hourly_temp = df.groupby(df['dateObserved'].dt.hour)['temperature'].mean()
        
        # calculate uhi
        uhi = calc_uhi(df_urban_station_mean_hourly_temp, df_rural_station_mean_hourly_temp)
        uhi.index.name = 'uhi'
        
        # write uhi to stations gdf
        stations.at[idx, 'uhi'] = uhi

    return stations.dropna(subset=['uhi'])  # Return only the stations with UHI data

# UHI function
def calc_uhi(temperatures_urban_station: pd.Series, temperatures_rural_station: pd.Series) -> pd.Series:
    """calculate temperature difference between two stations"""
    # calculate UHI
    uhi = temperatures_urban_station - temperatures_rural_station
    return uhi

In [126]:

# run uhi analysis
station_analysis = get_station_uhi_analysis(uhi_time_from, uhi_time_to)
station_analysis

Out[126]:

	name	stationId	geometry	uhi
1	Bundesplatz	11002	POINT (7.44353 46.94692)	uhi 0 1.372617 1 1.285352 2 1.2579...
9	Bümpliz Stöckacker	11010	POINT (7.39587 46.94302)	uhi 0 0.343092 1 0.389947 2 0.4306...
17	Breitenrainplatz 127	11019	POINT (7.45403 46.95825)	uhi 0 1.017292 1 1.022750 2 1.0620...
87	Köniz Blinzernplateau Strom-/Telefonmast	11090	POINT (7.42771 46.92253)	uhi 0 1.159601 1 1.166302 2 1.2812...
95	Ostermundigen Schermenweg	11098	POINT (7.48036 46.95778)	uhi 0 0.447994 1 0.485579 2 0.5590...
97	Umland Bolligen	11100	POINT (7.50376 46.96681)	uhi 0 -1.181658 1 -1.239773 2 -1.2348...
118	Zollikofen Zentrum Laterne	11121	POINT (7.45693 46.99266)	uhi 0 1.404620 1 1.419536 2 1.4175...

In [127]:

import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd

# Function to reorder the values in the series
def reorder_series(series):
    first_part = series.iloc[12:24].values
    second_part = series.iloc[0:12].values
    reordered_series = list(first_part) + list(second_part)
    return reordered_series

# Create a new figure and axis for the plot
fig, ax = plt.subplots(figsize=(12, 8))

# Plot each series in the 'uhi' column
for idx, row in station_analysis.iterrows():
    uhi_series = row['uhi']
    reordered_uhi = reorder_series(uhi_series)
    x_values = list(range(24))  # 0 to 23
    ax.plot(x_values, reordered_uhi, label=row['name'])

# Customize the plot
ax.set_xlabel('Tageszeit [CEST]')
ax.set_ylabel('°C')
fig.suptitle('Städtischer Wärmeinsel Effekt', fontsize=14)
ax.set_title('Temperatur an den Stationen minus Temperatur an der Station "Zollikofen 3m"', fontsize=10)
ax.legend(loc='lower right', fontsize="6")

# Adjust x-axis labels
ax.set_xticks(x_values)
ax.set_xticklabels([str((i+12)%24) for i in range(24)])

# Add solid black line at y=0
ax.axhline(y=0, color='black', linewidth=1)

# Add fine black dashed lines at every 0.5 interval on the y-axis
y_min, y_max = ax.get_ylim()
y_ticks = [y for y in range(int(y_min*2), int(y_max*2)+1)]
for y in y_ticks:
    if y != 0:  # Skip y=0 as it's already added as a solid line
        ax.axhline(y=y/2.0, color='black', linestyle='--', linewidth=0.5)

# Add slightly transparent fine grey solid lines at each x interval
for x in x_values:
    ax.axvline(x=x, color='grey', linestyle='-', linewidth=0.5, alpha=0.5)

# Add a text box in the top right corner
textstr = f"von {pd.to_datetime(uhi_time_from).strftime('%d.%m.%Y')} bis {pd.to_datetime(uhi_time_to).strftime('%d.%m.%Y')}"
props = dict(boxstyle='round', facecolor='white', alpha=0.5)
ax.text(0.95, 0.95, textstr, transform=ax.transAxes, fontsize=10,
        verticalalignment='top', horizontalalignment='right', bbox=props)

# Show the plot
plt.show()