<jupyter_start><jupyter_text>Temperature forecast for the general public (MAELSTROM-Yr dataset)This dataset contains temperature weather forecast for the Nordic region, and are used to produce public weather forecasts on the weather app Yr (www.yr.no). The goal of the prediction task is to generate a deterministic temperature forecast together with an uncertainty range (10% to 90%) as shown here: https://www.yr.no/en/details/graph/5-18700/Norway/Oslo/Oslo/Oslo%20(Blindern).The target field in the dataset is constructed using a high density network of citizen weather stations from [Netatmo](https://weathermap.netatmo.com/).The current operational implementation uses a very simple regression model based on only a subset of the predictors available in the dataset. It is described in this article: https://journals.ametsoc.org/view/journals/bams/101/1/bams-d-18-0237.1.xml PrerequisitesTo run the code in this notebook, you need the following packages:`pip install climetlab climetlab_maelstrom_yr keras tensorflow numpy matplotlib`  Loading the dataWe can use climetlab to load the dataset into an xarray dataset. There will be several datasets available of different sizes: 300 MB (not available yet), 5GB, and 5TB (not available yet). The 5TB dataset contains the entire Nordic domain at 1x1 km resolution for all 60 hour leadtimes. The 5GB dataset contains only a subset of grid points (128x128) surrounding Oslo, Norway and only for leadtimes 6, 12, ..., 42 hours. All datasets contain the same input predictors and time period (4 years).Currently, only "air_temperature" is available as the predictand parameter, however precipitation_amount will be added in the future.The entire 5GB dataset will take a few minutes to load, since the data must be downloaded from europeanweather.cloud. Climetlab caches files locally, so files need not be when rerunning the code later. To only load a subset, add a dates argument to load_dataset, e.g. `dates=['2017-01-01', '2017-01-02']` or `dates=pandas.date_range(start="2017-01-01", end="2017-03-01", freq="1D")`.<jupyter_code>import climetlab as cml
import pandas
cmlds = cml.load_dataset(
        'maelstrom-yr',
        size='5GB',
        parameter='air_temperature',
    )
ds = cmlds.to_xarray()<jupyter_output>By downloading data from this dataset, you agree to the terms and conditions defined at https://github.com/metno/maelstrom_yr/LICENSE. If you do not agree with such terms, do not download the data.<jupyter_text>This dataset contains the following dimensions and variables<jupyter_code>print(ds)<jupyter_output><xarray.Dataset>
Dimensions:          (record: 1457, x: 128, y: 128, leadtime: 8, predictor: 49, char: 50)
Coordinates:
  * x                (x) float32 -2.974e+05 -2.964e+05 ... -1.714e+05 -1.704e+05
  * y                (y) float32 -5.043e+05 -5.033e+05 ... -3.783e+05 -3.773e+05
  * leadtime         (leadtime) timedelta64[ns] 00:00:00 ... 1 days 18:00:00
    latitude         (y, x) float32 dask.array<chunksize=(128, 128), meta=np.ndarray>
    longitude        (y, x) float32 dask.array<chunksize=(128, 128), meta=np.ndarray>
Dimensions without coordinates: record, predictor, char
Data variables:
    time             (record) datetime64[ns] 2017-01-01 ... 2020-12-31
    projection       int32 -2147483647
    altitude         (y, x) float32 dask.array<chunksize=(128, 128), meta=np.ndarray>
    predictors       (record, y, x, predictor) float32 dask.array<chunksize=(1, 128, 128, 49), meta=np.ndarray>
    name_predictor   (predictor, char) |S1 dask.array<chunksize=(49, 50), meta=np.ndarray[...]<jupyter_text>The dataset is mostly self explanatory. The `record` dimension represent different samples. The `predictor` variable contains all predictors stacked one after the other, including values for different leadtimes. The `target` variable contain target values.  Plotting predictors and predictand (target)<jupyter_code>import matplotlib.pyplot as plt
import numpy as np

names = ds["name_predictor"].values
names = np.array([''.join([qq.decode('utf-8') for qq in names[p, :]]) for p in range(names.shape[0])])

num_leadtimes = len(ds["leadtime"])

unique_predictor_names = np.unique(names)
print("Available predictors:", unique_predictor_names)

index_date = 0
target = ds["target"].values
plt.rcParams['figure.dpi'] = 100
plt.rcParams['figure.figsize'] = [10, 6]
for i, name in enumerate(unique_predictor_names):
    plt.subplot(2, 4, i + 1)
    index = np.where(names == name)[0][0]
    plt.pcolormesh(ds["predictors"][index_date, :, :, index], shading="auto", rasterized=True)
    plt.gca().set_aspect(1)
    plt.title(name)
plt.subplot(2, 4, 8)
plt.pcolormesh(target[index_date, :, :, 0], shading="auto", rasterized=True)
plt.gca().set_aspect(1)
plt.title("Target")<jupyter_output><empty_output><jupyter_text>Example ML solution Normalizing the predictorsFirst we normalize the predictors, by subtractng the mean and dividing by the standard deviation:<jupyter_code>raw_forecast = np.copy(ds["predictors"][:, :, :, 0:num_leadtimes])
predictors = np.copy(ds["predictors"].values)
num_predictors = predictors.shape[3]
for p in range(num_predictors):
    predictors[:, :, :, p] -= np.nanmean(predictors[:, :, :, p])
    predictors[:, :, :, p] /= np.nanstd(predictors[:, :, :, p])<jupyter_output><empty_output><jupyter_text>Defining the loss functionWe use the quantile loss function, by scoring each of the three output quantiles of the model:<jupyter_code>import keras
import tensorflow as tf
import keras.backend as K

global num_leadtimes
def quantile_loss_function(y_true, y_pred):
    err0 = y_true - y_pred[:, :, :, 0:num_leadtimes]
    err1 = y_true - y_pred[:, :, :, num_leadtimes:(2*num_leadtimes)]
    err2 = y_true - y_pred[:, :, :, (2*num_leadtimes):(3*num_leadtimes)]
    qtloss0 = (0.5 - tf.cast((err0 < 0), tf.float32)) * err0
    qtloss1 = (0.1 - tf.cast((err1 < 0), tf.float32)) * err1
    qtloss2 = (0.9 - tf.cast((err2 < 0), tf.float32)) * err2
    return K.mean(qtloss0 + qtloss1 + qtloss2)<jupyter_output>2021-08-13 16:55:46.584476: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory
2021-08-13 16:55:46.584505: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.<jupyter_text>Setting up the modelThe model takes a gridded predictor set as input and outputs gridded fields for each leadtime and for three quantiles. The tempearture forecast on yr.no has both a deterministic best guess and a 10-90% confidence interval. We want the model to predict all three parameters simultaneously.<jupyter_code>num_quantiles = 3
num_outputs = num_quantiles * num_leadtimes

model = keras.Sequential()
model.add(keras.layers.InputLayer(predictors.shape[1:]))
model.add(keras.layers.Dense(num_outputs))

model.compile(optimizer = tf.keras.optimizers.Adam(learning_rate = 1e-2), loss = quantile_loss_function)
model.summary()<jupyter_output>Model: "sequential_1"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
dense_1 (Dense)              (None, 128, 128, 24)      1200      
=================================================================
Total params: 1,200
Trainable params: 1,200
Non-trainable params: 0
_________________________________________________________________<jupyter_text>Training the modelWe will split the dataset into a training and evaluation set, based on the record dimension.<jupyter_code>Itrain = range(predictors.shape[0]//2)
Ieval = range(predictors.shape[0]//2, predictors.shape[0])
num_epochs = 50
batch_size = 4
model.fit(predictors[Itrain, ...], target[Itrain, ...], epochs=num_epochs, batch_size=batch_size)<jupyter_output>2021-08-13 17:00:00.819926: W tensorflow/core/framework/cpu_allocator_impl.cc:80] Allocation of 2337800192 exceeds 10% of free system memory.
2021-08-13 17:00:02.380973: W tensorflow/core/framework/cpu_allocator_impl.cc:80] Allocation of 381681664 exceeds 10% of free system memory.
2021-08-13 17:00:02.716530: I tensorflow/compiler/mlir/mlir_graph_optimization_pass.cc:185] None of the MLIR Optimization Passes are enabled (registered 2)<jupyter_text>Predict output<jupyter_code>output = model.predict(predictors[Ieval, ...])<jupyter_output>2021-08-13 17:04:26.565755: W tensorflow/core/framework/cpu_allocator_impl.cc:80] Allocation of 2341011456 exceeds 10% of free system memory.
2021-08-13 17:04:28.901367: W tensorflow/core/framework/cpu_allocator_impl.cc:80] Allocation of 1146617856 exceeds 10% of free system memory.<jupyter_text>Model evaluation and visualization  Evaluating the modelFirst, lets compare the mean absolute error of the raw forecast and the ML-forecast of the median<jupyter_code>import numpy as np
print("Raw model MAE:", np.nanmean(np.abs(raw_forecast[Ieval, ...] - target[Ieval, ...])), "°C")
print("ML MAE:", np.nanmean(np.abs(output[:, :, :, 0:num_leadtimes] - target[Ieval, ...])), "°C")<jupyter_output>Raw model MAE: 1.1301652 °C
ML MAE: 0.8665128 °C<jupyter_text>Next, we can plot the MAE as a function of leadtime:<jupyter_code>x = ds["leadtime"].astype(float) / 3600 / 1e9
plt.plot(x, [np.nanmean(np.abs(output[:, :, :, i] - target[Ieval, :, :, i])) for i in range(num_leadtimes)], 'ro-', lw=2, label="Model")
plt.plot(x, [np.nanmean(np.abs(raw_forecast[Ieval, :, :, i] - target[Ieval, :, :, i])) for i in range(num_leadtimes)], 'yo-', lw=2, label="Raw")
plt.legend()
plt.xlabel("Lead time (hours)")
plt.ylabel("Mean absolute error (°C)")<jupyter_output><empty_output><jupyter_text>Visualizing the results as timeseriesWe can visualize the output as a timeseries. We will pick an example point (Oslo).<jupyter_code>Y = 55
X = 55
plt.plot(x, output[0, Y, X, 0:num_leadtimes], 'r-', lw=2, label="Median")
plt.plot(x, raw_forecast[Ieval[0], Y, X, 0:num_leadtimes], 'y-', lw=2, label="Raw")
lower = output[0, Y, X,num_leadtimes:2*num_leadtimes]
upper = output[0, Y, X, 2*num_leadtimes:3*num_leadtimes]
plt.plot(x, lower, 'r--', lw=2, label="10%")
plt.plot(x, upper, 'r--', lw=2, label="90%")
xx = np.concatenate((x, x[::-1]))
plt.fill(np.concatenate((x, x[::-1])), np.concatenate((lower, upper[::-1])), color='r', alpha=0.2, linewidth=0)
plt.plot(x, target[Ieval[0], Y, X, :], 'bo-', lw=2, label="Target")
plt.legend()
plt.xlabel("Lead time (hours)")
plt.ylabel("Air temperature (°C)")<jupyter_output><empty_output><jupyter_text>Visualizing the results on a map<jupyter_code>plt.subplot(1, 3, 1)
plt.pcolormesh(raw_forecast[Ieval[0], :, :, 0], rasterized=True)
plt.gca().set_aspect(1)
plt.title("Raw forecast")

plt.subplot(1, 3, 2)
plt.pcolormesh(output[0, :, :, 0], rasterized=True)
plt.gca().set_aspect(1)
plt.title("ML forecast")

plt.subplot(1, 3, 3)
plt.pcolormesh(target[Ieval[0], :, :, 0], rasterized=True)
plt.gca().set_aspect(1)
plt.title("Target (median)")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1413465883416177, 'fraction_numerical': 0.05503840939900587, 'mean_word_length': 5.199439775910364, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 47, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '30893', 'license': 'Apache-2.0', 'n_tokens_mistral': 4204, 'n_tokens_neox': 3732, 'n_words': 1037, 'path': 'notebooks/demo_yr.ipynb', 'repo_name': 'metno/maelstrom-yr', 'size': '11065'}