Preserve zero and positive values and global maxima with safeguards:¶
Example on precipitation from observations and reanalysis¶
This example explores the effects of applying three different lossy compressors (ZFP, SZ3, SPERR) on a time series (3 days, hourly intervals) of precipitation PR. The time series come from observations at Belém, Brazil and Helsinki, Finland, and the corresponding closest grid points in the reanalysis product ERA5. The compressors are applied to either the time series of the individual observations or the global grid of ERA5 (before extracting the observation-space values from the closest grid points). Finally, we apply safeguards to guarantee an absolute error bound and that zero values, positive values, and global extrema are preserved.
The meteorological relevant properties of precipitation we study are the time and magnitude of the maximum, the integral over the time series, the occurrences of (positive) precipitation (including false positives and negatives), and the occurrences of negative values after compression.
Copied!
import colorsys
from pathlib import Path
import humanize
import matplotlib.colors as mc
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
import colorsys
from pathlib import Path
import humanize
import matplotlib.colors as mc
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
Copied!
# Retrieve the data
Belem = pd.read_csv(Path() / "data" / "obs-pr" / "belem.csv")
Helsinki = pd.read_csv(Path() / "data" / "obs-pr" / "helsinki.csv")
ERA5 = xr.open_dataset(Path() / "data" / "era5-pr" / "data.nc")
# Retrieve the data
Belem = pd.read_csv(Path() / "data" / "obs-pr" / "belem.csv")
Helsinki = pd.read_csv(Path() / "data" / "obs-pr" / "helsinki.csv")
ERA5 = xr.open_dataset(Path() / "data" / "era5-pr" / "data.nc")
Copied!
# Extract the data variables
Time = ERA5["valid_time"].data
ERA5_PR = ERA5["tp"] * 1000
Belem_PR = Belem["PR"]
Helsinki_PR = Helsinki["PR"]
# Extract the data variables
Time = ERA5["valid_time"].data
ERA5_PR = ERA5["tp"] * 1000
Belem_PR = Belem["PR"]
Helsinki_PR = Helsinki["PR"]
Copied!
ERA5_Belem_PR = ERA5_PR.sel(
latitude=-1.4563, longitude=360 - 48.5013, method="nearest"
).data
ERA5_Helsinki_PR = ERA5_PR.sel(
latitude=60.1699, longitude=24.9384, method="nearest"
).data
ERA5_Belem_PR = ERA5_PR.sel(
latitude=-1.4563, longitude=360 - 48.5013, method="nearest"
).data
ERA5_Helsinki_PR = ERA5_PR.sel(
latitude=60.1699, longitude=24.9384, method="nearest"
).data
Copied!
ERA5_PR.shape, Belem_PR.shape, Helsinki_PR.shape
ERA5_PR.shape, Belem_PR.shape, Helsinki_PR.shape
((72, 721, 1440), (72,), (72,))
Copied!
def plot_positive_precipitation(pr, ax, y, c, lw=2, bars=False):
foox = []
fooy = []
lastx = None
wasgap = pr[0] <= 0.0
for t, p in zip(Time, pr):
if p > 0.0:
if wasgap:
wasgap = False
foox.append(lastx)
else:
foox.append(t)
fooy.append(y)
else:
if not wasgap:
wasgap = True
foox.append(lastx)
fooy.append(y)
ax.plot(foox, fooy, lw=lw, c=c, solid_capstyle="butt")
if bars:
ax.plot(
[foox[0], foox[0]],
[fooy[0] - 0.75, fooy[0] + 0.75],
lw=lw,
c=c,
solid_capstyle="butt",
)
ax.plot(
[foox[-1], foox[-1]],
[fooy[-1] - 0.75, fooy[-1] + 0.75],
lw=lw,
c=c,
solid_capstyle="butt",
)
foox = []
fooy = []
lastx = t
if len(foox) > 0:
foox.append(t + (t - lastx))
fooy.append(y)
ax.plot(foox, fooy, lw=lw, c=c, solid_capstyle="butt")
if bars:
ax.plot(
[foox[0], foox[0]],
[fooy[0] - 0.75, fooy[0] + 0.75],
lw=lw,
c=c,
solid_capstyle="butt",
)
ax.plot(
[foox[-1], foox[-1]],
[fooy[-1] - 0.75, fooy[-1] + 0.75],
lw=lw,
c=c,
solid_capstyle="butt",
)
def plot_positive_precipitation(pr, ax, y, c, lw=2, bars=False):
foox = []
fooy = []
lastx = None
wasgap = pr[0] <= 0.0
for t, p in zip(Time, pr):
if p > 0.0:
if wasgap:
wasgap = False
foox.append(lastx)
else:
foox.append(t)
fooy.append(y)
else:
if not wasgap:
wasgap = True
foox.append(lastx)
fooy.append(y)
ax.plot(foox, fooy, lw=lw, c=c, solid_capstyle="butt")
if bars:
ax.plot(
[foox[0], foox[0]],
[fooy[0] - 0.75, fooy[0] + 0.75],
lw=lw,
c=c,
solid_capstyle="butt",
)
ax.plot(
[foox[-1], foox[-1]],
[fooy[-1] - 0.75, fooy[-1] + 0.75],
lw=lw,
c=c,
solid_capstyle="butt",
)
foox = []
fooy = []
lastx = t
if len(foox) > 0:
foox.append(t + (t - lastx))
fooy.append(y)
ax.plot(foox, fooy, lw=lw, c=c, solid_capstyle="butt")
if bars:
ax.plot(
[foox[0], foox[0]],
[fooy[0] - 0.75, fooy[0] + 0.75],
lw=lw,
c=c,
solid_capstyle="butt",
)
ax.plot(
[foox[-1], foox[-1]],
[fooy[-1] - 0.75, fooy[-1] + 0.75],
lw=lw,
c=c,
solid_capstyle="butt",
)
Copied!
def plot_negative_precipitation(pr, ax, y):
foox = []
fooy = []
lastx = None
wasgap = pr[0] >= 0.0
for t, p in zip(Time, pr):
if p < 0.0:
if wasgap:
wasgap = False
foox.append(lastx)
else:
foox.append(t)
fooy.append(y)
else:
if not wasgap:
wasgap = True
foox.append(lastx)
fooy.append(y)
ax.plot(foox, fooy, lw=4, c="red", solid_capstyle="butt")
foox = []
fooy = []
lastx = t
if len(foox) > 0:
foox.append(t + (t - lastx))
fooy.append(y)
ax.plot(foox, fooy, lw=4, c="red", solid_capstyle="butt")
def plot_negative_precipitation(pr, ax, y):
foox = []
fooy = []
lastx = None
wasgap = pr[0] >= 0.0
for t, p in zip(Time, pr):
if p < 0.0:
if wasgap:
wasgap = False
foox.append(lastx)
else:
foox.append(t)
fooy.append(y)
else:
if not wasgap:
wasgap = True
foox.append(lastx)
fooy.append(y)
ax.plot(foox, fooy, lw=4, c="red", solid_capstyle="butt")
foox = []
fooy = []
lastx = t
if len(foox) > 0:
foox.append(t + (t - lastx))
fooy.append(y)
ax.plot(foox, fooy, lw=4, c="red", solid_capstyle="butt")
Copied!
# based on https://stackoverflow.com/a/49601444
def adjust_color_lightness(color):
try:
c = mc.cnames[color]
except Exception:
c = color
c = colorsys.rgb_to_hls(*mc.to_rgb(c))
if c[1] < 0.5:
return colorsys.hls_to_rgb(c[0], 0.66, c[2])
else:
return colorsys.hls_to_rgb(c[0], 0.33, c[2])
# based on https://stackoverflow.com/a/49601444
def adjust_color_lightness(color):
try:
c = mc.cnames[color]
except Exception:
c = color
c = colorsys.rgb_to_hls(*mc.to_rgb(c))
if c[1] < 0.5:
return colorsys.hls_to_rgb(c[0], 0.66, c[2])
else:
return colorsys.hls_to_rgb(c[0], 0.33, c[2])
Copied!
def compute_corrections_percentage(my_PR, orig_PR) -> float:
return np.mean(my_PR != orig_PR)
def compute_corrections_percentage(my_PR, orig_PR) -> float:
return np.mean(my_PR != orig_PR)
Copied!
def plot_precipitation(
ax1,
ax2,
my_ERA5_PR: xr.DataArray,
my_ERA5_PR_cr: float,
my_Belem_PR: pd.DataFrame,
my_Belem_PR_cr: float,
my_Helsinki_PR: pd.DataFrame,
my_Helsinki_PR_cr: float,
PR_eb_abs: float,
title: str,
reference: bool = False,
corr=None,
):
my_ERA5_Belem_PR = my_ERA5_PR.sel(
latitude=-1.4563, longitude=360 - 48.5013, method="nearest"
).data
my_ERA5_Helsinki_PR = my_ERA5_PR.sel(
latitude=60.1699, longitude=24.9384, method="nearest"
).data
pos = mdates.HourLocator(byhour=(0, 6, 12, 18))
fmt = mdates.DateFormatter("%d.%m %Hh")
if reference:
ax1.fill_between(Time, my_Belem_PR, alpha=0.5, step="pre")
ax1.fill_between(Time, my_Helsinki_PR, alpha=0.5, step="pre")
bp = ax1.plot(Time, my_Belem_PR, ds="steps-pre", ls=(0, (1, 1)))
hp = ax1.plot(Time, my_Helsinki_PR, ds="steps-pre", ls=(0, (1, 1)))
else:
ax1.fill_between(Time, my_Belem_PR - Belem_PR, alpha=0.5, step="pre")
ax1.fill_between(Time, my_Helsinki_PR - Helsinki_PR, alpha=0.5, step="pre")
bp = ax1.plot(Time, my_Belem_PR - Belem_PR, ds="steps-pre")
hp = ax1.plot(Time, my_Helsinki_PR - Helsinki_PR, ds="steps-pre")
if reference:
ax1.set_title("Original Observations\n")
ax2.set_title("Original Observation-space ERA5\n")
else:
err_Helsinki_v = np.mean(
~(
(np.abs(my_Helsinki_PR - Helsinki_PR) <= PR_eb_abs)
& (np.sign(my_Helsinki_PR) == np.sign(Helsinki_PR))
)
)
err_Helsinki_v = (
0
if err_Helsinki_v == 0
else np.format_float_positional(
100 * err_Helsinki_v, precision=1, min_digits=1
)
+ "%"
)
if err_Helsinki_v == "0.0%":
err_Helsinki_v = "<0.05%"
err_Belem_v = np.mean(
~(
(np.abs(my_Belem_PR - Belem_PR) <= PR_eb_abs)
& (np.sign(my_Belem_PR) == np.sign(Belem_PR))
)
)
err_Belem_v = (
0
if err_Belem_v == 0
else np.format_float_positional(
100 * err_Belem_v, precision=1, min_digits=1
)
+ "%"
)
if err_Belem_v == "0.0%":
err_Belem_v = "<0.05%"
err_era5_Helsinki_v = np.mean(
~(
(np.abs(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_Helsinki_PR) == np.sign(ERA5_Helsinki_PR))
)
)
err_era5_Helsinki_v = (
0
if err_era5_Helsinki_v == 0
else np.format_float_positional(
100 * err_era5_Helsinki_v, precision=1, min_digits=1
)
+ "%"
)
if err_era5_Helsinki_v == "0.0%":
err_era5_Helsinki_v = "<0.05%"
err_era5_Belem_v = np.mean(
~(
(np.abs(my_ERA5_Belem_PR - ERA5_Belem_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_Belem_PR) == np.sign(ERA5_Belem_PR))
)
)
err_era5_Belem_v = (
0
if err_era5_Belem_v == 0
else np.format_float_positional(
100 * err_era5_Belem_v, precision=1, min_digits=1
)
+ "%"
)
if err_era5_Belem_v == "0.0%":
err_era5_Belem_v = "<0.05%"
err_ERA5_v = np.mean(
~(
(np.abs(my_ERA5_PR - ERA5_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_PR) == np.sign(ERA5_PR))
)
)
err_ERA5_v = (
0
if err_ERA5_v == 0
else np.format_float_positional(100 * err_ERA5_v, precision=1, min_digits=1)
+ "%"
)
if err_ERA5_v == "0.0%":
err_ERA5_v = "<0.05%"
t = ax1.text(
0.05,
0.1,
f"V={err_Helsinki_v}",
ha="left",
va="bottom",
transform=ax1.transAxes,
color=adjust_color_lightness(hp[0].get_c()),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax1.text(
0.95,
0.1,
f"V={err_Belem_v}",
ha="right",
va="bottom",
transform=ax1.transAxes,
color=bp[0].get_c(),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax2.text(
0.05,
0.1,
f"V={err_era5_Helsinki_v}",
ha="left",
va="bottom",
transform=ax2.transAxes,
color=adjust_color_lightness(hp[0].get_c()),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax2.text(
0.95,
0.1,
f"V={err_era5_Belem_v}",
ha="right",
va="bottom",
transform=ax2.transAxes,
color=bp[0].get_c(),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax2.text(
0.5,
0.1,
f"V={err_ERA5_v}",
ha="center",
va="bottom",
transform=ax2.transAxes,
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
ax1.set_title(title)
ax2.set_title(title)
ax1p = ax1.inset_axes([0.0, -0.2, 1.0, 0.15], sharex=ax1)
ax1.tick_params(axis="x", labelbottom=False)
if reference:
ax1.set_ylabel("Precipitation (mm / 1h)")
else:
ax1.set_ylabel("Absolute error over precipitation (mm / 1h)")
ax1p.xaxis.set(major_locator=pos, major_formatter=fmt)
ax1p.set_xlim(Time.min(), Time.max())
ax1p.set_xticks(ax1p.get_xticks(), ax1p.get_xticklabels(), rotation=30, ha="right")
if reference:
ylim = (0, 25)
else:
ylim = max(abs(lim) for lim in ax1.get_ylim())
ylim = (-ylim, ylim)
ax1.set_ylim(*ylim)
plot_positive_precipitation(my_Belem_PR, ax1p, -1.75, bp[0].get_c(), lw=3)
plot_negative_precipitation(my_Belem_PR, ax1p, -1.75)
plot_positive_precipitation(my_Helsinki_PR, ax1p, 0.25, hp[0].get_c(), lw=3)
plot_negative_precipitation(my_Helsinki_PR, ax1p, 0.25)
plot_positive_precipitation(
Belem_PR, ax1p, -1.75, adjust_color_lightness(bp[0].get_c()), lw=1, bars=True
)
plot_positive_precipitation(
Helsinki_PR, ax1p, 0.25, adjust_color_lightness(hp[0].get_c()), lw=1, bars=True
)
ax1p.set_ylim(-3.5, 2.0)
ax1p.set_yticks([])
if reference:
ax1.text(
0.05,
1.09,
"Helsinki, Finland",
ha="left",
va="top",
transform=ax1.transAxes,
color=hp[0].get_c(),
size="large",
)
ax1.text(
0.95,
1.09,
"Belém, Brazil",
ha="right",
va="top",
transform=ax1.transAxes,
color=bp[0].get_c(),
size="large",
)
ax2.text(
0.05,
1.09,
"Helsinki, Finland",
ha="left",
va="top",
transform=ax2.transAxes,
color=hp[0].get_c(),
size="large",
)
ax2.text(
0.5,
1.09,
"Total",
ha="center",
va="top",
transform=ax2.transAxes,
size="large",
)
ax2.text(
0.95,
1.09,
"Belém, Brazil",
ha="right",
va="top",
transform=ax2.transAxes,
color=bp[0].get_c(),
size="large",
)
t = ax1.text(
0.05,
0.9,
humanize.naturalsize(Helsinki_PR.nbytes, binary=True)
if reference
else rf"$\times$ {np.round(my_Helsinki_PR_cr, 2)}",
ha="left",
va="top",
transform=ax1.transAxes,
color=adjust_color_lightness(hp[0].get_c()),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax1.text(
0.95,
0.9,
humanize.naturalsize(Belem_PR.nbytes, binary=True)
if reference
else rf"$\times$ {np.round(my_Belem_PR_cr, 2)}",
ha="right",
va="top",
transform=ax1.transAxes,
color=bp[0].get_c(),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
ax1.scatter(
[Time[np.argmax(my_Belem_PR)]],
[my_Belem_PR[np.argmax(my_Belem_PR)]]
if reference
else [(my_Belem_PR - Belem_PR)[np.argmax(my_Belem_PR)]],
marker="o",
facecolors="none",
edgecolors=bp[0].get_c(),
zorder=5,
)
ax1.scatter(
[Time[np.argmax(my_Helsinki_PR)]],
[my_Helsinki_PR[np.argmax(my_Helsinki_PR)]]
if reference
else [(my_Helsinki_PR - Helsinki_PR)[np.argmax(my_Helsinki_PR)]],
marker="o",
facecolors="none",
edgecolors=hp[0].get_c(),
zorder=5,
)
cax1 = ax1.inset_axes([1.025, 0.0, 0.1, 1.0], sharey=ax1)
cax1.hist(
[Helsinki_PR, Belem_PR]
if reference
else [my_Helsinki_PR - Helsinki_PR, my_Belem_PR - Belem_PR],
range=ylim,
bins=21,
orientation="horizontal",
color=[hp[0].get_c(), bp[0].get_c()],
histtype="stepfilled",
alpha=0.5,
)
cax1.hist(
[Helsinki_PR, Belem_PR]
if reference
else [my_Helsinki_PR - Helsinki_PR, my_Belem_PR - Belem_PR],
range=ylim,
bins=21,
orientation="horizontal",
color=[hp[0].get_c(), bp[0].get_c()],
histtype="step",
)
cax1.set_xticks([])
cax1.tick_params(axis="y", labelleft=False)
cax1.spines[:].set_visible(False)
cax1.patch.set_alpha(0.0)
q1, q2, q3 = np.quantile(my_Belem_PR - Belem_PR, [0.25, 0.5, 0.75])
cax1.axhline(
(my_Belem_PR - Belem_PR).mean(), ls=":", xmin=0.55, xmax=1.0, c="w", lw=2
)
cax1.axhline(q1, xmin=0.55, xmax=1.0, c="w", lw=2)
cax1.axhline(q2, xmin=0.75, xmax=1.0, c="w", lw=2)
cax1.axhline(q3, xmin=0.55, xmax=1.0, c="w", lw=2)
cax1.axhline(
(my_Belem_PR - Belem_PR).mean(),
xmin=0.55,
xmax=1.0,
ls=":",
c=bp[0].get_c(),
lw=1,
)
cax1.axhline(q1, xmin=0.55, xmax=1.0, c=bp[0].get_c(), lw=1)
cax1.axhline(q2, xmin=0.75, xmax=1.0, c=bp[0].get_c(), lw=1)
cax1.axhline(q3, xmin=0.55, xmax=1.0, c=bp[0].get_c(), lw=1)
q1, q2, q3 = np.quantile(my_Helsinki_PR - Helsinki_PR, [0.25, 0.5, 0.75])
cax1.axhline(
(my_Helsinki_PR - Helsinki_PR).mean(), ls=":", xmin=0.0, xmax=0.45, c="w", lw=2
)
cax1.axhline(q1, xmin=0.0, xmax=0.25, c="w", lw=2)
cax1.axhline(q2, xmin=0.0, xmax=0.45, c="w", lw=2)
cax1.axhline(q3, xmin=0.0, xmax=0.25, c="w", lw=2)
cax1.axhline(
(my_Helsinki_PR - Helsinki_PR).mean(),
xmin=0.0,
xmax=0.45,
ls=":",
c=hp[0].get_c(),
lw=1,
)
cax1.axhline(q1, xmin=0.0, xmax=0.25, c=hp[0].get_c(), lw=1)
cax1.axhline(q2, xmin=0.0, xmax=0.45, c=hp[0].get_c(), lw=1)
cax1.axhline(q3, xmin=0.0, xmax=0.25, c=hp[0].get_c(), lw=1)
if reference:
ax2.fill_between(Time, my_ERA5_Belem_PR, alpha=0.5, step="pre")
ax2.fill_between(Time, my_ERA5_Helsinki_PR, alpha=0.5, step="pre")
bp = ax2.plot(Time, my_ERA5_Belem_PR, ds="steps-pre", ls=(0, (1, 1)))
hp = ax2.plot(Time, my_ERA5_Helsinki_PR, ds="steps-pre", ls=(0, (1, 1)))
else:
ax2.fill_between(Time, my_ERA5_Belem_PR - ERA5_Belem_PR, alpha=0.5, step="pre")
ax2.fill_between(
Time, my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, alpha=0.5, step="pre"
)
bp = ax2.plot(Time, my_ERA5_Belem_PR - ERA5_Belem_PR, ds="steps-pre")
hp = ax2.plot(Time, my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, ds="steps-pre")
ax2p = ax2.inset_axes([0.0, -0.2, 1.0, 0.15], sharex=ax2)
ax2.tick_params(axis="x", labelbottom=False)
if reference:
ax2.set_ylabel("Precipitation (mm / 1h)")
else:
ax2.set_ylabel("Absolute error over precipitation (mm / 1h)")
ax2p.xaxis.set(major_locator=pos, major_formatter=fmt)
ax2p.set_xlim(Time.min(), Time.max())
ax2p.set_xticks(ax2p.get_xticks(), ax2p.get_xticklabels(), rotation=30, ha="right")
if reference:
ylim = (0, 25)
else:
ylim = max(abs(lim) for lim in ax2.get_ylim())
ylim = (-ylim, ylim)
ax2.set_ylim(*ylim)
plot_positive_precipitation(my_ERA5_Belem_PR, ax2p, -1.75, bp[0].get_c())
plot_negative_precipitation(my_ERA5_Belem_PR, ax2p, -1.75)
plot_positive_precipitation(my_ERA5_Helsinki_PR, ax2p, 0.25, hp[0].get_c())
plot_negative_precipitation(my_ERA5_Helsinki_PR, ax2p, 0.25)
plot_positive_precipitation(
ERA5_Belem_PR,
ax2p,
-1.75,
adjust_color_lightness(bp[0].get_c()),
lw=1,
bars=True,
)
plot_positive_precipitation(
ERA5_Helsinki_PR,
ax2p,
0.25,
adjust_color_lightness(hp[0].get_c()),
lw=1,
bars=True,
)
ax2p.set_ylim(-3, 1.5)
ax2p.set_yticks([])
t = ax2.text(
0.5,
0.9,
humanize.naturalsize(ERA5_PR.nbytes, binary=True)
if reference
else rf"$\times$ {np.round(my_ERA5_PR_cr, 2)}",
ha="center",
va="top",
transform=ax2.transAxes,
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
ax2.scatter(
[Time[np.argmax(my_ERA5_Belem_PR)]],
[my_ERA5_Belem_PR[np.argmax(my_ERA5_Belem_PR)]]
if reference
else [(my_ERA5_Belem_PR - ERA5_Belem_PR)[np.argmax(my_ERA5_Belem_PR)]],
marker="o",
facecolors="none",
edgecolors=bp[0].get_c(),
zorder=5,
)
ax2.scatter(
[Time[np.argmax(my_ERA5_Helsinki_PR)]],
[my_ERA5_Helsinki_PR[np.argmax(my_ERA5_Helsinki_PR)]]
if reference
else [(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR)[np.argmax(my_ERA5_Helsinki_PR)]],
marker="o",
facecolors="none",
edgecolors=hp[0].get_c(),
zorder=5,
)
cax2 = ax2.inset_axes([1.025, 0.0, 0.1, 1.0], sharey=ax2)
cax2.hist(
[ERA5_Helsinki_PR, ERA5_Belem_PR]
if reference
else [my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, my_ERA5_Belem_PR - ERA5_Belem_PR],
range=ylim,
bins=21,
orientation="horizontal",
color=[hp[0].get_c(), bp[0].get_c()],
histtype="stepfilled",
alpha=0.5,
)
cax2.hist(
[ERA5_Helsinki_PR, ERA5_Belem_PR]
if reference
else [my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, my_ERA5_Belem_PR - ERA5_Belem_PR],
range=ylim,
bins=21,
orientation="horizontal",
color=[hp[0].get_c(), bp[0].get_c()],
histtype="step",
)
cax2.set_xticks([])
cax2.tick_params(axis="y", labelleft=False)
cax2.spines[:].set_visible(False)
cax2.patch.set_alpha(0.0)
q1, q2, q3 = np.quantile(my_ERA5_Belem_PR - ERA5_Belem_PR, [0.25, 0.5, 0.75])
cax2.axhline(
(my_ERA5_Belem_PR - ERA5_Belem_PR).mean(),
ls=":",
xmin=0.55,
xmax=1.0,
c="w",
lw=2,
)
cax2.axhline(q1, xmin=0.75, xmax=1.0, c="w", lw=2)
cax2.axhline(q2, xmin=0.55, xmax=1.0, c="w", lw=2)
cax2.axhline(q3, xmin=0.75, xmax=1.0, c="w", lw=2)
cax2.axhline(
(my_ERA5_Belem_PR - ERA5_Belem_PR).mean(),
xmin=0.55,
xmax=1.0,
ls=":",
c=bp[0].get_c(),
lw=1,
)
cax2.axhline(q1, xmin=0.75, xmax=1.0, c=bp[0].get_c(), lw=1)
cax2.axhline(q2, xmin=0.55, xmax=1.0, c=bp[0].get_c(), lw=1)
cax2.axhline(q3, xmin=0.75, xmax=1.0, c=bp[0].get_c(), lw=1)
q1, q2, q3 = np.quantile(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, [0.25, 0.5, 0.75])
cax2.axhline(
(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR).mean(),
ls=":",
xmin=0.0,
xmax=0.45,
c="w",
lw=2,
)
cax2.axhline(q1, xmin=0.0, xmax=0.25, c="w", lw=2)
cax2.axhline(q2, xmin=0.0, xmax=0.45, c="w", lw=2)
cax2.axhline(q3, xmin=0.0, xmax=0.25, c="w", lw=2)
cax2.axhline(
(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR).mean(),
xmin=0.0,
xmax=0.45,
ls=":",
c=hp[0].get_c(),
lw=1,
)
cax2.axhline(q1, xmin=0.0, xmax=0.25, c=hp[0].get_c(), lw=1)
cax2.axhline(q2, xmin=0.0, xmax=0.45, c=hp[0].get_c(), lw=1)
cax2.axhline(q3, xmin=0.0, xmax=0.25, c=hp[0].get_c(), lw=1)
def plot_precipitation(
ax1,
ax2,
my_ERA5_PR: xr.DataArray,
my_ERA5_PR_cr: float,
my_Belem_PR: pd.DataFrame,
my_Belem_PR_cr: float,
my_Helsinki_PR: pd.DataFrame,
my_Helsinki_PR_cr: float,
PR_eb_abs: float,
title: str,
reference: bool = False,
corr=None,
):
my_ERA5_Belem_PR = my_ERA5_PR.sel(
latitude=-1.4563, longitude=360 - 48.5013, method="nearest"
).data
my_ERA5_Helsinki_PR = my_ERA5_PR.sel(
latitude=60.1699, longitude=24.9384, method="nearest"
).data
pos = mdates.HourLocator(byhour=(0, 6, 12, 18))
fmt = mdates.DateFormatter("%d.%m %Hh")
if reference:
ax1.fill_between(Time, my_Belem_PR, alpha=0.5, step="pre")
ax1.fill_between(Time, my_Helsinki_PR, alpha=0.5, step="pre")
bp = ax1.plot(Time, my_Belem_PR, ds="steps-pre", ls=(0, (1, 1)))
hp = ax1.plot(Time, my_Helsinki_PR, ds="steps-pre", ls=(0, (1, 1)))
else:
ax1.fill_between(Time, my_Belem_PR - Belem_PR, alpha=0.5, step="pre")
ax1.fill_between(Time, my_Helsinki_PR - Helsinki_PR, alpha=0.5, step="pre")
bp = ax1.plot(Time, my_Belem_PR - Belem_PR, ds="steps-pre")
hp = ax1.plot(Time, my_Helsinki_PR - Helsinki_PR, ds="steps-pre")
if reference:
ax1.set_title("Original Observations\n")
ax2.set_title("Original Observation-space ERA5\n")
else:
err_Helsinki_v = np.mean(
~(
(np.abs(my_Helsinki_PR - Helsinki_PR) <= PR_eb_abs)
& (np.sign(my_Helsinki_PR) == np.sign(Helsinki_PR))
)
)
err_Helsinki_v = (
0
if err_Helsinki_v == 0
else np.format_float_positional(
100 * err_Helsinki_v, precision=1, min_digits=1
)
+ "%"
)
if err_Helsinki_v == "0.0%":
err_Helsinki_v = "<0.05%"
err_Belem_v = np.mean(
~(
(np.abs(my_Belem_PR - Belem_PR) <= PR_eb_abs)
& (np.sign(my_Belem_PR) == np.sign(Belem_PR))
)
)
err_Belem_v = (
0
if err_Belem_v == 0
else np.format_float_positional(
100 * err_Belem_v, precision=1, min_digits=1
)
+ "%"
)
if err_Belem_v == "0.0%":
err_Belem_v = "<0.05%"
err_era5_Helsinki_v = np.mean(
~(
(np.abs(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_Helsinki_PR) == np.sign(ERA5_Helsinki_PR))
)
)
err_era5_Helsinki_v = (
0
if err_era5_Helsinki_v == 0
else np.format_float_positional(
100 * err_era5_Helsinki_v, precision=1, min_digits=1
)
+ "%"
)
if err_era5_Helsinki_v == "0.0%":
err_era5_Helsinki_v = "<0.05%"
err_era5_Belem_v = np.mean(
~(
(np.abs(my_ERA5_Belem_PR - ERA5_Belem_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_Belem_PR) == np.sign(ERA5_Belem_PR))
)
)
err_era5_Belem_v = (
0
if err_era5_Belem_v == 0
else np.format_float_positional(
100 * err_era5_Belem_v, precision=1, min_digits=1
)
+ "%"
)
if err_era5_Belem_v == "0.0%":
err_era5_Belem_v = "<0.05%"
err_ERA5_v = np.mean(
~(
(np.abs(my_ERA5_PR - ERA5_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_PR) == np.sign(ERA5_PR))
)
)
err_ERA5_v = (
0
if err_ERA5_v == 0
else np.format_float_positional(100 * err_ERA5_v, precision=1, min_digits=1)
+ "%"
)
if err_ERA5_v == "0.0%":
err_ERA5_v = "<0.05%"
t = ax1.text(
0.05,
0.1,
f"V={err_Helsinki_v}",
ha="left",
va="bottom",
transform=ax1.transAxes,
color=adjust_color_lightness(hp[0].get_c()),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax1.text(
0.95,
0.1,
f"V={err_Belem_v}",
ha="right",
va="bottom",
transform=ax1.transAxes,
color=bp[0].get_c(),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax2.text(
0.05,
0.1,
f"V={err_era5_Helsinki_v}",
ha="left",
va="bottom",
transform=ax2.transAxes,
color=adjust_color_lightness(hp[0].get_c()),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax2.text(
0.95,
0.1,
f"V={err_era5_Belem_v}",
ha="right",
va="bottom",
transform=ax2.transAxes,
color=bp[0].get_c(),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax2.text(
0.5,
0.1,
f"V={err_ERA5_v}",
ha="center",
va="bottom",
transform=ax2.transAxes,
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
ax1.set_title(title)
ax2.set_title(title)
ax1p = ax1.inset_axes([0.0, -0.2, 1.0, 0.15], sharex=ax1)
ax1.tick_params(axis="x", labelbottom=False)
if reference:
ax1.set_ylabel("Precipitation (mm / 1h)")
else:
ax1.set_ylabel("Absolute error over precipitation (mm / 1h)")
ax1p.xaxis.set(major_locator=pos, major_formatter=fmt)
ax1p.set_xlim(Time.min(), Time.max())
ax1p.set_xticks(ax1p.get_xticks(), ax1p.get_xticklabels(), rotation=30, ha="right")
if reference:
ylim = (0, 25)
else:
ylim = max(abs(lim) for lim in ax1.get_ylim())
ylim = (-ylim, ylim)
ax1.set_ylim(*ylim)
plot_positive_precipitation(my_Belem_PR, ax1p, -1.75, bp[0].get_c(), lw=3)
plot_negative_precipitation(my_Belem_PR, ax1p, -1.75)
plot_positive_precipitation(my_Helsinki_PR, ax1p, 0.25, hp[0].get_c(), lw=3)
plot_negative_precipitation(my_Helsinki_PR, ax1p, 0.25)
plot_positive_precipitation(
Belem_PR, ax1p, -1.75, adjust_color_lightness(bp[0].get_c()), lw=1, bars=True
)
plot_positive_precipitation(
Helsinki_PR, ax1p, 0.25, adjust_color_lightness(hp[0].get_c()), lw=1, bars=True
)
ax1p.set_ylim(-3.5, 2.0)
ax1p.set_yticks([])
if reference:
ax1.text(
0.05,
1.09,
"Helsinki, Finland",
ha="left",
va="top",
transform=ax1.transAxes,
color=hp[0].get_c(),
size="large",
)
ax1.text(
0.95,
1.09,
"Belém, Brazil",
ha="right",
va="top",
transform=ax1.transAxes,
color=bp[0].get_c(),
size="large",
)
ax2.text(
0.05,
1.09,
"Helsinki, Finland",
ha="left",
va="top",
transform=ax2.transAxes,
color=hp[0].get_c(),
size="large",
)
ax2.text(
0.5,
1.09,
"Total",
ha="center",
va="top",
transform=ax2.transAxes,
size="large",
)
ax2.text(
0.95,
1.09,
"Belém, Brazil",
ha="right",
va="top",
transform=ax2.transAxes,
color=bp[0].get_c(),
size="large",
)
t = ax1.text(
0.05,
0.9,
humanize.naturalsize(Helsinki_PR.nbytes, binary=True)
if reference
else rf"$\times$ {np.round(my_Helsinki_PR_cr, 2)}",
ha="left",
va="top",
transform=ax1.transAxes,
color=adjust_color_lightness(hp[0].get_c()),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
t = ax1.text(
0.95,
0.9,
humanize.naturalsize(Belem_PR.nbytes, binary=True)
if reference
else rf"$\times$ {np.round(my_Belem_PR_cr, 2)}",
ha="right",
va="top",
transform=ax1.transAxes,
color=bp[0].get_c(),
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
ax1.scatter(
[Time[np.argmax(my_Belem_PR)]],
[my_Belem_PR[np.argmax(my_Belem_PR)]]
if reference
else [(my_Belem_PR - Belem_PR)[np.argmax(my_Belem_PR)]],
marker="o",
facecolors="none",
edgecolors=bp[0].get_c(),
zorder=5,
)
ax1.scatter(
[Time[np.argmax(my_Helsinki_PR)]],
[my_Helsinki_PR[np.argmax(my_Helsinki_PR)]]
if reference
else [(my_Helsinki_PR - Helsinki_PR)[np.argmax(my_Helsinki_PR)]],
marker="o",
facecolors="none",
edgecolors=hp[0].get_c(),
zorder=5,
)
cax1 = ax1.inset_axes([1.025, 0.0, 0.1, 1.0], sharey=ax1)
cax1.hist(
[Helsinki_PR, Belem_PR]
if reference
else [my_Helsinki_PR - Helsinki_PR, my_Belem_PR - Belem_PR],
range=ylim,
bins=21,
orientation="horizontal",
color=[hp[0].get_c(), bp[0].get_c()],
histtype="stepfilled",
alpha=0.5,
)
cax1.hist(
[Helsinki_PR, Belem_PR]
if reference
else [my_Helsinki_PR - Helsinki_PR, my_Belem_PR - Belem_PR],
range=ylim,
bins=21,
orientation="horizontal",
color=[hp[0].get_c(), bp[0].get_c()],
histtype="step",
)
cax1.set_xticks([])
cax1.tick_params(axis="y", labelleft=False)
cax1.spines[:].set_visible(False)
cax1.patch.set_alpha(0.0)
q1, q2, q3 = np.quantile(my_Belem_PR - Belem_PR, [0.25, 0.5, 0.75])
cax1.axhline(
(my_Belem_PR - Belem_PR).mean(), ls=":", xmin=0.55, xmax=1.0, c="w", lw=2
)
cax1.axhline(q1, xmin=0.55, xmax=1.0, c="w", lw=2)
cax1.axhline(q2, xmin=0.75, xmax=1.0, c="w", lw=2)
cax1.axhline(q3, xmin=0.55, xmax=1.0, c="w", lw=2)
cax1.axhline(
(my_Belem_PR - Belem_PR).mean(),
xmin=0.55,
xmax=1.0,
ls=":",
c=bp[0].get_c(),
lw=1,
)
cax1.axhline(q1, xmin=0.55, xmax=1.0, c=bp[0].get_c(), lw=1)
cax1.axhline(q2, xmin=0.75, xmax=1.0, c=bp[0].get_c(), lw=1)
cax1.axhline(q3, xmin=0.55, xmax=1.0, c=bp[0].get_c(), lw=1)
q1, q2, q3 = np.quantile(my_Helsinki_PR - Helsinki_PR, [0.25, 0.5, 0.75])
cax1.axhline(
(my_Helsinki_PR - Helsinki_PR).mean(), ls=":", xmin=0.0, xmax=0.45, c="w", lw=2
)
cax1.axhline(q1, xmin=0.0, xmax=0.25, c="w", lw=2)
cax1.axhline(q2, xmin=0.0, xmax=0.45, c="w", lw=2)
cax1.axhline(q3, xmin=0.0, xmax=0.25, c="w", lw=2)
cax1.axhline(
(my_Helsinki_PR - Helsinki_PR).mean(),
xmin=0.0,
xmax=0.45,
ls=":",
c=hp[0].get_c(),
lw=1,
)
cax1.axhline(q1, xmin=0.0, xmax=0.25, c=hp[0].get_c(), lw=1)
cax1.axhline(q2, xmin=0.0, xmax=0.45, c=hp[0].get_c(), lw=1)
cax1.axhline(q3, xmin=0.0, xmax=0.25, c=hp[0].get_c(), lw=1)
if reference:
ax2.fill_between(Time, my_ERA5_Belem_PR, alpha=0.5, step="pre")
ax2.fill_between(Time, my_ERA5_Helsinki_PR, alpha=0.5, step="pre")
bp = ax2.plot(Time, my_ERA5_Belem_PR, ds="steps-pre", ls=(0, (1, 1)))
hp = ax2.plot(Time, my_ERA5_Helsinki_PR, ds="steps-pre", ls=(0, (1, 1)))
else:
ax2.fill_between(Time, my_ERA5_Belem_PR - ERA5_Belem_PR, alpha=0.5, step="pre")
ax2.fill_between(
Time, my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, alpha=0.5, step="pre"
)
bp = ax2.plot(Time, my_ERA5_Belem_PR - ERA5_Belem_PR, ds="steps-pre")
hp = ax2.plot(Time, my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, ds="steps-pre")
ax2p = ax2.inset_axes([0.0, -0.2, 1.0, 0.15], sharex=ax2)
ax2.tick_params(axis="x", labelbottom=False)
if reference:
ax2.set_ylabel("Precipitation (mm / 1h)")
else:
ax2.set_ylabel("Absolute error over precipitation (mm / 1h)")
ax2p.xaxis.set(major_locator=pos, major_formatter=fmt)
ax2p.set_xlim(Time.min(), Time.max())
ax2p.set_xticks(ax2p.get_xticks(), ax2p.get_xticklabels(), rotation=30, ha="right")
if reference:
ylim = (0, 25)
else:
ylim = max(abs(lim) for lim in ax2.get_ylim())
ylim = (-ylim, ylim)
ax2.set_ylim(*ylim)
plot_positive_precipitation(my_ERA5_Belem_PR, ax2p, -1.75, bp[0].get_c())
plot_negative_precipitation(my_ERA5_Belem_PR, ax2p, -1.75)
plot_positive_precipitation(my_ERA5_Helsinki_PR, ax2p, 0.25, hp[0].get_c())
plot_negative_precipitation(my_ERA5_Helsinki_PR, ax2p, 0.25)
plot_positive_precipitation(
ERA5_Belem_PR,
ax2p,
-1.75,
adjust_color_lightness(bp[0].get_c()),
lw=1,
bars=True,
)
plot_positive_precipitation(
ERA5_Helsinki_PR,
ax2p,
0.25,
adjust_color_lightness(hp[0].get_c()),
lw=1,
bars=True,
)
ax2p.set_ylim(-3, 1.5)
ax2p.set_yticks([])
t = ax2.text(
0.5,
0.9,
humanize.naturalsize(ERA5_PR.nbytes, binary=True)
if reference
else rf"$\times$ {np.round(my_ERA5_PR_cr, 2)}",
ha="center",
va="top",
transform=ax2.transAxes,
)
t.set_bbox(dict(facecolor="white", alpha=0.75, edgecolor="black"))
ax2.scatter(
[Time[np.argmax(my_ERA5_Belem_PR)]],
[my_ERA5_Belem_PR[np.argmax(my_ERA5_Belem_PR)]]
if reference
else [(my_ERA5_Belem_PR - ERA5_Belem_PR)[np.argmax(my_ERA5_Belem_PR)]],
marker="o",
facecolors="none",
edgecolors=bp[0].get_c(),
zorder=5,
)
ax2.scatter(
[Time[np.argmax(my_ERA5_Helsinki_PR)]],
[my_ERA5_Helsinki_PR[np.argmax(my_ERA5_Helsinki_PR)]]
if reference
else [(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR)[np.argmax(my_ERA5_Helsinki_PR)]],
marker="o",
facecolors="none",
edgecolors=hp[0].get_c(),
zorder=5,
)
cax2 = ax2.inset_axes([1.025, 0.0, 0.1, 1.0], sharey=ax2)
cax2.hist(
[ERA5_Helsinki_PR, ERA5_Belem_PR]
if reference
else [my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, my_ERA5_Belem_PR - ERA5_Belem_PR],
range=ylim,
bins=21,
orientation="horizontal",
color=[hp[0].get_c(), bp[0].get_c()],
histtype="stepfilled",
alpha=0.5,
)
cax2.hist(
[ERA5_Helsinki_PR, ERA5_Belem_PR]
if reference
else [my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, my_ERA5_Belem_PR - ERA5_Belem_PR],
range=ylim,
bins=21,
orientation="horizontal",
color=[hp[0].get_c(), bp[0].get_c()],
histtype="step",
)
cax2.set_xticks([])
cax2.tick_params(axis="y", labelleft=False)
cax2.spines[:].set_visible(False)
cax2.patch.set_alpha(0.0)
q1, q2, q3 = np.quantile(my_ERA5_Belem_PR - ERA5_Belem_PR, [0.25, 0.5, 0.75])
cax2.axhline(
(my_ERA5_Belem_PR - ERA5_Belem_PR).mean(),
ls=":",
xmin=0.55,
xmax=1.0,
c="w",
lw=2,
)
cax2.axhline(q1, xmin=0.75, xmax=1.0, c="w", lw=2)
cax2.axhline(q2, xmin=0.55, xmax=1.0, c="w", lw=2)
cax2.axhline(q3, xmin=0.75, xmax=1.0, c="w", lw=2)
cax2.axhline(
(my_ERA5_Belem_PR - ERA5_Belem_PR).mean(),
xmin=0.55,
xmax=1.0,
ls=":",
c=bp[0].get_c(),
lw=1,
)
cax2.axhline(q1, xmin=0.75, xmax=1.0, c=bp[0].get_c(), lw=1)
cax2.axhline(q2, xmin=0.55, xmax=1.0, c=bp[0].get_c(), lw=1)
cax2.axhline(q3, xmin=0.75, xmax=1.0, c=bp[0].get_c(), lw=1)
q1, q2, q3 = np.quantile(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR, [0.25, 0.5, 0.75])
cax2.axhline(
(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR).mean(),
ls=":",
xmin=0.0,
xmax=0.45,
c="w",
lw=2,
)
cax2.axhline(q1, xmin=0.0, xmax=0.25, c="w", lw=2)
cax2.axhline(q2, xmin=0.0, xmax=0.45, c="w", lw=2)
cax2.axhline(q3, xmin=0.0, xmax=0.25, c="w", lw=2)
cax2.axhline(
(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR).mean(),
xmin=0.0,
xmax=0.45,
ls=":",
c=hp[0].get_c(),
lw=1,
)
cax2.axhline(q1, xmin=0.0, xmax=0.25, c=hp[0].get_c(), lw=1)
cax2.axhline(q2, xmin=0.0, xmax=0.45, c=hp[0].get_c(), lw=1)
cax2.axhline(q3, xmin=0.0, xmax=0.25, c=hp[0].get_c(), lw=1)
Copied!
def table_precipitation(
my_ERA5_PR: xr.DataArray,
my_ERA5_PR_cr: float,
my_Belem_PR: pd.DataFrame,
my_Belem_PR_cr: float,
my_Helsinki_PR: pd.DataFrame,
my_Helsinki_PR_cr: float,
title: tuple[str, str, str],
PR_eb_abs: float,
corr: None | tuple[xr.DataArray, xr.DataArray, xr.DataArray],
reference: bool = False,
) -> pd.DataFrame:
my_ERA5_Belem_PR = my_ERA5_PR.sel(
latitude=-1.4563, longitude=360 - 48.5013, method="nearest"
).data
my_ERA5_Helsinki_PR = my_ERA5_PR.sel(
latitude=60.1699, longitude=24.9384, method="nearest"
).data
if reference:
err_Helsinki_inf = ""
err_Helsinki_2 = ""
err_Helsinki_v = ""
err_Belem_inf = ""
err_Belem_2 = ""
err_Belem_v = ""
err_era5_Helsinki_inf = ""
err_era5_Helsinki_2 = ""
err_era5_Helsinki_v = ""
err_era5_Belem_inf = ""
err_era5_Belem_2 = ""
err_era5_Belem_v = ""
err_ERA5_inf = ""
err_ERA5_2 = ""
err_ERA5_v = ""
corr_Helsinki = ""
corr_Belem = ""
corr_era5_Helsinki = ""
corr_era5_Belem = ""
corr_era5 = ""
else:
err_Helsinki_inf = np.amax(np.abs(my_Helsinki_PR - Helsinki_PR))
err_Helsinki_inf = f"{err_Helsinki_inf:.02}"
err_Helsinki_2 = np.sqrt(np.mean(np.square(my_Helsinki_PR - Helsinki_PR)))
err_Helsinki_2 = f"{err_Helsinki_2:.02}"
err_Helsinki_v = np.mean(
~(
(np.abs(my_Helsinki_PR - Helsinki_PR) <= PR_eb_abs)
& (np.sign(my_Helsinki_PR) == np.sign(Helsinki_PR))
)
)
err_Helsinki_v = (
0
if err_Helsinki_v == 0
else np.format_float_positional(
100 * err_Helsinki_v, precision=1, min_digits=1
)
+ "%"
)
if err_Helsinki_v == "0.0%":
err_Helsinki_v = "<0.05%"
err_Belem_inf = np.amax(np.abs(my_Belem_PR - Belem_PR))
err_Belem_inf = f"{err_Belem_inf:.02}"
err_Belem_2 = np.sqrt(np.mean(np.square(my_Belem_PR - Belem_PR)))
err_Belem_2 = f"{err_Belem_2:.02}"
err_Belem_v = np.mean(
~(
(np.abs(my_Belem_PR - Belem_PR) <= PR_eb_abs)
& (np.sign(my_Belem_PR) == np.sign(Belem_PR))
)
)
err_Belem_v = (
0
if err_Belem_v == 0
else np.format_float_positional(
100 * err_Belem_v, precision=1, min_digits=1
)
+ "%"
)
if err_Belem_v == "0.0%":
err_Belem_v = "<0.05%"
err_era5_Helsinki_inf = np.amax(np.abs(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR))
err_era5_Helsinki_inf = f"{err_era5_Helsinki_inf:.02}"
err_era5_Helsinki_2 = np.sqrt(
np.mean(np.square(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR))
)
err_era5_Helsinki_2 = f"{err_era5_Helsinki_2:.02}"
err_era5_Helsinki_v = np.mean(
~(
(np.abs(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_Helsinki_PR) == np.sign(ERA5_Helsinki_PR))
)
)
err_era5_Helsinki_v = (
0
if err_era5_Helsinki_v == 0
else np.format_float_positional(
100 * err_era5_Helsinki_v, precision=1, min_digits=1
)
+ "%"
)
if err_era5_Helsinki_v == "0.0%":
err_era5_Helsinki_v = "<0.05%"
err_era5_Belem_inf = np.amax(np.abs(my_ERA5_Belem_PR - ERA5_Belem_PR))
err_era5_Belem_inf = f"{err_era5_Belem_inf:.02}"
err_era5_Belem_2 = np.sqrt(np.mean(np.square(my_ERA5_Belem_PR - ERA5_Belem_PR)))
err_era5_Belem_2 = f"{err_era5_Belem_2:.02}"
err_era5_Belem_v = np.mean(
~(
(np.abs(my_ERA5_Belem_PR - ERA5_Belem_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_Belem_PR) == np.sign(ERA5_Belem_PR))
)
)
err_era5_Belem_v = (
0
if err_era5_Belem_v == 0
else np.format_float_positional(
100 * err_era5_Belem_v, precision=1, min_digits=1
)
+ "%"
)
if err_era5_Belem_v == "0.0%":
err_era5_Belem_v = "<0.05%"
err_ERA5_inf = np.amax(np.abs(my_ERA5_PR - ERA5_PR))
err_ERA5_inf = f"{err_ERA5_inf:.02}"
err_ERA5_2 = np.sqrt(np.mean(np.square(my_ERA5_PR - ERA5_PR)))
err_ERA5_2 = f"{err_ERA5_2:.02}"
err_ERA5_v = np.mean(
~(
(np.abs(my_ERA5_PR - ERA5_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_PR) == np.sign(ERA5_PR))
)
)
err_ERA5_v = (
0
if err_ERA5_v == 0
else np.format_float_positional(100 * err_ERA5_v, precision=1, min_digits=1)
+ "%"
)
if err_ERA5_v == "0.0%":
err_ERA5_v = "<0.05%"
if corr is not None:
corr_ERA5, corr_Belem, corr_Helsinki = corr
corr_Helsinki = compute_corrections_percentage(
my_Helsinki_PR, corr_Helsinki
)
corr_Helsinki = (
0
if corr_Helsinki == 0
else np.format_float_positional(
100 * corr_Helsinki, precision=1, min_digits=1
)
+ "%"
)
if corr_Helsinki == "0.0%":
corr_Helsinki = "<0.05%"
corr_Belem = compute_corrections_percentage(my_Belem_PR, corr_Belem)
corr_Belem = (
0
if corr_Belem == 0
else np.format_float_positional(
100 * corr_Belem, precision=1, min_digits=1
)
+ "%"
)
if corr_Belem == "0.0%":
corr_Belem = "<0.05%"
corr_era5_Helsinki = compute_corrections_percentage(
my_ERA5_Helsinki_PR,
corr_ERA5.sel(
latitude=60.1699, longitude=24.9384, method="nearest"
).data,
)
corr_era5_Helsinki = (
0
if corr_era5_Helsinki == 0
else np.format_float_positional(
100 * corr_era5_Helsinki, precision=1, min_digits=1
)
+ "%"
)
if corr_era5_Helsinki == "0.0%":
corr_era5_Helsinki = "<0.05%"
corr_era5_Belem = compute_corrections_percentage(
my_ERA5_Belem_PR,
corr_ERA5.sel(
latitude=-1.4563, longitude=360 - 48.5013, method="nearest"
).data,
)
corr_era5_Belem = (
0
if corr_era5_Belem == 0
else np.format_float_positional(
100 * corr_era5_Belem, precision=1, min_digits=1
)
+ "%"
)
if corr_era5_Belem == "0.0%":
corr_era5_Belem = "<0.05%"
corr_era5 = compute_corrections_percentage(my_ERA5_PR, corr_ERA5)
corr_era5 = (
0
if corr_era5 == 0
else np.format_float_positional(
100 * corr_era5, precision=1, min_digits=1
)
+ "%"
)
if corr_era5 == "0.0%":
corr_era5 = "<0.05%"
else:
corr_Helsinki = ""
corr_Belem = ""
corr_era5_Helsinki = ""
corr_era5_Belem = ""
corr_era5 = ""
if reference:
integral_belem = (
np.format_float_positional(np.sum(Belem_PR), precision=1, min_digits=1)
+ "mm"
)
integral_era5_belem = (
np.format_float_positional(np.sum(ERA5_Belem_PR), precision=1, min_digits=1)
+ "mm"
)
integral_helsinki = (
np.format_float_positional(np.sum(Helsinki_PR), precision=1, min_digits=1)
+ "mm"
)
integral_era5_helsinki = (
np.format_float_positional(
np.sum(ERA5_Helsinki_PR), precision=1, min_digits=1
)
+ "mm"
)
integral_era5 = (
np.format_float_positional(
np.mean(ERA5_PR.sum("valid_time")), precision=1, min_digits=1
)
+ "mm"
)
else:
integral_belem = (
100 * (np.sum(my_Belem_PR) - np.sum(Belem_PR)) / np.sum(Belem_PR)
)
integral_belem = (
"+0mm"
if integral_belem == 0
else np.format_float_positional(
integral_belem, precision=1, min_digits=1, sign=True
)
+ "%"
)
integral_era5_belem = (
100
* (np.sum(my_ERA5_Belem_PR) - np.sum(ERA5_Belem_PR))
/ np.sum(ERA5_Belem_PR)
)
integral_era5_belem = (
"+0mm"
if integral_era5_belem == 0
else np.format_float_positional(
integral_era5_belem, precision=1, min_digits=1, sign=True
)
+ "%"
)
integral_helsinki = (
100 * (np.sum(my_Helsinki_PR) - np.sum(Helsinki_PR)) / np.sum(Helsinki_PR)
)
integral_helsinki = (
"+0mm"
if integral_helsinki == 0
else np.format_float_positional(
integral_helsinki, precision=1, min_digits=1, sign=True
)
+ "%"
)
integral_era5_helsinki = (
100
* (np.sum(my_ERA5_Helsinki_PR) - np.sum(ERA5_Helsinki_PR))
/ np.sum(ERA5_Helsinki_PR)
)
integral_era5_helsinki = (
"+0mm"
if integral_era5_helsinki == 0
else np.format_float_positional(
integral_era5_helsinki, precision=1, min_digits=1, sign=True
)
+ "%"
)
integral_era5 = 100 * (np.sum(my_ERA5_PR) - np.sum(ERA5_PR)) / np.sum(ERA5_PR)
integral_era5 = (
"+0mm"
if integral_era5 == 0
else np.format_float_positional(
integral_era5, precision=1, min_digits=1, sign=True
)
+ "%"
)
if reference:
maxv_belem = (
np.format_float_positional(np.amax(Belem_PR), precision=2, min_digits=2)
+ "mm"
)
maxv_era5_belem = (
np.format_float_positional(
np.amax(ERA5_Belem_PR), precision=2, min_digits=2
)
+ "mm"
)
maxv_helsinki = (
np.format_float_positional(np.amax(Helsinki_PR), precision=2, min_digits=2)
+ "mm"
)
maxv_era5_helsinki = (
np.format_float_positional(
np.amax(ERA5_Helsinki_PR), precision=2, min_digits=2
)
+ "mm"
)
maxv_era5 = (
np.format_float_positional(
np.mean(ERA5_PR.max("valid_time")), precision=2, min_digits=2
)
+ "mm"
)
else:
maxv_belem = np.amax(my_Belem_PR) - np.amax(Belem_PR)
maxv_belem = (
"+0mm"
if maxv_belem == 0
else np.format_float_positional(
maxv_belem, precision=2, min_digits=2, sign=True
)
+ "mm"
)
maxv_era5_belem = np.amax(my_ERA5_Belem_PR) - np.amax(ERA5_Belem_PR)
maxv_era5_belem = (
"+0mm"
if maxv_era5_belem == 0
else np.format_float_positional(
maxv_era5_belem, precision=2, min_digits=2, sign=True
)
+ "mm"
)
maxv_helsinki = np.amax(my_Helsinki_PR) - np.amax(Helsinki_PR)
maxv_helsinki = (
"+0mm"
if maxv_helsinki == 0
else np.format_float_positional(
maxv_helsinki, precision=2, min_digits=2, sign=True
)
+ "mm"
)
maxv_era5_helsinki = np.amax(my_ERA5_Helsinki_PR) - np.amax(ERA5_Helsinki_PR)
maxv_era5_helsinki = (
"+0mm"
if maxv_era5_helsinki == 0
else np.format_float_positional(
maxv_era5_helsinki, precision=2, min_digits=2, sign=True
)
+ "mm"
)
maxv_era5 = np.mean(my_ERA5_PR.max("valid_time") - ERA5_PR.max("valid_time"))
maxv_era5 = (
"+0mm"
if maxv_era5 == 0
else np.format_float_positional(
maxv_era5, precision=2, min_digits=2, sign=True
)
)
if reference:
argmax_belem = pd.to_datetime(str(Time[np.argmax(Belem_PR)])).strftime(
"%d.%m %Hh"
)
argmax_era5_belem = pd.to_datetime(
str(Time[np.argmax(ERA5_Belem_PR)])
).strftime("%d.%m %Hh")
argmax_helsinki = pd.to_datetime(str(Time[np.argmax(Helsinki_PR)])).strftime(
"%d.%m %Hh"
)
argmax_era5_helsinki = pd.to_datetime(
str(Time[np.argmax(ERA5_Helsinki_PR)])
).strftime("%d.%m %Hh")
argmax_era5 = ""
else:
argmax_belem = (
Time[np.argmax(my_Belem_PR)] - Time[np.argmax(Belem_PR)]
) / np.timedelta64(1, "h")
argmax_belem = (
np.format_float_positional(argmax_belem, precision=0, sign=True, trim="-")
+ "h"
)
argmax_era5_belem = (
Time[np.argmax(my_ERA5_Belem_PR)] - Time[np.argmax(ERA5_Belem_PR)]
) / np.timedelta64(1, "h")
argmax_era5_belem = (
np.format_float_positional(
argmax_era5_belem, precision=0, sign=True, trim="-"
)
+ "h"
)
argmax_helsinki = (
Time[np.argmax(my_Helsinki_PR)] - Time[np.argmax(Helsinki_PR)]
) / np.timedelta64(1, "h")
argmax_helsinki = (
np.format_float_positional(
argmax_helsinki, precision=0, sign=True, trim="-"
)
+ "h"
)
argmax_era5_helsinki = (
Time[np.argmax(my_ERA5_Helsinki_PR)] - Time[np.argmax(ERA5_Helsinki_PR)]
) / np.timedelta64(1, "h")
argmax_era5_helsinki = (
np.format_float_positional(
argmax_era5_helsinki, precision=0, sign=True, trim="-"
)
+ "h"
)
argmax_era5 = (
np.mean(my_ERA5_PR.argmax("valid_time") - ERA5_PR.argmax("valid_time"))
* (Time[1] - Time[0])
/ np.timedelta64(1, "h")
)
argmax_era5 = (
np.format_float_positional(
argmax_era5, precision=1, min_digits=1, sign=True, trim="-"
)
+ "h"
)
if reference:
positive_belem = f"{np.sum(Belem_PR > 0)}" + "h"
positive_era5_belem = f"{np.sum(ERA5_Belem_PR > 0)}" + "h"
positive_helsinki = f"{np.sum(Helsinki_PR > 0)}" + "h"
positive_era5_helsinki = f"{np.sum(ERA5_Helsinki_PR > 0)}" + "h"
positive_era5 = (
np.format_float_positional(np.mean(ERA5_PR > 0), precision=2, min_digits=2)
+ "%"
)
else:
positive_belem = (
100
* (np.sum(my_Belem_PR > 0) - np.sum(Belem_PR > 0))
/ np.sum(Belem_PR > 0)
)
positive_belem = (
"+0h"
if positive_belem == 0
else np.format_float_positional(
positive_belem, precision=1, min_digits=1, sign=True
)
+ "%"
)
positive_era5_belem = (
100
* (np.sum(my_ERA5_Belem_PR > 0) - np.sum(ERA5_Belem_PR > 0))
/ np.sum(ERA5_Belem_PR > 0)
)
positive_era5_belem = (
"+0h"
if positive_era5_belem == 0
else np.format_float_positional(
positive_era5_belem, precision=1, min_digits=1, sign=True
)
+ "%"
)
positive_helsinki = (
100
* (np.sum(my_Helsinki_PR > 0) - np.sum(Helsinki_PR > 0))
/ np.sum(Helsinki_PR > 0)
)
positive_helsinki = (
"+0h"
if positive_helsinki == 0
else np.format_float_positional(
positive_helsinki, precision=1, min_digits=1, sign=True
)
+ "%"
)
positive_era5_helsinki = (
100
* (np.sum(my_ERA5_Helsinki_PR > 0) - np.sum(ERA5_Helsinki_PR > 0))
/ np.sum(ERA5_Helsinki_PR > 0)
)
positive_era5_helsinki = (
"+0h"
if positive_era5_helsinki == 0
else np.format_float_positional(
positive_era5_helsinki, precision=1, min_digits=1, sign=True
)
+ "%"
)
positive_era5 = (
np.format_float_positional(
np.mean(my_ERA5_PR > 0), precision=2, min_digits=2
)
+ "%"
)
if reference:
fpfn_belem = "0h"
fpfn_era5_belem = "0h"
fpfn_helsinki = "0h"
fpfn_era5_helsinki = "0h"
fpfn_era5 = "0h"
else:
fpfn_belem = 100 * np.mean((my_Belem_PR > 0) != (Belem_PR > 0))
fpfn_belem = (
"0h"
if fpfn_belem == 0
else np.format_float_positional(fpfn_belem, precision=1, min_digits=1) + "%"
)
fpfn_era5_belem = 100 * np.mean((my_ERA5_Belem_PR > 0) != (ERA5_Belem_PR > 0))
fpfn_era5_belem = (
"0h"
if fpfn_era5_belem == 0
else np.format_float_positional(fpfn_era5_belem, precision=1, min_digits=1)
+ "%"
)
fpfn_helsinki = 100 * np.mean((my_Helsinki_PR > 0) != (Helsinki_PR > 0))
fpfn_helsinki = (
"0h"
if fpfn_helsinki == 0
else np.format_float_positional(fpfn_helsinki, precision=1, min_digits=1)
+ "%"
)
fpfn_era5_helsinki = 100 * np.mean(
(my_ERA5_Helsinki_PR > 0) != (ERA5_Helsinki_PR > 0)
)
fpfn_era5_helsinki = (
"0h"
if fpfn_era5_helsinki == 0
else np.format_float_positional(
fpfn_era5_helsinki, precision=1, min_digits=1
)
+ "%"
)
fpfn_era5 = 100 * np.mean((my_ERA5_PR > 0) != (ERA5_PR > 0))
fpfn_era5 = (
"0h"
if fpfn_era5 == 0
else np.format_float_positional(fpfn_era5, precision=1, min_digits=1) + "%"
)
if reference:
negative_belem = f"{np.sum(Belem_PR < 0)}" + "h"
negative_era5_belem = f"{np.sum(ERA5_Belem_PR < 0)}" + "h"
negative_helsinki = f"{np.sum(Helsinki_PR < 0)}" + "h"
negative_era5_helsinki = f"{np.sum(ERA5_Helsinki_PR < 0)}" + "h"
negative_era5 = f"{int(np.sum(ERA5_PR < 0))}" + "h"
else:
negative_belem = 100 * np.mean(my_Belem_PR < 0)
negative_belem = (
"0h"
if negative_belem == 0
else np.format_float_positional(negative_belem, precision=1, min_digits=1)
+ "%"
)
negative_era5_belem = 100 * np.mean(my_ERA5_Belem_PR < 0)
negative_era5_belem = (
"0h"
if negative_era5_belem == 0
else np.format_float_positional(
negative_era5_belem, precision=1, min_digits=1
)
+ "%"
)
negative_helsinki = 100 * np.mean(my_Helsinki_PR < 0)
negative_helsinki = (
"0h"
if negative_helsinki == 0
else np.format_float_positional(
negative_helsinki, precision=1, min_digits=1
)
+ "%"
)
negative_era5_helsinki = 100 * np.mean(my_ERA5_Helsinki_PR < 0)
negative_era5_helsinki = (
"0h"
if negative_era5_helsinki == 0
else np.format_float_positional(
negative_era5_helsinki, precision=1, min_digits=1
)
+ "%"
)
negative_era5 = 100 * np.mean(my_ERA5_PR < 0)
negative_era5 = (
"0h"
if negative_era5 == 0
else np.format_float_positional(negative_era5, precision=1, min_digits=1)
+ "%"
)
if reference:
cr_helsinki = r"$\times$ 1"
cr_belem = r"$\times$ 1"
cr_era5 = r"$\times$ 1"
else:
cr_helsinki = rf"$\times$ {np.round(my_Helsinki_PR_cr, 2)}"
cr_belem = rf"$\times$ {np.round(my_Belem_PR_cr, 2)}"
cr_era5 = rf"$\times$ {np.round(my_ERA5_PR_cr, 2)}"
return pd.DataFrame(
{
"Compressor": [title[0], title[0], title[0], title[0], title[0]],
"Safeguarded": [title[1], title[1], title[1], title[1], title[1]],
"Corrections": [title[2], title[2], title[2], title[2], title[2]],
"Station": ["Helsinki", "Helsinki", "Belém", "Belém", "ERA5"],
"Source": ["obs", "ERA5", "obs", "ERA5", "avg"],
r"$L_{\infty}(\hat{PR})$": [
err_Helsinki_inf,
err_era5_Helsinki_inf,
err_Belem_inf,
err_era5_Belem_inf,
err_ERA5_inf,
],
r"$L_{2}(\hat{PR})$": [
err_Helsinki_2,
err_era5_Helsinki_2,
err_Belem_2,
err_era5_Belem_2,
err_ERA5_2,
],
"Integral": [
integral_helsinki,
integral_era5_helsinki,
integral_belem,
integral_era5_belem,
integral_era5,
],
"max(PR)": [
maxv_helsinki,
maxv_era5_helsinki,
maxv_belem,
maxv_era5_belem,
maxv_era5,
],
"argmax(PR)": [
argmax_helsinki,
argmax_era5_helsinki,
argmax_belem,
argmax_era5_belem,
argmax_era5,
],
"PR > 0": [
positive_helsinki,
positive_era5_helsinki,
positive_belem,
positive_era5_belem,
positive_era5,
],
"FP + FN": [
fpfn_helsinki,
fpfn_era5_helsinki,
fpfn_belem,
fpfn_era5_belem,
fpfn_era5,
],
"PR < 0": [
negative_helsinki,
negative_era5_helsinki,
negative_belem,
negative_era5_belem,
negative_era5,
],
"V": [
err_Helsinki_v,
err_era5_Helsinki_v,
err_Belem_v,
err_era5_Belem_v,
err_ERA5_v,
],
"C": [
corr_Helsinki,
corr_era5_Helsinki,
corr_Belem,
corr_era5_Belem,
corr_era5,
],
"CR": [cr_helsinki, "", cr_belem, "", cr_era5],
}
)
def table_precipitation(
my_ERA5_PR: xr.DataArray,
my_ERA5_PR_cr: float,
my_Belem_PR: pd.DataFrame,
my_Belem_PR_cr: float,
my_Helsinki_PR: pd.DataFrame,
my_Helsinki_PR_cr: float,
title: tuple[str, str, str],
PR_eb_abs: float,
corr: None | tuple[xr.DataArray, xr.DataArray, xr.DataArray],
reference: bool = False,
) -> pd.DataFrame:
my_ERA5_Belem_PR = my_ERA5_PR.sel(
latitude=-1.4563, longitude=360 - 48.5013, method="nearest"
).data
my_ERA5_Helsinki_PR = my_ERA5_PR.sel(
latitude=60.1699, longitude=24.9384, method="nearest"
).data
if reference:
err_Helsinki_inf = ""
err_Helsinki_2 = ""
err_Helsinki_v = ""
err_Belem_inf = ""
err_Belem_2 = ""
err_Belem_v = ""
err_era5_Helsinki_inf = ""
err_era5_Helsinki_2 = ""
err_era5_Helsinki_v = ""
err_era5_Belem_inf = ""
err_era5_Belem_2 = ""
err_era5_Belem_v = ""
err_ERA5_inf = ""
err_ERA5_2 = ""
err_ERA5_v = ""
corr_Helsinki = ""
corr_Belem = ""
corr_era5_Helsinki = ""
corr_era5_Belem = ""
corr_era5 = ""
else:
err_Helsinki_inf = np.amax(np.abs(my_Helsinki_PR - Helsinki_PR))
err_Helsinki_inf = f"{err_Helsinki_inf:.02}"
err_Helsinki_2 = np.sqrt(np.mean(np.square(my_Helsinki_PR - Helsinki_PR)))
err_Helsinki_2 = f"{err_Helsinki_2:.02}"
err_Helsinki_v = np.mean(
~(
(np.abs(my_Helsinki_PR - Helsinki_PR) <= PR_eb_abs)
& (np.sign(my_Helsinki_PR) == np.sign(Helsinki_PR))
)
)
err_Helsinki_v = (
0
if err_Helsinki_v == 0
else np.format_float_positional(
100 * err_Helsinki_v, precision=1, min_digits=1
)
+ "%"
)
if err_Helsinki_v == "0.0%":
err_Helsinki_v = "<0.05%"
err_Belem_inf = np.amax(np.abs(my_Belem_PR - Belem_PR))
err_Belem_inf = f"{err_Belem_inf:.02}"
err_Belem_2 = np.sqrt(np.mean(np.square(my_Belem_PR - Belem_PR)))
err_Belem_2 = f"{err_Belem_2:.02}"
err_Belem_v = np.mean(
~(
(np.abs(my_Belem_PR - Belem_PR) <= PR_eb_abs)
& (np.sign(my_Belem_PR) == np.sign(Belem_PR))
)
)
err_Belem_v = (
0
if err_Belem_v == 0
else np.format_float_positional(
100 * err_Belem_v, precision=1, min_digits=1
)
+ "%"
)
if err_Belem_v == "0.0%":
err_Belem_v = "<0.05%"
err_era5_Helsinki_inf = np.amax(np.abs(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR))
err_era5_Helsinki_inf = f"{err_era5_Helsinki_inf:.02}"
err_era5_Helsinki_2 = np.sqrt(
np.mean(np.square(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR))
)
err_era5_Helsinki_2 = f"{err_era5_Helsinki_2:.02}"
err_era5_Helsinki_v = np.mean(
~(
(np.abs(my_ERA5_Helsinki_PR - ERA5_Helsinki_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_Helsinki_PR) == np.sign(ERA5_Helsinki_PR))
)
)
err_era5_Helsinki_v = (
0
if err_era5_Helsinki_v == 0
else np.format_float_positional(
100 * err_era5_Helsinki_v, precision=1, min_digits=1
)
+ "%"
)
if err_era5_Helsinki_v == "0.0%":
err_era5_Helsinki_v = "<0.05%"
err_era5_Belem_inf = np.amax(np.abs(my_ERA5_Belem_PR - ERA5_Belem_PR))
err_era5_Belem_inf = f"{err_era5_Belem_inf:.02}"
err_era5_Belem_2 = np.sqrt(np.mean(np.square(my_ERA5_Belem_PR - ERA5_Belem_PR)))
err_era5_Belem_2 = f"{err_era5_Belem_2:.02}"
err_era5_Belem_v = np.mean(
~(
(np.abs(my_ERA5_Belem_PR - ERA5_Belem_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_Belem_PR) == np.sign(ERA5_Belem_PR))
)
)
err_era5_Belem_v = (
0
if err_era5_Belem_v == 0
else np.format_float_positional(
100 * err_era5_Belem_v, precision=1, min_digits=1
)
+ "%"
)
if err_era5_Belem_v == "0.0%":
err_era5_Belem_v = "<0.05%"
err_ERA5_inf = np.amax(np.abs(my_ERA5_PR - ERA5_PR))
err_ERA5_inf = f"{err_ERA5_inf:.02}"
err_ERA5_2 = np.sqrt(np.mean(np.square(my_ERA5_PR - ERA5_PR)))
err_ERA5_2 = f"{err_ERA5_2:.02}"
err_ERA5_v = np.mean(
~(
(np.abs(my_ERA5_PR - ERA5_PR) <= PR_eb_abs)
& (np.sign(my_ERA5_PR) == np.sign(ERA5_PR))
)
)
err_ERA5_v = (
0
if err_ERA5_v == 0
else np.format_float_positional(100 * err_ERA5_v, precision=1, min_digits=1)
+ "%"
)
if err_ERA5_v == "0.0%":
err_ERA5_v = "<0.05%"
if corr is not None:
corr_ERA5, corr_Belem, corr_Helsinki = corr
corr_Helsinki = compute_corrections_percentage(
my_Helsinki_PR, corr_Helsinki
)
corr_Helsinki = (
0
if corr_Helsinki == 0
else np.format_float_positional(
100 * corr_Helsinki, precision=1, min_digits=1
)
+ "%"
)
if corr_Helsinki == "0.0%":
corr_Helsinki = "<0.05%"
corr_Belem = compute_corrections_percentage(my_Belem_PR, corr_Belem)
corr_Belem = (
0
if corr_Belem == 0
else np.format_float_positional(
100 * corr_Belem, precision=1, min_digits=1
)
+ "%"
)
if corr_Belem == "0.0%":
corr_Belem = "<0.05%"
corr_era5_Helsinki = compute_corrections_percentage(
my_ERA5_Helsinki_PR,
corr_ERA5.sel(
latitude=60.1699, longitude=24.9384, method="nearest"
).data,
)
corr_era5_Helsinki = (
0
if corr_era5_Helsinki == 0
else np.format_float_positional(
100 * corr_era5_Helsinki, precision=1, min_digits=1
)
+ "%"
)
if corr_era5_Helsinki == "0.0%":
corr_era5_Helsinki = "<0.05%"
corr_era5_Belem = compute_corrections_percentage(
my_ERA5_Belem_PR,
corr_ERA5.sel(
latitude=-1.4563, longitude=360 - 48.5013, method="nearest"
).data,
)
corr_era5_Belem = (
0
if corr_era5_Belem == 0
else np.format_float_positional(
100 * corr_era5_Belem, precision=1, min_digits=1
)
+ "%"
)
if corr_era5_Belem == "0.0%":
corr_era5_Belem = "<0.05%"
corr_era5 = compute_corrections_percentage(my_ERA5_PR, corr_ERA5)
corr_era5 = (
0
if corr_era5 == 0
else np.format_float_positional(
100 * corr_era5, precision=1, min_digits=1
)
+ "%"
)
if corr_era5 == "0.0%":
corr_era5 = "<0.05%"
else:
corr_Helsinki = ""
corr_Belem = ""
corr_era5_Helsinki = ""
corr_era5_Belem = ""
corr_era5 = ""
if reference:
integral_belem = (
np.format_float_positional(np.sum(Belem_PR), precision=1, min_digits=1)
+ "mm"
)
integral_era5_belem = (
np.format_float_positional(np.sum(ERA5_Belem_PR), precision=1, min_digits=1)
+ "mm"
)
integral_helsinki = (
np.format_float_positional(np.sum(Helsinki_PR), precision=1, min_digits=1)
+ "mm"
)
integral_era5_helsinki = (
np.format_float_positional(
np.sum(ERA5_Helsinki_PR), precision=1, min_digits=1
)
+ "mm"
)
integral_era5 = (
np.format_float_positional(
np.mean(ERA5_PR.sum("valid_time")), precision=1, min_digits=1
)
+ "mm"
)
else:
integral_belem = (
100 * (np.sum(my_Belem_PR) - np.sum(Belem_PR)) / np.sum(Belem_PR)
)
integral_belem = (
"+0mm"
if integral_belem == 0
else np.format_float_positional(
integral_belem, precision=1, min_digits=1, sign=True
)
+ "%"
)
integral_era5_belem = (
100
* (np.sum(my_ERA5_Belem_PR) - np.sum(ERA5_Belem_PR))
/ np.sum(ERA5_Belem_PR)
)
integral_era5_belem = (
"+0mm"
if integral_era5_belem == 0
else np.format_float_positional(
integral_era5_belem, precision=1, min_digits=1, sign=True
)
+ "%"
)
integral_helsinki = (
100 * (np.sum(my_Helsinki_PR) - np.sum(Helsinki_PR)) / np.sum(Helsinki_PR)
)
integral_helsinki = (
"+0mm"
if integral_helsinki == 0
else np.format_float_positional(
integral_helsinki, precision=1, min_digits=1, sign=True
)
+ "%"
)
integral_era5_helsinki = (
100
* (np.sum(my_ERA5_Helsinki_PR) - np.sum(ERA5_Helsinki_PR))
/ np.sum(ERA5_Helsinki_PR)
)
integral_era5_helsinki = (
"+0mm"
if integral_era5_helsinki == 0
else np.format_float_positional(
integral_era5_helsinki, precision=1, min_digits=1, sign=True
)
+ "%"
)
integral_era5 = 100 * (np.sum(my_ERA5_PR) - np.sum(ERA5_PR)) / np.sum(ERA5_PR)
integral_era5 = (
"+0mm"
if integral_era5 == 0
else np.format_float_positional(
integral_era5, precision=1, min_digits=1, sign=True
)
+ "%"
)
if reference:
maxv_belem = (
np.format_float_positional(np.amax(Belem_PR), precision=2, min_digits=2)
+ "mm"
)
maxv_era5_belem = (
np.format_float_positional(
np.amax(ERA5_Belem_PR), precision=2, min_digits=2
)
+ "mm"
)
maxv_helsinki = (
np.format_float_positional(np.amax(Helsinki_PR), precision=2, min_digits=2)
+ "mm"
)
maxv_era5_helsinki = (
np.format_float_positional(
np.amax(ERA5_Helsinki_PR), precision=2, min_digits=2
)
+ "mm"
)
maxv_era5 = (
np.format_float_positional(
np.mean(ERA5_PR.max("valid_time")), precision=2, min_digits=2
)
+ "mm"
)
else:
maxv_belem = np.amax(my_Belem_PR) - np.amax(Belem_PR)
maxv_belem = (
"+0mm"
if maxv_belem == 0
else np.format_float_positional(
maxv_belem, precision=2, min_digits=2, sign=True
)
+ "mm"
)
maxv_era5_belem = np.amax(my_ERA5_Belem_PR) - np.amax(ERA5_Belem_PR)
maxv_era5_belem = (
"+0mm"
if maxv_era5_belem == 0
else np.format_float_positional(
maxv_era5_belem, precision=2, min_digits=2, sign=True
)
+ "mm"
)
maxv_helsinki = np.amax(my_Helsinki_PR) - np.amax(Helsinki_PR)
maxv_helsinki = (
"+0mm"
if maxv_helsinki == 0
else np.format_float_positional(
maxv_helsinki, precision=2, min_digits=2, sign=True
)
+ "mm"
)
maxv_era5_helsinki = np.amax(my_ERA5_Helsinki_PR) - np.amax(ERA5_Helsinki_PR)
maxv_era5_helsinki = (
"+0mm"
if maxv_era5_helsinki == 0
else np.format_float_positional(
maxv_era5_helsinki, precision=2, min_digits=2, sign=True
)
+ "mm"
)
maxv_era5 = np.mean(my_ERA5_PR.max("valid_time") - ERA5_PR.max("valid_time"))
maxv_era5 = (
"+0mm"
if maxv_era5 == 0
else np.format_float_positional(
maxv_era5, precision=2, min_digits=2, sign=True
)
)
if reference:
argmax_belem = pd.to_datetime(str(Time[np.argmax(Belem_PR)])).strftime(
"%d.%m %Hh"
)
argmax_era5_belem = pd.to_datetime(
str(Time[np.argmax(ERA5_Belem_PR)])
).strftime("%d.%m %Hh")
argmax_helsinki = pd.to_datetime(str(Time[np.argmax(Helsinki_PR)])).strftime(
"%d.%m %Hh"
)
argmax_era5_helsinki = pd.to_datetime(
str(Time[np.argmax(ERA5_Helsinki_PR)])
).strftime("%d.%m %Hh")
argmax_era5 = ""
else:
argmax_belem = (
Time[np.argmax(my_Belem_PR)] - Time[np.argmax(Belem_PR)]
) / np.timedelta64(1, "h")
argmax_belem = (
np.format_float_positional(argmax_belem, precision=0, sign=True, trim="-")
+ "h"
)
argmax_era5_belem = (
Time[np.argmax(my_ERA5_Belem_PR)] - Time[np.argmax(ERA5_Belem_PR)]
) / np.timedelta64(1, "h")
argmax_era5_belem = (
np.format_float_positional(
argmax_era5_belem, precision=0, sign=True, trim="-"
)
+ "h"
)
argmax_helsinki = (
Time[np.argmax(my_Helsinki_PR)] - Time[np.argmax(Helsinki_PR)]
) / np.timedelta64(1, "h")
argmax_helsinki = (
np.format_float_positional(
argmax_helsinki, precision=0, sign=True, trim="-"
)
+ "h"
)
argmax_era5_helsinki = (
Time[np.argmax(my_ERA5_Helsinki_PR)] - Time[np.argmax(ERA5_Helsinki_PR)]
) / np.timedelta64(1, "h")
argmax_era5_helsinki = (
np.format_float_positional(
argmax_era5_helsinki, precision=0, sign=True, trim="-"
)
+ "h"
)
argmax_era5 = (
np.mean(my_ERA5_PR.argmax("valid_time") - ERA5_PR.argmax("valid_time"))
* (Time[1] - Time[0])
/ np.timedelta64(1, "h")
)
argmax_era5 = (
np.format_float_positional(
argmax_era5, precision=1, min_digits=1, sign=True, trim="-"
)
+ "h"
)
if reference:
positive_belem = f"{np.sum(Belem_PR > 0)}" + "h"
positive_era5_belem = f"{np.sum(ERA5_Belem_PR > 0)}" + "h"
positive_helsinki = f"{np.sum(Helsinki_PR > 0)}" + "h"
positive_era5_helsinki = f"{np.sum(ERA5_Helsinki_PR > 0)}" + "h"
positive_era5 = (
np.format_float_positional(np.mean(ERA5_PR > 0), precision=2, min_digits=2)
+ "%"
)
else:
positive_belem = (
100
* (np.sum(my_Belem_PR > 0) - np.sum(Belem_PR > 0))
/ np.sum(Belem_PR > 0)
)
positive_belem = (
"+0h"
if positive_belem == 0
else np.format_float_positional(
positive_belem, precision=1, min_digits=1, sign=True
)
+ "%"
)
positive_era5_belem = (
100
* (np.sum(my_ERA5_Belem_PR > 0) - np.sum(ERA5_Belem_PR > 0))
/ np.sum(ERA5_Belem_PR > 0)
)
positive_era5_belem = (
"+0h"
if positive_era5_belem == 0
else np.format_float_positional(
positive_era5_belem, precision=1, min_digits=1, sign=True
)
+ "%"
)
positive_helsinki = (
100
* (np.sum(my_Helsinki_PR > 0) - np.sum(Helsinki_PR > 0))
/ np.sum(Helsinki_PR > 0)
)
positive_helsinki = (
"+0h"
if positive_helsinki == 0
else np.format_float_positional(
positive_helsinki, precision=1, min_digits=1, sign=True
)
+ "%"
)
positive_era5_helsinki = (
100
* (np.sum(my_ERA5_Helsinki_PR > 0) - np.sum(ERA5_Helsinki_PR > 0))
/ np.sum(ERA5_Helsinki_PR > 0)
)
positive_era5_helsinki = (
"+0h"
if positive_era5_helsinki == 0
else np.format_float_positional(
positive_era5_helsinki, precision=1, min_digits=1, sign=True
)
+ "%"
)
positive_era5 = (
np.format_float_positional(
np.mean(my_ERA5_PR > 0), precision=2, min_digits=2
)
+ "%"
)
if reference:
fpfn_belem = "0h"
fpfn_era5_belem = "0h"
fpfn_helsinki = "0h"
fpfn_era5_helsinki = "0h"
fpfn_era5 = "0h"
else:
fpfn_belem = 100 * np.mean((my_Belem_PR > 0) != (Belem_PR > 0))
fpfn_belem = (
"0h"
if fpfn_belem == 0
else np.format_float_positional(fpfn_belem, precision=1, min_digits=1) + "%"
)
fpfn_era5_belem = 100 * np.mean((my_ERA5_Belem_PR > 0) != (ERA5_Belem_PR > 0))
fpfn_era5_belem = (
"0h"
if fpfn_era5_belem == 0
else np.format_float_positional(fpfn_era5_belem, precision=1, min_digits=1)
+ "%"
)
fpfn_helsinki = 100 * np.mean((my_Helsinki_PR > 0) != (Helsinki_PR > 0))
fpfn_helsinki = (
"0h"
if fpfn_helsinki == 0
else np.format_float_positional(fpfn_helsinki, precision=1, min_digits=1)
+ "%"
)
fpfn_era5_helsinki = 100 * np.mean(
(my_ERA5_Helsinki_PR > 0) != (ERA5_Helsinki_PR > 0)
)
fpfn_era5_helsinki = (
"0h"
if fpfn_era5_helsinki == 0
else np.format_float_positional(
fpfn_era5_helsinki, precision=1, min_digits=1
)
+ "%"
)
fpfn_era5 = 100 * np.mean((my_ERA5_PR > 0) != (ERA5_PR > 0))
fpfn_era5 = (
"0h"
if fpfn_era5 == 0
else np.format_float_positional(fpfn_era5, precision=1, min_digits=1) + "%"
)
if reference:
negative_belem = f"{np.sum(Belem_PR < 0)}" + "h"
negative_era5_belem = f"{np.sum(ERA5_Belem_PR < 0)}" + "h"
negative_helsinki = f"{np.sum(Helsinki_PR < 0)}" + "h"
negative_era5_helsinki = f"{np.sum(ERA5_Helsinki_PR < 0)}" + "h"
negative_era5 = f"{int(np.sum(ERA5_PR < 0))}" + "h"
else:
negative_belem = 100 * np.mean(my_Belem_PR < 0)
negative_belem = (
"0h"
if negative_belem == 0
else np.format_float_positional(negative_belem, precision=1, min_digits=1)
+ "%"
)
negative_era5_belem = 100 * np.mean(my_ERA5_Belem_PR < 0)
negative_era5_belem = (
"0h"
if negative_era5_belem == 0
else np.format_float_positional(
negative_era5_belem, precision=1, min_digits=1
)
+ "%"
)
negative_helsinki = 100 * np.mean(my_Helsinki_PR < 0)
negative_helsinki = (
"0h"
if negative_helsinki == 0
else np.format_float_positional(
negative_helsinki, precision=1, min_digits=1
)
+ "%"
)
negative_era5_helsinki = 100 * np.mean(my_ERA5_Helsinki_PR < 0)
negative_era5_helsinki = (
"0h"
if negative_era5_helsinki == 0
else np.format_float_positional(
negative_era5_helsinki, precision=1, min_digits=1
)
+ "%"
)
negative_era5 = 100 * np.mean(my_ERA5_PR < 0)
negative_era5 = (
"0h"
if negative_era5 == 0
else np.format_float_positional(negative_era5, precision=1, min_digits=1)
+ "%"
)
if reference:
cr_helsinki = r"$\times$ 1"
cr_belem = r"$\times$ 1"
cr_era5 = r"$\times$ 1"
else:
cr_helsinki = rf"$\times$ {np.round(my_Helsinki_PR_cr, 2)}"
cr_belem = rf"$\times$ {np.round(my_Belem_PR_cr, 2)}"
cr_era5 = rf"$\times$ {np.round(my_ERA5_PR_cr, 2)}"
return pd.DataFrame(
{
"Compressor": [title[0], title[0], title[0], title[0], title[0]],
"Safeguarded": [title[1], title[1], title[1], title[1], title[1]],
"Corrections": [title[2], title[2], title[2], title[2], title[2]],
"Station": ["Helsinki", "Helsinki", "Belém", "Belém", "ERA5"],
"Source": ["obs", "ERA5", "obs", "ERA5", "avg"],
r"$L_{\infty}(\hat{PR})$": [
err_Helsinki_inf,
err_era5_Helsinki_inf,
err_Belem_inf,
err_era5_Belem_inf,
err_ERA5_inf,
],
r"$L_{2}(\hat{PR})$": [
err_Helsinki_2,
err_era5_Helsinki_2,
err_Belem_2,
err_era5_Belem_2,
err_ERA5_2,
],
"Integral": [
integral_helsinki,
integral_era5_helsinki,
integral_belem,
integral_era5_belem,
integral_era5,
],
"max(PR)": [
maxv_helsinki,
maxv_era5_helsinki,
maxv_belem,
maxv_era5_belem,
maxv_era5,
],
"argmax(PR)": [
argmax_helsinki,
argmax_era5_helsinki,
argmax_belem,
argmax_era5_belem,
argmax_era5,
],
"PR > 0": [
positive_helsinki,
positive_era5_helsinki,
positive_belem,
positive_era5_belem,
positive_era5,
],
"FP + FN": [
fpfn_helsinki,
fpfn_era5_helsinki,
fpfn_belem,
fpfn_era5_belem,
fpfn_era5,
],
"PR < 0": [
negative_helsinki,
negative_era5_helsinki,
negative_belem,
negative_era5_belem,
negative_era5,
],
"V": [
err_Helsinki_v,
err_era5_Helsinki_v,
err_Belem_v,
err_era5_Belem_v,
err_ERA5_v,
],
"C": [
corr_Helsinki,
corr_era5_Helsinki,
corr_Belem,
corr_era5_Belem,
corr_era5,
],
"CR": [cr_helsinki, "", cr_belem, "", cr_era5],
}
)
Copied!
import observe
observations = []
import observe
observations = []
Lossless compression¶
We first compress the data losslessly with ZStandard at level 22, which gives maximum compression, to provide a baseline.
Copied!
def encode_decode(codec):
with observe.observe(codec, observations):
ERA5_PR_codec_enc = codec.encode(ERA5_PR.values)
ERA5_PR_codec = ERA5_PR.copy(data=codec.decode(ERA5_PR_codec_enc))
ERA5_PR_codec_cr = ERA5_PR.nbytes / np.asarray(ERA5_PR_codec_enc).nbytes
Belem_PR_codec = Belem_PR.copy(deep=True)
with observe.observe(codec, observations):
Belem_PR_codec_enc = codec.encode(Belem_PR_codec.values)
Belem_PR_codec.values[:] = codec.decode(Belem_PR_codec_enc)
Belem_PR_codec_cr = Belem_PR.nbytes / np.asarray(Belem_PR_codec_enc).nbytes
Helsinki_PR_codec = Helsinki_PR.copy(deep=True)
with observe.observe(codec, observations):
Helsinki_PR_codec_enc = codec.encode(Helsinki_PR_codec.values)
Helsinki_PR_codec.values[:] = codec.decode(Helsinki_PR_codec_enc)
Helsinki_PR_codec_cr = Helsinki_PR.nbytes / np.asarray(Helsinki_PR_codec_enc).nbytes
return (
ERA5_PR_codec,
ERA5_PR_codec_cr,
Belem_PR_codec,
Belem_PR_codec_cr,
Helsinki_PR_codec,
Helsinki_PR_codec_cr,
)
def encode_decode(codec):
with observe.observe(codec, observations):
ERA5_PR_codec_enc = codec.encode(ERA5_PR.values)
ERA5_PR_codec = ERA5_PR.copy(data=codec.decode(ERA5_PR_codec_enc))
ERA5_PR_codec_cr = ERA5_PR.nbytes / np.asarray(ERA5_PR_codec_enc).nbytes
Belem_PR_codec = Belem_PR.copy(deep=True)
with observe.observe(codec, observations):
Belem_PR_codec_enc = codec.encode(Belem_PR_codec.values)
Belem_PR_codec.values[:] = codec.decode(Belem_PR_codec_enc)
Belem_PR_codec_cr = Belem_PR.nbytes / np.asarray(Belem_PR_codec_enc).nbytes
Helsinki_PR_codec = Helsinki_PR.copy(deep=True)
with observe.observe(codec, observations):
Helsinki_PR_codec_enc = codec.encode(Helsinki_PR_codec.values)
Helsinki_PR_codec.values[:] = codec.decode(Helsinki_PR_codec_enc)
Helsinki_PR_codec_cr = Helsinki_PR.nbytes / np.asarray(Helsinki_PR_codec_enc).nbytes
return (
ERA5_PR_codec,
ERA5_PR_codec_cr,
Belem_PR_codec,
Belem_PR_codec_cr,
Helsinki_PR_codec,
Helsinki_PR_codec_cr,
)
Copied!
from numcodecs_wasm_zstd import Zstd
zstd = Zstd(level=22)
(
ERA5_PR_zstd,
ERA5_PR_zstd_cr,
Belem_PR_zstd,
Belem_PR_zstd_cr,
Helsinki_PR_zstd,
Helsinki_PR_zstd_cr,
) = encode_decode(zstd)
from numcodecs_wasm_zstd import Zstd
zstd = Zstd(level=22)
(
ERA5_PR_zstd,
ERA5_PR_zstd_cr,
Belem_PR_zstd,
Belem_PR_zstd_cr,
Helsinki_PR_zstd,
Helsinki_PR_zstd_cr,
) = encode_decode(zstd)
Compressing precipitation with lossy compressors¶
We configure each compressor with an absolute error bound of 0.1mm. The observation time series are compressed independently, for ERA5 we compress the entire three day dataset before extracting the observation-space time series.
Copied!
from numcodecs_wasm_sperr import Sperr
from numcodecs_wasm_sz3 import Sz3
from numcodecs_wasm_zfp import Zfp
from numcodecs_zero import ZeroCodec
from numcodecs_wasm_sperr import Sperr
from numcodecs_wasm_sz3 import Sz3
from numcodecs_wasm_zfp import Zfp
from numcodecs_zero import ZeroCodec
Copied!
eb_abs = 0.1
eb_abs = 0.1
Copied!
zfp = Zfp(mode="fixed-accuracy", tolerance=eb_abs)
(
ERA5_PR_zfp,
ERA5_PR_zfp_cr,
Belem_PR_zfp,
Belem_PR_zfp_cr,
Helsinki_PR_zfp,
Helsinki_PR_zfp_cr,
) = encode_decode(zfp)
zfp = Zfp(mode="fixed-accuracy", tolerance=eb_abs)
(
ERA5_PR_zfp,
ERA5_PR_zfp_cr,
Belem_PR_zfp,
Belem_PR_zfp_cr,
Helsinki_PR_zfp,
Helsinki_PR_zfp_cr,
) = encode_decode(zfp)
Copied!
sz3 = Sz3(eb_mode="abs", eb_abs=eb_abs)
(
ERA5_PR_sz3,
ERA5_PR_sz3_cr,
Belem_PR_sz3,
Belem_PR_sz3_cr,
Helsinki_PR_sz3,
Helsinki_PR_sz3_cr,
) = encode_decode(sz3)
sz3 = Sz3(eb_mode="abs", eb_abs=eb_abs)
(
ERA5_PR_sz3,
ERA5_PR_sz3_cr,
Belem_PR_sz3,
Belem_PR_sz3_cr,
Helsinki_PR_sz3,
Helsinki_PR_sz3_cr,
) = encode_decode(sz3)
Copied!
sperr = Sperr(mode="pwe", pwe=eb_abs)
(
ERA5_PR_sperr,
ERA5_PR_sperr_cr,
Belem_PR_sperr,
Belem_PR_sperr_cr,
Helsinki_PR_sperr,
Helsinki_PR_sperr_cr,
) = encode_decode(sperr)
sperr = Sperr(mode="pwe", pwe=eb_abs)
(
ERA5_PR_sperr,
ERA5_PR_sperr_cr,
Belem_PR_sperr,
Belem_PR_sperr_cr,
Helsinki_PR_sperr,
Helsinki_PR_sperr_cr,
) = encode_decode(sperr)
Copied!
zero = ZeroCodec()
(
ERA5_PR_zero,
_,
Belem_PR_zero,
_,
Helsinki_PR_zero,
_,
) = encode_decode(zero)
zero = ZeroCodec()
(
ERA5_PR_zero,
_,
Belem_PR_zero,
_,
Helsinki_PR_zero,
_,
) = encode_decode(zero)
Compressing precipitation with the safeguarded lossy compressors¶
We configure the safeguards with an absolute error bound of 0.1mm, and to preserve:
- zero and positive values by preserving the sign of the data
- the maximum by preserving the sign relative to the maximum (the maximum retains its exact values, other values cannot exceed it)
For compressing the observation time series, the maximum is the global maximum, which the numcodecs-safeguards frontend exposes as the built-in $x_max constant. For compressing ERA5, we want to preserve the observed maxima at Helsinki and Belém, so we compute them beforehand and preserve the two maxima separately.
Copied!
from numcodecs_safeguards import SafeguardedCodec
from numcodecs_safeguards import SafeguardedCodec
Copied!
ERA5_PR_sg = dict()
ERA5_PR_sg_cr = dict()
for codec in [zero, zfp, sz3, sperr]:
sg_era5 = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs),
dict(kind="sign"),
dict(kind="sign", offset=float(np.amax(ERA5_Belem_PR))),
dict(kind="sign", offset=float(np.amax(ERA5_Helsinki_PR))),
],
)
with observe.observe(sg_era5, observations):
ERA5_PR_sg_enc = sg_era5.encode(ERA5_PR.values)
ERA5_PR_sg[codec.codec_id] = ERA5_PR.copy(data=sg_era5.decode(ERA5_PR_sg_enc))
ERA5_PR_sg_cr[codec.codec_id] = ERA5_PR.nbytes / np.asarray(ERA5_PR_sg_enc).nbytes
ERA5_PR_sg = dict()
ERA5_PR_sg_cr = dict()
for codec in [zero, zfp, sz3, sperr]:
sg_era5 = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs),
dict(kind="sign"),
dict(kind="sign", offset=float(np.amax(ERA5_Belem_PR))),
dict(kind="sign", offset=float(np.amax(ERA5_Helsinki_PR))),
],
)
with observe.observe(sg_era5, observations):
ERA5_PR_sg_enc = sg_era5.encode(ERA5_PR.values)
ERA5_PR_sg[codec.codec_id] = ERA5_PR.copy(data=sg_era5.decode(ERA5_PR_sg_enc))
ERA5_PR_sg_cr[codec.codec_id] = ERA5_PR.nbytes / np.asarray(ERA5_PR_sg_enc).nbytes
Copied!
Belem_PR_sg = dict()
Belem_PR_sg_cr = dict()
Helsinki_PR_sg = dict()
Helsinki_PR_sg_cr = dict()
for codec in [ZeroCodec(), zfp, sz3, sperr]:
sg_obs = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs),
dict(kind="sign"),
# we only want to preserve *the* global maximum
dict(kind="sign", offset="$x_max"),
],
)
Belem_PR_sg[codec.codec_id] = Belem_PR.copy(deep=True)
with observe.observe(sg_obs, observations):
Belem_PR_sg_enc = sg_obs.encode(Belem_PR_sg[codec.codec_id].values)
Belem_PR_sg[codec.codec_id].values[:] = sg_obs.decode(Belem_PR_sg_enc)
Belem_PR_sg_cr[codec.codec_id] = (
Belem_PR.nbytes / np.asarray(Belem_PR_sg_enc).nbytes
)
Helsinki_PR_sg[codec.codec_id] = Helsinki_PR.copy(deep=True)
with observe.observe(sg_obs, observations):
Helsinki_PR_sg_enc = sg_obs.encode(Helsinki_PR_sg[codec.codec_id].values)
Helsinki_PR_sg[codec.codec_id].values[:] = sg_obs.decode(Helsinki_PR_sg_enc)
Helsinki_PR_sg_cr[codec.codec_id] = (
Helsinki_PR.nbytes / np.asarray(Helsinki_PR_sg_enc).nbytes
)
Belem_PR_sg = dict()
Belem_PR_sg_cr = dict()
Helsinki_PR_sg = dict()
Helsinki_PR_sg_cr = dict()
for codec in [ZeroCodec(), zfp, sz3, sperr]:
sg_obs = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs),
dict(kind="sign"),
# we only want to preserve *the* global maximum
dict(kind="sign", offset="$x_max"),
],
)
Belem_PR_sg[codec.codec_id] = Belem_PR.copy(deep=True)
with observe.observe(sg_obs, observations):
Belem_PR_sg_enc = sg_obs.encode(Belem_PR_sg[codec.codec_id].values)
Belem_PR_sg[codec.codec_id].values[:] = sg_obs.decode(Belem_PR_sg_enc)
Belem_PR_sg_cr[codec.codec_id] = (
Belem_PR.nbytes / np.asarray(Belem_PR_sg_enc).nbytes
)
Helsinki_PR_sg[codec.codec_id] = Helsinki_PR.copy(deep=True)
with observe.observe(sg_obs, observations):
Helsinki_PR_sg_enc = sg_obs.encode(Helsinki_PR_sg[codec.codec_id].values)
Helsinki_PR_sg[codec.codec_id].values[:] = sg_obs.decode(Helsinki_PR_sg_enc)
Helsinki_PR_sg_cr[codec.codec_id] = (
Helsinki_PR.nbytes / np.asarray(Helsinki_PR_sg_enc).nbytes
)
Copied!
ERA5_PR_sg_lossless = dict()
ERA5_PR_sg_lossless_cr = dict()
for codec in [zero, zfp, sz3, sperr]:
sg_era5 = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs),
dict(kind="sign"),
dict(kind="sign", offset=float(np.amax(ERA5_Belem_PR))),
dict(kind="sign", offset=float(np.amax(ERA5_Helsinki_PR))),
],
# produce lossless corrections and refine them with iteration
compute=dict(unstable_iterative=True, unstable_lossless_corrections=True),
)
with observe.observe(sg_era5, observations):
ERA5_PR_sg_lossless_enc = sg_era5.encode(ERA5_PR.values)
ERA5_PR_sg_lossless[codec.codec_id] = ERA5_PR.copy(
data=sg_era5.decode(ERA5_PR_sg_lossless_enc)
)
ERA5_PR_sg_lossless_cr[codec.codec_id] = (
ERA5_PR.nbytes / np.asarray(ERA5_PR_sg_lossless_enc).nbytes
)
ERA5_PR_sg_lossless = dict()
ERA5_PR_sg_lossless_cr = dict()
for codec in [zero, zfp, sz3, sperr]:
sg_era5 = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs),
dict(kind="sign"),
dict(kind="sign", offset=float(np.amax(ERA5_Belem_PR))),
dict(kind="sign", offset=float(np.amax(ERA5_Helsinki_PR))),
],
# produce lossless corrections and refine them with iteration
compute=dict(unstable_iterative=True, unstable_lossless_corrections=True),
)
with observe.observe(sg_era5, observations):
ERA5_PR_sg_lossless_enc = sg_era5.encode(ERA5_PR.values)
ERA5_PR_sg_lossless[codec.codec_id] = ERA5_PR.copy(
data=sg_era5.decode(ERA5_PR_sg_lossless_enc)
)
ERA5_PR_sg_lossless_cr[codec.codec_id] = (
ERA5_PR.nbytes / np.asarray(ERA5_PR_sg_lossless_enc).nbytes
)
Copied!
Belem_PR_sg_lossless = dict()
Belem_PR_sg_lossless_cr = dict()
Helsinki_PR_sg_lossless = dict()
Helsinki_PR_sg_lossless_cr = dict()
for codec in [ZeroCodec(), zfp, sz3, sperr]:
sg_obs = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs),
dict(kind="sign"),
# we only want to preserve *the* global maximum
dict(kind="sign", offset="$x_max"),
],
# produce lossless corrections and refine them with iteration
compute=dict(unstable_iterative=True, unstable_lossless_corrections=True),
)
Belem_PR_sg_lossless[codec.codec_id] = Belem_PR.copy(deep=True)
with observe.observe(sg_obs, observations):
Belem_PR_sg_lossless_enc = sg_obs.encode(Belem_PR_sg[codec.codec_id].values)
Belem_PR_sg_lossless[codec.codec_id].values[:] = sg_obs.decode(
Belem_PR_sg_lossless_enc
)
Belem_PR_sg_lossless_cr[codec.codec_id] = (
Belem_PR.nbytes / np.asarray(Belem_PR_sg_lossless_enc).nbytes
)
Helsinki_PR_sg_lossless[codec.codec_id] = Helsinki_PR.copy(deep=True)
with observe.observe(sg_obs, observations):
Helsinki_PR_sg_lossless_enc = sg_obs.encode(
Helsinki_PR_sg[codec.codec_id].values
)
Helsinki_PR_sg_lossless[codec.codec_id].values[:] = sg_obs.decode(
Helsinki_PR_sg_lossless_enc
)
Helsinki_PR_sg_lossless_cr[codec.codec_id] = (
Helsinki_PR.nbytes / np.asarray(Helsinki_PR_sg_lossless_enc).nbytes
)
Belem_PR_sg_lossless = dict()
Belem_PR_sg_lossless_cr = dict()
Helsinki_PR_sg_lossless = dict()
Helsinki_PR_sg_lossless_cr = dict()
for codec in [ZeroCodec(), zfp, sz3, sperr]:
sg_obs = SafeguardedCodec(
codec=codec,
safeguards=[
dict(kind="eb", type="abs", eb=eb_abs),
dict(kind="sign"),
# we only want to preserve *the* global maximum
dict(kind="sign", offset="$x_max"),
],
# produce lossless corrections and refine them with iteration
compute=dict(unstable_iterative=True, unstable_lossless_corrections=True),
)
Belem_PR_sg_lossless[codec.codec_id] = Belem_PR.copy(deep=True)
with observe.observe(sg_obs, observations):
Belem_PR_sg_lossless_enc = sg_obs.encode(Belem_PR_sg[codec.codec_id].values)
Belem_PR_sg_lossless[codec.codec_id].values[:] = sg_obs.decode(
Belem_PR_sg_lossless_enc
)
Belem_PR_sg_lossless_cr[codec.codec_id] = (
Belem_PR.nbytes / np.asarray(Belem_PR_sg_lossless_enc).nbytes
)
Helsinki_PR_sg_lossless[codec.codec_id] = Helsinki_PR.copy(deep=True)
with observe.observe(sg_obs, observations):
Helsinki_PR_sg_lossless_enc = sg_obs.encode(
Helsinki_PR_sg[codec.codec_id].values
)
Helsinki_PR_sg_lossless[codec.codec_id].values[:] = sg_obs.decode(
Helsinki_PR_sg_lossless_enc
)
Helsinki_PR_sg_lossless_cr[codec.codec_id] = (
Helsinki_PR.nbytes / np.asarray(Helsinki_PR_sg_lossless_enc).nbytes
)
Compressing precipitation using OptZConfig¶
We configure OptZConfig with a custom safety violations metric, implemented in Python, that computes the percentage of violations $V$. We then maximise the score
$$ \textrm{score} = \begin{cases} -\textrm{V} \quad &\text{if } \textrm{V} > 0 \\ \textrm{CR} \quad &\text{otherwise} \end{cases} $$
using the FRAZ search algorithm with 25 iterations, where CR is the achieved compression ratio. Since FRAZ seems to struggle with finding sufficient absolute error bounds spread across several orders of magnitude, we search for bounds in logarithmic space by wrapping each codec in an Exponential<CODEC> meta-codec.
Copied!
import numcodecs
class SafetyViolationsMetric(numcodecs.abc.Codec):
codec_id = "safety-violations-metric"
def __init__(self):
self._data = None
def encode(self, buf):
# store the original data for later
self._data = np.array(buf, copy=True)
# return no metric
return np.empty(0, dtype=np.float64)
def decode(self, buf, out=None):
# compute the violations
violations = np.mean(
~(
(np.abs(buf - self._data) <= eb_abs)
& (np.sign(buf) == np.sign(self._data))
)
)
self._data = None
# return the violations score metric
return numcodecs.compat.ndarray_copy(np.float64(violations), out)
numcodecs.registry.register_codec(SafetyViolationsMetric)
import numcodecs
class SafetyViolationsMetric(numcodecs.abc.Codec):
codec_id = "safety-violations-metric"
def __init__(self):
self._data = None
def encode(self, buf):
# store the original data for later
self._data = np.array(buf, copy=True)
# return no metric
return np.empty(0, dtype=np.float64)
def decode(self, buf, out=None):
# compute the violations
violations = np.mean(
~(
(np.abs(buf - self._data) <= eb_abs)
& (np.sign(buf) == np.sign(self._data))
)
)
self._data = None
# return the violations score metric
return numcodecs.compat.ndarray_copy(np.float64(violations), out)
numcodecs.registry.register_codec(SafetyViolationsMetric)
Copied!
class ExponentialZfp(Zfp):
codec_id = "e-zfp.rs"
def __new__(cls, tolerance: float, **kwargs):
codec = super().__new__(cls, tolerance=np.exp(tolerance), **kwargs)
codec._tolerance = tolerance
return codec
def get_config(self):
return {
**super().get_config(),
"id": type(self).codec_id,
"tolerance": self._tolerance,
}
class ExponentialSz3(Sz3):
codec_id = "e-sz3.rs"
def __new__(cls, eb_abs: float, **kwargs):
codec = super().__new__(cls, eb_abs=np.exp(eb_abs), **kwargs)
codec._eb_abs = eb_abs
return codec
def get_config(self):
return {
**super().get_config(),
"id": type(self).codec_id,
"eb_abs": self._eb_abs,
}
class ExponentialSperr(Sperr):
codec_id = "e-sperr.rs"
def __new__(cls, pwe: float, **kwargs):
codec = super().__new__(cls, pwe=np.exp(pwe), **kwargs)
codec._pwe = pwe
return codec
def get_config(self):
return {**super().get_config(), "id": type(self).codec_id, "pwe": self._pwe}
numcodecs.registry.register_codec(ExponentialZfp)
numcodecs.registry.register_codec(ExponentialSz3)
numcodecs.registry.register_codec(ExponentialSperr)
class ExponentialZfp(Zfp):
codec_id = "e-zfp.rs"
def __new__(cls, tolerance: float, **kwargs):
codec = super().__new__(cls, tolerance=np.exp(tolerance), **kwargs)
codec._tolerance = tolerance
return codec
def get_config(self):
return {
**super().get_config(),
"id": type(self).codec_id,
"tolerance": self._tolerance,
}
class ExponentialSz3(Sz3):
codec_id = "e-sz3.rs"
def __new__(cls, eb_abs: float, **kwargs):
codec = super().__new__(cls, eb_abs=np.exp(eb_abs), **kwargs)
codec._eb_abs = eb_abs
return codec
def get_config(self):
return {
**super().get_config(),
"id": type(self).codec_id,
"eb_abs": self._eb_abs,
}
class ExponentialSperr(Sperr):
codec_id = "e-sperr.rs"
def __new__(cls, pwe: float, **kwargs):
codec = super().__new__(cls, pwe=np.exp(pwe), **kwargs)
codec._pwe = pwe
return codec
def get_config(self):
return {**super().get_config(), "id": type(self).codec_id, "pwe": self._pwe}
numcodecs.registry.register_codec(ExponentialZfp)
numcodecs.registry.register_codec(ExponentialSz3)
numcodecs.registry.register_codec(ExponentialSperr)
Copied!
from numcodecs_wasm_pressio import Pressio
ERA5_PR_optzconfig = dict()
ERA5_PR_optzconfig_cr = dict()
for codec, parameter, lower_bound in [
(zfp, "tolerance", 1e-20), # initial guess
(sz3, "eb_abs", 1e-15), # initial guess
(sperr, "pwe", 1e-15), # initial guess
]:
optzconfig_era5 = Pressio(
compressor_id="opt",
compressor_config={
"opt:output": ["composite:score"],
"opt:inputs": [f"numcodecs.rs:{parameter}"],
"opt:lower_bound": np.log(lower_bound),
"opt:upper_bound": np.log(eb_abs),
"opt:max_iterations": 25,
"opt:objective_mode_name": "max",
},
early_config={
"opt:compressor": "pressio",
"pressio:compressor": "numcodecs.rs",
**{
f"numcodecs.rs:{k}": f"e-{v}" if k == "id" else v
for k, v in codec.get_config().items()
},
"opt:search": "fraz",
"pressio:metric": "composite",
"composite:plugins": ["size", "numcodecs.rs-metric"],
"composite:scripts": [
"""
violations = metrics["numcodecs.rs-metric:decompression"]
if violations > 0 then
return "score", -violations
else
return "score", metrics["size:compression_ratio"]
end
"""
],
"numcodecs.rs-metric:id": "safety-violations-metric",
},
)
with observe.observe(optzconfig_era5, observations):
ERA5_PR_optzconfig_enc = optzconfig_era5.encode(ERA5_PR.values)
ERA5_PR_optzconfig[codec.codec_id] = ERA5_PR.copy(
data=optzconfig_era5.decode(ERA5_PR_optzconfig_enc)
)
ERA5_PR_optzconfig_cr[codec.codec_id] = (
ERA5_PR.nbytes / np.asarray(ERA5_PR_optzconfig_enc).nbytes
)
from numcodecs_wasm_pressio import Pressio
ERA5_PR_optzconfig = dict()
ERA5_PR_optzconfig_cr = dict()
for codec, parameter, lower_bound in [
(zfp, "tolerance", 1e-20), # initial guess
(sz3, "eb_abs", 1e-15), # initial guess
(sperr, "pwe", 1e-15), # initial guess
]:
optzconfig_era5 = Pressio(
compressor_id="opt",
compressor_config={
"opt:output": ["composite:score"],
"opt:inputs": [f"numcodecs.rs:{parameter}"],
"opt:lower_bound": np.log(lower_bound),
"opt:upper_bound": np.log(eb_abs),
"opt:max_iterations": 25,
"opt:objective_mode_name": "max",
},
early_config={
"opt:compressor": "pressio",
"pressio:compressor": "numcodecs.rs",
**{
f"numcodecs.rs:{k}": f"e-{v}" if k == "id" else v
for k, v in codec.get_config().items()
},
"opt:search": "fraz",
"pressio:metric": "composite",
"composite:plugins": ["size", "numcodecs.rs-metric"],
"composite:scripts": [
"""
violations = metrics["numcodecs.rs-metric:decompression"]
if violations > 0 then
return "score", -violations
else
return "score", metrics["size:compression_ratio"]
end
"""
],
"numcodecs.rs-metric:id": "safety-violations-metric",
},
)
with observe.observe(optzconfig_era5, observations):
ERA5_PR_optzconfig_enc = optzconfig_era5.encode(ERA5_PR.values)
ERA5_PR_optzconfig[codec.codec_id] = ERA5_PR.copy(
data=optzconfig_era5.decode(ERA5_PR_optzconfig_enc)
)
ERA5_PR_optzconfig_cr[codec.codec_id] = (
ERA5_PR.nbytes / np.asarray(ERA5_PR_optzconfig_enc).nbytes
)
rank={0,1,} iter={0} input={-24.1771,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={1} input={-35.8462,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={2} input={-12.6864,} output={-0.163744,} objective={-0.163744}
rank={0,1,} iter={3} input={-36.7947,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={4} input={-5.28951,} output={-0.144939,} objective={-0.144939}
rank={0,1,} iter={5} input={-43.1005,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={6} input={-6.65975,} output={-0.162324,} objective={-0.162324}
rank={0,1,} iter={7} input={-5.79875,} output={-0.155607,} objective={-0.155607}
rank={0,1,} iter={8} input={-24.5978,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={9} input={-39.1015,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={10} input={-10.546,} output={-0.164972,} objective={-0.164972}
rank={0,1,} iter={11} input={-24.3295,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={12} input={-39.6063,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={13} input={-9.31911,} output={-0.16487,} objective={-0.16487}
rank={0,1,} iter={14} input={-42.3905,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={15} input={-36.1439,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={16} input={-30.1969,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={17} input={-2.89648,} output={-0.125638,} objective={-0.125638}
rank={0,1,} iter={18} input={-23.4451,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={19} input={-6.06246,} output={-0.155607,} objective={-0.155607}
rank={0,1,} iter={20} input={-45.698,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={21} input={-39.0272,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={22} input={-10.4072,} output={-0.164972,} objective={-0.164972}
rank={0,1,} iter={23} input={-42.5172,} output={-4.27005e-05,} objective={-4.27005e-05}
rank={0,1,} iter={24} input={-3.60391,} output={-0.12807,} objective={-0.12807}
final_iter={25} inputs={-24.1771,} output={-4.27005e-05,}
rank={0,1,} iter={0} input={-18.4207,} output={-0.128425,} objective={-0.128425}
rank={0,1,} iter={1} input={-27.019,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={2} input={-9.95384,} output={-0.318886,} objective={-0.318886}
rank={0,1,} iter={3} input={-34.5388,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={4} input={-30.7838,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={5} input={-28.9018,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={6} input={-32.6581,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={7} input={-27.9601,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={8} input={-33.6007,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={9} input={-29.8388,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={10} input={-31.7232,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={11} input={-33.1279,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={12} input={-29.3705,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={13} input={-27.4901,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={14} input={-28.4303,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={15} input={-31.2541,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={16} input={-30.3114,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={17} input={-34.0634,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={18} input={-32.1941,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={19} input={-28.6658,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={20} input={-34.3016,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={21} input={-30.5471,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={22} input={-30.0768,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={23} input={-33.3646,} output={3.84178,} objective={3.84178}
rank={0,1,} iter={24} input={-27.2549,} output={3.84178,} objective={3.84178}
final_iter={25} inputs={-30.7838,} output={3.84178,}
rank={0,1,} iter={0} input={-18.4207,} output={-0.34835,} objective={-0.34835}
rank={0,1,} iter={1} input={-27.019,} output={-0.348347,} objective={-0.348347}
rank={0,1,} iter={2} input={-9.95384,} output={-0.34835,} objective={-0.34835}
rank={0,1,} iter={3} input={-27.7178,} output={-0.348344,} objective={-0.348344}
rank={0,1,} iter={4} input={-4.50347,} output={-0.379522,} objective={-0.379522}
rank={0,1,} iter={5} input={-32.3642,} output={-0.347709,} objective={-0.347709}
rank={0,1,} iter={6} input={-5.51313,} output={-0.363145,} objective={-0.363145}
rank={0,1,} iter={7} input={-4.87871,} output={-0.372818,} objective={-0.372818}
rank={0,1,} iter={8} input={-18.7306,} output={-0.34835,} objective={-0.34835}
rank={0,1,} iter={9} input={-29.4176,} output={-0.348316,} objective={-0.348316}
rank={0,1,} iter={10} input={-8.37672,} output={-0.34835,} objective={-0.34835}
rank={0,1,} iter={11} input={-18.5329,} output={-0.34835,} objective={-0.34835}
rank={0,1,} iter={12} input={-29.7896,} output={-0.348301,} objective={-0.348301}
rank={0,1,} iter={13} input={-7.47266,} output={-0.348611,} objective={-0.348611}
rank={0,1,} iter={14} input={-31.841,} output={-0.347973,} objective={-0.347973}
rank={0,1,} iter={15} input={-27.2383,} output={-0.348347,} objective={-0.348347}
rank={0,1,} iter={16} input={-22.8563,} output={-0.34835,} objective={-0.34835}
rank={0,1,} iter={17} input={-2.74019,} output={-0.41689,} objective={-0.41689}
rank={0,1,} iter={18} input={-17.8813,} output={-0.34835,} objective={-0.34835}
rank={0,1,} iter={19} input={-5.07302,} output={-0.369629,} objective={-0.369629}
rank={0,1,} iter={20} input={-34.2781,} output={-0.344363,} objective={-0.344363}
rank={0,1,} iter={21} input={-29.3628,} output={-0.348319,} objective={-0.348319}
rank={0,1,} iter={22} input={-8.2744,} output={-0.34835,} objective={-0.34835}
rank={0,1,} iter={23} input={-31.9344,} output={-0.347937,} objective={-0.347937}
rank={0,1,} iter={24} input={-3.26146,} output={-0.40547,} objective={-0.40547}
final_iter={25} inputs={-18.4207,} output={-0.34835,}
Copied!
Belem_PR_optzconfig = dict()
Belem_PR_optzconfig_cr = dict()
Helsinki_PR_optzconfig = dict()
Helsinki_PR_optzconfig_cr = dict()
for codec, parameter, lower_bound in [
(zfp, "tolerance", 1e-20), # tiny bound
(sz3, "eb_abs", 1e-6), # decent guess
(sperr, "pwe", 1e-16), # crash for lower error bounds
]:
optzconfig_obs = Pressio(
compressor_id="opt",
compressor_config={
"opt:output": ["composite:score"],
"opt:inputs": [f"numcodecs.rs:{parameter}"],
"opt:lower_bound": np.log(lower_bound),
"opt:upper_bound": np.log(eb_abs),
"opt:max_iterations": 25,
"opt:objective_mode_name": "max",
},
early_config={
"opt:compressor": "pressio",
"pressio:compressor": "numcodecs.rs",
**{
f"numcodecs.rs:{k}": f"e-{v}" if k == "id" else v
for k, v in codec.get_config().items()
},
"opt:search": "fraz",
"pressio:metric": "composite",
"composite:plugins": ["size", "numcodecs.rs-metric"],
"composite:scripts": [
"""
violations = metrics["numcodecs.rs-metric:decompression"]
if violations > 0 then
return "score", -violations
else
return "score", metrics["size:compression_ratio"]
end
"""
],
"numcodecs.rs-metric:id": "safety-violations-metric",
},
)
Belem_PR_optzconfig[codec.codec_id] = Belem_PR.copy(deep=True)
with observe.observe(optzconfig_obs, observations):
Belem_PR_optzconfig_enc = optzconfig_obs.encode(
Belem_PR_optzconfig[codec.codec_id].values
)
Belem_PR_optzconfig[codec.codec_id].values[:] = optzconfig_obs.decode(
Belem_PR_optzconfig_enc
)
Belem_PR_optzconfig_cr[codec.codec_id] = (
Belem_PR.nbytes / np.asarray(Belem_PR_optzconfig_enc).nbytes
)
Helsinki_PR_optzconfig[codec.codec_id] = Helsinki_PR.copy(deep=True)
with observe.observe(optzconfig_obs, observations):
Helsinki_PR_optzconfig_enc = optzconfig_obs.encode(
Helsinki_PR_optzconfig[codec.codec_id].values
)
Helsinki_PR_optzconfig[codec.codec_id].values[:] = optzconfig_obs.decode(
Helsinki_PR_optzconfig_enc
)
Helsinki_PR_optzconfig_cr[codec.codec_id] = (
Helsinki_PR.nbytes / np.asarray(Helsinki_PR_optzconfig_enc).nbytes
)
Belem_PR_optzconfig = dict()
Belem_PR_optzconfig_cr = dict()
Helsinki_PR_optzconfig = dict()
Helsinki_PR_optzconfig_cr = dict()
for codec, parameter, lower_bound in [
(zfp, "tolerance", 1e-20), # tiny bound
(sz3, "eb_abs", 1e-6), # decent guess
(sperr, "pwe", 1e-16), # crash for lower error bounds
]:
optzconfig_obs = Pressio(
compressor_id="opt",
compressor_config={
"opt:output": ["composite:score"],
"opt:inputs": [f"numcodecs.rs:{parameter}"],
"opt:lower_bound": np.log(lower_bound),
"opt:upper_bound": np.log(eb_abs),
"opt:max_iterations": 25,
"opt:objective_mode_name": "max",
},
early_config={
"opt:compressor": "pressio",
"pressio:compressor": "numcodecs.rs",
**{
f"numcodecs.rs:{k}": f"e-{v}" if k == "id" else v
for k, v in codec.get_config().items()
},
"opt:search": "fraz",
"pressio:metric": "composite",
"composite:plugins": ["size", "numcodecs.rs-metric"],
"composite:scripts": [
"""
violations = metrics["numcodecs.rs-metric:decompression"]
if violations > 0 then
return "score", -violations
else
return "score", metrics["size:compression_ratio"]
end
"""
],
"numcodecs.rs-metric:id": "safety-violations-metric",
},
)
Belem_PR_optzconfig[codec.codec_id] = Belem_PR.copy(deep=True)
with observe.observe(optzconfig_obs, observations):
Belem_PR_optzconfig_enc = optzconfig_obs.encode(
Belem_PR_optzconfig[codec.codec_id].values
)
Belem_PR_optzconfig[codec.codec_id].values[:] = optzconfig_obs.decode(
Belem_PR_optzconfig_enc
)
Belem_PR_optzconfig_cr[codec.codec_id] = (
Belem_PR.nbytes / np.asarray(Belem_PR_optzconfig_enc).nbytes
)
Helsinki_PR_optzconfig[codec.codec_id] = Helsinki_PR.copy(deep=True)
with observe.observe(optzconfig_obs, observations):
Helsinki_PR_optzconfig_enc = optzconfig_obs.encode(
Helsinki_PR_optzconfig[codec.codec_id].values
)
Helsinki_PR_optzconfig[codec.codec_id].values[:] = optzconfig_obs.decode(
Helsinki_PR_optzconfig_enc
)
Helsinki_PR_optzconfig_cr[codec.codec_id] = (
Helsinki_PR.nbytes / np.asarray(Helsinki_PR_optzconfig_enc).nbytes
)
rank={0,1,} iter={0} input={-24.1771,} output={-0.125,} objective={-0.125}
rank={0,1,} iter={1} input={-35.8462,} output={-0.125,} objective={-0.125}
rank={0,1,} iter={2} input={-12.6864,} output={-0.125,} objective={-0.125}
rank={0,1,} iter={3} input={-36.7947,} output={-0.0694444,} objective={-0.0694444}
rank={0,1,} iter={4} input={-5.28951,} output={-0.125,} objective={-0.125}
rank={0,1,} iter={5} input={-43.1005,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={6} input={-46.0517,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={7} input={-44.5683,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={8} input={-45.308,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={9} input={-43.8341,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={10} input={-45.6812,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={11} input={-44.9298,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={12} input={-44.2051,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={13} input={-43.4633,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={14} input={-45.115,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={15} input={-45.4943,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={16} input={-43.6472,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={17} input={-44.0203,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={18} input={-45.8655,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={19} input={-44.3882,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={20} input={-44.7497,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={21} input={-43.2885,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={22} input={-45.2118,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={23} input={-45.5877,} output={3.76471,} objective={3.76471}
rank={0,1,} iter={24} input={-43.1952,} output={3.76471,} objective={3.76471}
final_iter={25} inputs={-43.1005,} output={3.76471,}
rank={0,1,} iter={0} input={-24.1771,} output={-0.0694444,} objective={-0.0694444}
rank={0,1,} iter={1} input={-35.8462,} output={-0.236111,} objective={-0.236111}
rank={0,1,} iter={2} input={-12.6864,} output={-0.236111,} objective={-0.236111}
rank={0,1,} iter={3} input={-36.7947,} output={-0.180556,} objective={-0.180556}
rank={0,1,} iter={4} input={-5.28951,} output={-0.236111,} objective={-0.236111}
rank={0,1,} iter={5} input={-43.1005,} output={1.02674,} objective={1.02674}
rank={0,1,} iter={6} input={-46.0517,} output={0.998267,} objective={0.998267}
rank={0,1,} iter={7} input={-44.4883,} output={1.0123,} objective={1.0123}
rank={0,1,} iter={8} input={-39.6547,} output={-0.236111,} objective={-0.236111}
rank={0,1,} iter={9} input={-30.0126,} output={-0.180556,} objective={-0.180556}
rank={0,1,} iter={10} input={-41.3776,} output={-0.236111,} objective={-0.236111}
rank={0,1,} iter={11} input={-18.4353,} output={-0.236111,} objective={-0.236111}
rank={0,1,} iter={12} input={-43.7677,} output={1.0123,} objective={1.0123}
rank={0,1,} iter={13} input={-8.9921,} output={-0.236111,} objective={-0.236111}
rank={0,1,} iter={14} input={-42.6698,} output={1.04159,} objective={1.04159}
rank={0,1,} iter={15} input={-2.32243,} output={-0.180556,} objective={-0.180556}
rank={0,1,} iter={16} input={-41.8084,} output={-0.166667,} objective={-0.166667}
rank={0,1,} iter={17} input={-27.092,} output={-0.208333,} objective={-0.208333}
rank={0,1,} iter={18} input={-42.8673,} output={1.04159,} objective={1.04159}
rank={0,1,} iter={19} input={-32.9292,} output={-0.236111,} objective={-0.236111}
rank={0,1,} iter={20} input={-15.5603,} output={-0.236111,} objective={-0.236111}
rank={0,1,} iter={21} input={-21.3068,} output={-0.222222,} objective={-0.222222}
rank={0,1,} iter={22} input={-7.14094,} output={-0.222222,} objective={-0.222222}
rank={0,1,} iter={23} input={-10.8504,} output={-0.0694444,} objective={-0.0694444}
rank={0,1,} iter={24} input={-45.2612,} output={0.998267,} objective={0.998267}
final_iter={25} inputs={-42.6698,} output={1.04159,}
rank={0,1,} iter={0} input={-8.05905,} output={4.608,} objective={4.608}
rank={0,1,} iter={1} input={-11.1299,} output={4.608,} objective={4.608}
rank={0,1,} iter={2} input={-5.03517,} output={4.608,} objective={4.608}
rank={0,1,} iter={3} input={-11.3795,} output={4.608,} objective={4.608}
rank={0,1,} iter={4} input={-3.08862,} output={4.608,} objective={4.608}
rank={0,1,} iter={5} input={-13.0389,} output={4.608,} objective={4.608}
rank={0,1,} iter={6} input={-3.44921,} output={4.608,} objective={4.608}
rank={0,1,} iter={7} input={-3.22263,} output={4.608,} objective={4.608}
rank={0,1,} iter={8} input={-8.16974,} output={4.608,} objective={4.608}
rank={0,1,} iter={9} input={-11.9865,} output={4.608,} objective={4.608}
rank={0,1,} iter={10} input={-4.47192,} output={4.608,} objective={4.608}
rank={0,1,} iter={11} input={-8.09914,} output={4.608,} objective={4.608}
rank={0,1,} iter={12} input={-12.1194,} output={4.608,} objective={4.608}
rank={0,1,} iter={13} input={-4.14904,} output={4.608,} objective={4.608}
rank={0,1,} iter={14} input={-12.852,} output={4.608,} objective={4.608}
rank={0,1,} iter={15} input={-11.2082,} output={4.608,} objective={4.608}
rank={0,1,} iter={16} input={-9.64318,} output={4.608,} objective={4.608}
rank={0,1,} iter={17} input={-2.45887,} output={4.608,} objective={4.608}
rank={0,1,} iter={18} input={-7.86641,} output={4.608,} objective={4.608}
rank={0,1,} iter={19} input={-3.29203,} output={4.608,} objective={4.608}
rank={0,1,} iter={20} input={-13.7224,} output={-0.0277778,} objective={-0.0277778}
rank={0,1,} iter={21} input={-6.45527,} output={4.608,} objective={4.608}
rank={0,1,} iter={22} input={-10.3863,} output={4.608,} objective={4.608}
rank={0,1,} iter={23} input={-8.90622,} output={4.608,} objective={4.608}
rank={0,1,} iter={24} input={-5.74466,} output={4.608,} objective={4.608}
final_iter={25} inputs={-8.05905,} output={4.608,}
rank={0,1,} iter={0} input={-8.05905,} output={2,} objective={2}
rank={0,1,} iter={1} input={-11.1299,} output={2,} objective={2}
rank={0,1,} iter={2} input={-5.03517,} output={2,} objective={2}
rank={0,1,} iter={3} input={-11.3795,} output={2,} objective={2}
rank={0,1,} iter={4} input={-3.08862,} output={2,} objective={2}
rank={0,1,} iter={5} input={-13.0389,} output={2,} objective={2}
rank={0,1,} iter={6} input={-3.44921,} output={2,} objective={2}
rank={0,1,} iter={7} input={-3.22263,} output={2,} objective={2}
rank={0,1,} iter={8} input={-8.16974,} output={2,} objective={2}
rank={0,1,} iter={9} input={-11.9865,} output={2,} objective={2}
rank={0,1,} iter={10} input={-4.47192,} output={2,} objective={2}
rank={0,1,} iter={11} input={-8.09914,} output={2,} objective={2}
rank={0,1,} iter={12} input={-12.1194,} output={2,} objective={2}
rank={0,1,} iter={13} input={-4.14904,} output={2,} objective={2}
rank={0,1,} iter={14} input={-12.852,} output={2,} objective={2}
rank={0,1,} iter={15} input={-11.2082,} output={2,} objective={2}
rank={0,1,} iter={16} input={-9.64318,} output={2,} objective={2}
rank={0,1,} iter={17} input={-2.45887,} output={-0.361111,} objective={-0.361111}
rank={0,1,} iter={18} input={-6.54644,} output={2,} objective={2}
rank={0,1,} iter={19} input={-13.8132,} output={2,} objective={2}
rank={0,1,} iter={20} input={-7.30218,} output={2,} objective={2}
rank={0,1,} iter={21} input={-5.7914,} output={2,} objective={2}
rank={0,1,} iter={22} input={-10.3863,} output={2,} objective={2}
rank={0,1,} iter={23} input={-8.90622,} output={2,} objective={2}
rank={0,1,} iter={24} input={-13.4257,} output={2,} objective={2}
final_iter={25} inputs={-8.05905,} output={2,}
rank={0,1,} iter={0} input={-19.572,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={1} input={-28.7844,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={2} input={-10.5004,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={3} input={-29.5332,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={4} input={-4.66068,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={5} input={-34.5115,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={6} input={-5.74245,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={7} input={-5.06272,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={8} input={-19.904,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={9} input={-31.3544,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={10} input={-8.81058,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={11} input={-19.6922,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={12} input={-31.7529,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={13} input={-7.84195,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={14} input={-33.9509,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={15} input={-29.0194,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={16} input={-24.3244,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={17} input={-2.77145,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={18} input={-18.9941,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={19} input={-5.27091,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={20} input={-36.5621,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={21} input={-31.2957,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={22} input={-8.70096,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={23} input={-34.0509,} output={-0.902778,} objective={-0.902778}
rank={0,1,} iter={24} input={-3.32995,} output={-0.902778,} objective={-0.902778}
final_iter={25} inputs={-19.572,} output={-0.902778,}
rank={0,1,} iter={0} input={-19.572,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={1} input={-28.7844,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={2} input={-10.5004,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={3} input={-29.5332,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={4} input={-4.66068,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={5} input={-34.5115,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={6} input={-5.74245,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={7} input={-5.06272,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={8} input={-19.904,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={9} input={-31.3544,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={10} input={-8.81058,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={11} input={-19.6922,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={12} input={-31.7529,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={13} input={-7.84195,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={14} input={-33.9509,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={15} input={-29.0194,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={16} input={-24.3244,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={17} input={-2.77145,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={18} input={-18.9941,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={19} input={-5.27091,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={20} input={-36.5621,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={21} input={-31.2957,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={22} input={-8.70096,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={23} input={-34.0509,} output={-0.291667,} objective={-0.291667}
rank={0,1,} iter={24} input={-3.32995,} output={-0.416667,} objective={-0.416667}
final_iter={25} inputs={-19.572,} output={-0.291667,}
Visual comparison of the precipitation time series¶
In the following analyses, we plot the time series (error) of the observed and modelled precipitation in Helsinki, Finland and Belém, Brazil. On each time series, we show the time and magnitude of the maximum as a dot, the filled integral under the curve. Below the time series, we also highlight times at which precipitation occurred, i.e. was positive. Since negative precipitation does not occur but lossy compression can interpolate negative precipitation values, we also highlight time steps with these errors in red.
It is worth noting that the ERA5 time series is more smoothed out than the original observations and includes much longer precipitation events with lower extrema.
All three compressors preserve the global maximum well (only SZ3 changes the time of maxima by one hour for the Helsinki observations) and produce intervals without systematic bias. However, all compressors struggle with preserving zero and positive values as such and with not producing negative values. Since ZFP is a blockwise compressor, it has the least artefacts. SZ3 and SPERR, which are predictive and wavelet-transform compressors, respectively, produce significant negative precipitation values, with SPERR performing worse since its transform is global.
The safeguards correctly preserve all properties of interest. However, the total precipitation integral has a consistent negative bias when safeguarding a constant-zero approximation.
Copied!
fig1, axs1 = plt.subplots(6, 2, figsize=(12, 23))
fig2, axs2 = plt.subplots(6, 2, figsize=(12, 23))
plot_precipitation(
axs1[0, 0],
axs2[0, 0],
ERA5_PR,
1.0,
Belem_PR,
1.0,
Helsinki_PR,
1.0,
eb_abs,
"Original",
reference=True,
)
plot_precipitation(
axs1[1, 0],
axs2[1, 0],
ERA5_PR_zfp,
ERA5_PR_zfp_cr,
Belem_PR_zfp,
Belem_PR_zfp_cr,
Helsinki_PR_zfp,
Helsinki_PR_zfp_cr,
eb_abs,
r"ZFP($\epsilon_{{abs}}$)",
)
plot_precipitation(
axs1[2, 0],
axs2[2, 0],
ERA5_PR_sz3,
ERA5_PR_sz3_cr,
Belem_PR_sz3,
Belem_PR_sz3_cr,
Helsinki_PR_sz3,
Helsinki_PR_sz3_cr,
eb_abs,
r"SZ3($\epsilon_{{abs}}$)",
)
plot_precipitation(
axs1[3, 0],
axs2[3, 0],
ERA5_PR_sperr,
ERA5_PR_sperr_cr,
Belem_PR_sperr,
Belem_PR_sperr_cr,
Helsinki_PR_sperr,
Helsinki_PR_sperr_cr,
eb_abs,
r"SPERR($\epsilon_{{abs}}$)",
)
plot_precipitation(
axs1[0, 1],
axs2[0, 1],
ERA5_PR_sg["zero"],
ERA5_PR_sg_cr["zero"],
Belem_PR_sg["zero"],
Belem_PR_sg_cr["zero"],
Helsinki_PR_sg["zero"],
Helsinki_PR_sg_cr["zero"],
eb_abs,
r"Safeguarded(0, $\epsilon_{{abs}} \cup \text{sgn}(PR) \cup \max(PR)$)",
corr=(ERA5_PR_zero, Belem_PR_zero, Helsinki_PR_zero),
)
plot_precipitation(
axs1[1, 1],
axs2[1, 1],
ERA5_PR_sg["zfp.rs"],
ERA5_PR_sg_cr["zfp.rs"],
Belem_PR_sg["zfp.rs"],
Belem_PR_sg_cr["zfp.rs"],
Helsinki_PR_sg["zfp.rs"],
Helsinki_PR_sg_cr["zfp.rs"],
eb_abs,
r"Safeguarded(ZFP, $\epsilon_{{abs}} \cup \text{sgn}(PR) \cup \max(PR)$)",
corr=(ERA5_PR_zfp, Belem_PR_zfp, Helsinki_PR_zfp),
)
plot_precipitation(
axs1[2, 1],
axs2[2, 1],
ERA5_PR_sg["sz3.rs"],
ERA5_PR_sg_cr["sz3.rs"],
Belem_PR_sg["sz3.rs"],
Belem_PR_sg_cr["sz3.rs"],
Helsinki_PR_sg["sz3.rs"],
Helsinki_PR_sg_cr["sz3.rs"],
eb_abs,
r"Safeguarded(SZ3, $\epsilon_{{abs}} \cup \text{sgn}(PR) \cup \max(PR)$)",
corr=(ERA5_PR_sz3, Belem_PR_sz3, Helsinki_PR_sz3),
)
plot_precipitation(
axs1[3, 1],
axs2[3, 1],
ERA5_PR_sg["sperr.rs"],
ERA5_PR_sg_cr["sperr.rs"],
Belem_PR_sg["sperr.rs"],
Belem_PR_sg_cr["sperr.rs"],
Helsinki_PR_sg["sperr.rs"],
Helsinki_PR_sg_cr["sperr.rs"],
eb_abs,
r"Safeguarded(SPERR, $\epsilon_{{abs}} \cup \text{sgn}(PR) \cup \max(PR)$)",
corr=(ERA5_PR_sperr, Belem_PR_sperr, Helsinki_PR_sperr),
)
plot_precipitation(
axs1[4, 0],
axs2[4, 0],
ERA5_PR_optzconfig["zfp.rs"],
ERA5_PR_optzconfig_cr["zfp.rs"],
Belem_PR_optzconfig["zfp.rs"],
Belem_PR_optzconfig_cr["zfp.rs"],
Helsinki_PR_optzconfig["zfp.rs"],
Helsinki_PR_optzconfig_cr["zfp.rs"],
eb_abs,
r"OptZConfig(ZFP, $\epsilon_{{abs}} \cup \text{sgn}(PR)$)",
)
plot_precipitation(
axs1[4, 1],
axs2[4, 1],
ERA5_PR_optzconfig["sz3.rs"],
ERA5_PR_optzconfig_cr["sz3.rs"],
Belem_PR_optzconfig["sz3.rs"],
Belem_PR_optzconfig_cr["sz3.rs"],
Helsinki_PR_optzconfig["sz3.rs"],
Helsinki_PR_optzconfig_cr["sz3.rs"],
eb_abs,
r"OptZConfig(SZ3, $\epsilon_{{abs}} \cup \text{sgn}(PR)$)",
)
plot_precipitation(
axs1[5, 0],
axs2[5, 0],
ERA5_PR_optzconfig["sperr.rs"],
ERA5_PR_optzconfig_cr["sperr.rs"],
Belem_PR_optzconfig["sperr.rs"],
Belem_PR_optzconfig_cr["sperr.rs"],
Helsinki_PR_optzconfig["sperr.rs"],
Helsinki_PR_optzconfig_cr["sperr.rs"],
eb_abs,
r"OptZConfig(SPERR, $\epsilon_{{abs}} \cup \text{sgn}(PR)$)",
)
axs1[5, 1].set_axis_off()
axs2[5, 1].set_axis_off()
fig1.tight_layout()
fig2.tight_layout()
fig1.savefig(Path("plots") / "precipitation-obs.pdf")
fig2.savefig(Path("plots") / "precipitation-era5.pdf")
plt.show()
fig1, axs1 = plt.subplots(6, 2, figsize=(12, 23))
fig2, axs2 = plt.subplots(6, 2, figsize=(12, 23))
plot_precipitation(
axs1[0, 0],
axs2[0, 0],
ERA5_PR,
1.0,
Belem_PR,
1.0,
Helsinki_PR,
1.0,
eb_abs,
"Original",
reference=True,
)
plot_precipitation(
axs1[1, 0],
axs2[1, 0],
ERA5_PR_zfp,
ERA5_PR_zfp_cr,
Belem_PR_zfp,
Belem_PR_zfp_cr,
Helsinki_PR_zfp,
Helsinki_PR_zfp_cr,
eb_abs,
r"ZFP($\epsilon_{{abs}}$)",
)
plot_precipitation(
axs1[2, 0],
axs2[2, 0],
ERA5_PR_sz3,
ERA5_PR_sz3_cr,
Belem_PR_sz3,
Belem_PR_sz3_cr,
Helsinki_PR_sz3,
Helsinki_PR_sz3_cr,
eb_abs,
r"SZ3($\epsilon_{{abs}}$)",
)
plot_precipitation(
axs1[3, 0],
axs2[3, 0],
ERA5_PR_sperr,
ERA5_PR_sperr_cr,
Belem_PR_sperr,
Belem_PR_sperr_cr,
Helsinki_PR_sperr,
Helsinki_PR_sperr_cr,
eb_abs,
r"SPERR($\epsilon_{{abs}}$)",
)
plot_precipitation(
axs1[0, 1],
axs2[0, 1],
ERA5_PR_sg["zero"],
ERA5_PR_sg_cr["zero"],
Belem_PR_sg["zero"],
Belem_PR_sg_cr["zero"],
Helsinki_PR_sg["zero"],
Helsinki_PR_sg_cr["zero"],
eb_abs,
r"Safeguarded(0, $\epsilon_{{abs}} \cup \text{sgn}(PR) \cup \max(PR)$)",
corr=(ERA5_PR_zero, Belem_PR_zero, Helsinki_PR_zero),
)
plot_precipitation(
axs1[1, 1],
axs2[1, 1],
ERA5_PR_sg["zfp.rs"],
ERA5_PR_sg_cr["zfp.rs"],
Belem_PR_sg["zfp.rs"],
Belem_PR_sg_cr["zfp.rs"],
Helsinki_PR_sg["zfp.rs"],
Helsinki_PR_sg_cr["zfp.rs"],
eb_abs,
r"Safeguarded(ZFP, $\epsilon_{{abs}} \cup \text{sgn}(PR) \cup \max(PR)$)",
corr=(ERA5_PR_zfp, Belem_PR_zfp, Helsinki_PR_zfp),
)
plot_precipitation(
axs1[2, 1],
axs2[2, 1],
ERA5_PR_sg["sz3.rs"],
ERA5_PR_sg_cr["sz3.rs"],
Belem_PR_sg["sz3.rs"],
Belem_PR_sg_cr["sz3.rs"],
Helsinki_PR_sg["sz3.rs"],
Helsinki_PR_sg_cr["sz3.rs"],
eb_abs,
r"Safeguarded(SZ3, $\epsilon_{{abs}} \cup \text{sgn}(PR) \cup \max(PR)$)",
corr=(ERA5_PR_sz3, Belem_PR_sz3, Helsinki_PR_sz3),
)
plot_precipitation(
axs1[3, 1],
axs2[3, 1],
ERA5_PR_sg["sperr.rs"],
ERA5_PR_sg_cr["sperr.rs"],
Belem_PR_sg["sperr.rs"],
Belem_PR_sg_cr["sperr.rs"],
Helsinki_PR_sg["sperr.rs"],
Helsinki_PR_sg_cr["sperr.rs"],
eb_abs,
r"Safeguarded(SPERR, $\epsilon_{{abs}} \cup \text{sgn}(PR) \cup \max(PR)$)",
corr=(ERA5_PR_sperr, Belem_PR_sperr, Helsinki_PR_sperr),
)
plot_precipitation(
axs1[4, 0],
axs2[4, 0],
ERA5_PR_optzconfig["zfp.rs"],
ERA5_PR_optzconfig_cr["zfp.rs"],
Belem_PR_optzconfig["zfp.rs"],
Belem_PR_optzconfig_cr["zfp.rs"],
Helsinki_PR_optzconfig["zfp.rs"],
Helsinki_PR_optzconfig_cr["zfp.rs"],
eb_abs,
r"OptZConfig(ZFP, $\epsilon_{{abs}} \cup \text{sgn}(PR)$)",
)
plot_precipitation(
axs1[4, 1],
axs2[4, 1],
ERA5_PR_optzconfig["sz3.rs"],
ERA5_PR_optzconfig_cr["sz3.rs"],
Belem_PR_optzconfig["sz3.rs"],
Belem_PR_optzconfig_cr["sz3.rs"],
Helsinki_PR_optzconfig["sz3.rs"],
Helsinki_PR_optzconfig_cr["sz3.rs"],
eb_abs,
r"OptZConfig(SZ3, $\epsilon_{{abs}} \cup \text{sgn}(PR)$)",
)
plot_precipitation(
axs1[5, 0],
axs2[5, 0],
ERA5_PR_optzconfig["sperr.rs"],
ERA5_PR_optzconfig_cr["sperr.rs"],
Belem_PR_optzconfig["sperr.rs"],
Belem_PR_optzconfig_cr["sperr.rs"],
Helsinki_PR_optzconfig["sperr.rs"],
Helsinki_PR_optzconfig_cr["sperr.rs"],
eb_abs,
r"OptZConfig(SPERR, $\epsilon_{{abs}} \cup \text{sgn}(PR)$)",
)
axs1[5, 1].set_axis_off()
axs2[5, 1].set_axis_off()
fig1.tight_layout()
fig2.tight_layout()
fig1.savefig(Path("plots") / "precipitation-obs.pdf")
fig2.savefig(Path("plots") / "precipitation-era5.pdf")
plt.show()
Copied!
pr_table = pd.concat(
[
table_precipitation(
ERA5_PR,
ERA5_PR.nbytes,
Belem_PR,
Belem_PR.nbytes,
Helsinki_PR,
Helsinki_PR.nbytes,
["Original", "-", ""],
eb_abs,
None,
reference=True,
),
table_precipitation(
ERA5_PR_sg_lossless["zero"],
ERA5_PR_sg_lossless_cr["zero"],
Belem_PR_sg_lossless["zero"],
Belem_PR_sg_lossless_cr["zero"],
Helsinki_PR_sg_lossless["zero"],
Helsinki_PR_sg_lossless_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$", "lossless"],
eb_abs,
(ERA5_PR_zero, Belem_PR_zero, Helsinki_PR_zero),
),
table_precipitation(
ERA5_PR_sg["zero"],
ERA5_PR_sg_cr["zero"],
Belem_PR_sg["zero"],
Belem_PR_sg_cr["zero"],
Helsinki_PR_sg["zero"],
Helsinki_PR_sg_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$", "one-shot"],
eb_abs,
(ERA5_PR_zero, Belem_PR_zero, Helsinki_PR_zero),
),
table_precipitation(
ERA5_PR_zfp,
ERA5_PR_zfp_cr,
Belem_PR_zfp,
Belem_PR_zfp_cr,
Helsinki_PR_zfp,
Helsinki_PR_zfp_cr,
[r"ZFP($\epsilon_{abs}$)", "-", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sg_lossless["zfp.rs"],
ERA5_PR_sg_lossless_cr["zfp.rs"],
Belem_PR_sg_lossless["zfp.rs"],
Belem_PR_sg_lossless_cr["zfp.rs"],
Helsinki_PR_sg_lossless["zfp.rs"],
Helsinki_PR_sg_lossless_cr["zfp.rs"],
[
r"ZFP($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"lossless",
],
eb_abs,
(ERA5_PR_zfp, Belem_PR_zfp, Helsinki_PR_zfp),
),
table_precipitation(
ERA5_PR_sg["zfp.rs"],
ERA5_PR_sg_cr["zfp.rs"],
Belem_PR_sg["zfp.rs"],
Belem_PR_sg_cr["zfp.rs"],
Helsinki_PR_sg["zfp.rs"],
Helsinki_PR_sg_cr["zfp.rs"],
[
r"ZFP($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"one-shot",
],
eb_abs,
(ERA5_PR_zfp, Belem_PR_zfp, Helsinki_PR_zfp),
),
table_precipitation(
ERA5_PR_optzconfig["zfp.rs"],
ERA5_PR_optzconfig_cr["zfp.rs"],
Belem_PR_optzconfig["zfp.rs"],
Belem_PR_optzconfig_cr["zfp.rs"],
Helsinki_PR_optzconfig["zfp.rs"],
Helsinki_PR_optzconfig_cr["zfp.rs"],
["OptZConfig(ZFP)", r"$\epsilon_{abs} \cup \text{sgn}(PR)$", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sz3,
ERA5_PR_sz3_cr,
Belem_PR_sz3,
Belem_PR_sz3_cr,
Helsinki_PR_sz3,
Helsinki_PR_sz3_cr,
[r"SZ3($\epsilon_{abs}$)", "-", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sg_lossless["sz3.rs"],
ERA5_PR_sg_lossless_cr["sz3.rs"],
Belem_PR_sg_lossless["sz3.rs"],
Belem_PR_sg_lossless_cr["sz3.rs"],
Helsinki_PR_sg_lossless["sz3.rs"],
Helsinki_PR_sg_lossless_cr["sz3.rs"],
[
r"SZ3($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"lossless",
],
eb_abs,
(ERA5_PR_sz3, Belem_PR_sz3, Helsinki_PR_sz3),
),
table_precipitation(
ERA5_PR_sg["sz3.rs"],
ERA5_PR_sg_cr["sz3.rs"],
Belem_PR_sg["sz3.rs"],
Belem_PR_sg_cr["sz3.rs"],
Helsinki_PR_sg["sz3.rs"],
Helsinki_PR_sg_cr["sz3.rs"],
[
r"SZ3($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"one-shot",
],
eb_abs,
(ERA5_PR_sz3, Belem_PR_sz3, Helsinki_PR_sz3),
),
table_precipitation(
ERA5_PR_optzconfig["sz3.rs"],
ERA5_PR_optzconfig_cr["sz3.rs"],
Belem_PR_optzconfig["sz3.rs"],
Belem_PR_optzconfig_cr["sz3.rs"],
Helsinki_PR_optzconfig["sz3.rs"],
Helsinki_PR_optzconfig_cr["sz3.rs"],
["OptZConfig(SZ3)", r"$\epsilon_{abs} \cup \text{sgn}(PR)$", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sperr,
ERA5_PR_sperr_cr,
Belem_PR_sperr,
Belem_PR_sperr_cr,
Helsinki_PR_sperr,
Helsinki_PR_sperr_cr,
[r"SPERR($\epsilon_{abs}$)", "-", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sg_lossless["sperr.rs"],
ERA5_PR_sg_lossless_cr["sperr.rs"],
Belem_PR_sg_lossless["sperr.rs"],
Belem_PR_sg_lossless_cr["sperr.rs"],
Helsinki_PR_sg_lossless["sperr.rs"],
Helsinki_PR_sg_lossless_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"lossless",
],
eb_abs,
(ERA5_PR_sperr, Belem_PR_sperr, Helsinki_PR_sperr),
),
table_precipitation(
ERA5_PR_sg["sperr.rs"],
ERA5_PR_sg_cr["sperr.rs"],
Belem_PR_sg["sperr.rs"],
Belem_PR_sg_cr["sperr.rs"],
Helsinki_PR_sg["sperr.rs"],
Helsinki_PR_sg_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"one-shot",
],
eb_abs,
(ERA5_PR_sperr, Belem_PR_sperr, Helsinki_PR_sperr),
),
table_precipitation(
ERA5_PR_optzconfig["sperr.rs"],
ERA5_PR_optzconfig_cr["sperr.rs"],
Belem_PR_optzconfig["sperr.rs"],
Belem_PR_optzconfig_cr["sperr.rs"],
Helsinki_PR_optzconfig["sperr.rs"],
Helsinki_PR_optzconfig_cr["sperr.rs"],
["OptZConfig(SPERR)", r"$\epsilon_{abs} \cup \text{sgn}(PR)$", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_zstd,
ERA5_PR_zstd_cr,
Belem_PR_zstd,
Belem_PR_zstd_cr,
Helsinki_PR_zstd,
Helsinki_PR_zstd_cr,
["ZSTD(22)", "-", ""],
eb_abs,
None,
),
]
).set_index(["Compressor", "Safeguarded", "Corrections", "Station", "Source"])
Path("tables").joinpath("precipitation.tex").write_text(
pr_table.to_latex(escape=False)
.replace("%", r"\%")
.replace(
"\\cline{1-16} \\cline{2-16} \\cline{3-16} \\cline{4-16}\n\\bottomrule",
"\\bottomrule",
)
)
pr_table
pr_table = pd.concat(
[
table_precipitation(
ERA5_PR,
ERA5_PR.nbytes,
Belem_PR,
Belem_PR.nbytes,
Helsinki_PR,
Helsinki_PR.nbytes,
["Original", "-", ""],
eb_abs,
None,
reference=True,
),
table_precipitation(
ERA5_PR_sg_lossless["zero"],
ERA5_PR_sg_lossless_cr["zero"],
Belem_PR_sg_lossless["zero"],
Belem_PR_sg_lossless_cr["zero"],
Helsinki_PR_sg_lossless["zero"],
Helsinki_PR_sg_lossless_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$", "lossless"],
eb_abs,
(ERA5_PR_zero, Belem_PR_zero, Helsinki_PR_zero),
),
table_precipitation(
ERA5_PR_sg["zero"],
ERA5_PR_sg_cr["zero"],
Belem_PR_sg["zero"],
Belem_PR_sg_cr["zero"],
Helsinki_PR_sg["zero"],
Helsinki_PR_sg_cr["zero"],
["0", r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$", "one-shot"],
eb_abs,
(ERA5_PR_zero, Belem_PR_zero, Helsinki_PR_zero),
),
table_precipitation(
ERA5_PR_zfp,
ERA5_PR_zfp_cr,
Belem_PR_zfp,
Belem_PR_zfp_cr,
Helsinki_PR_zfp,
Helsinki_PR_zfp_cr,
[r"ZFP($\epsilon_{abs}$)", "-", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sg_lossless["zfp.rs"],
ERA5_PR_sg_lossless_cr["zfp.rs"],
Belem_PR_sg_lossless["zfp.rs"],
Belem_PR_sg_lossless_cr["zfp.rs"],
Helsinki_PR_sg_lossless["zfp.rs"],
Helsinki_PR_sg_lossless_cr["zfp.rs"],
[
r"ZFP($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"lossless",
],
eb_abs,
(ERA5_PR_zfp, Belem_PR_zfp, Helsinki_PR_zfp),
),
table_precipitation(
ERA5_PR_sg["zfp.rs"],
ERA5_PR_sg_cr["zfp.rs"],
Belem_PR_sg["zfp.rs"],
Belem_PR_sg_cr["zfp.rs"],
Helsinki_PR_sg["zfp.rs"],
Helsinki_PR_sg_cr["zfp.rs"],
[
r"ZFP($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"one-shot",
],
eb_abs,
(ERA5_PR_zfp, Belem_PR_zfp, Helsinki_PR_zfp),
),
table_precipitation(
ERA5_PR_optzconfig["zfp.rs"],
ERA5_PR_optzconfig_cr["zfp.rs"],
Belem_PR_optzconfig["zfp.rs"],
Belem_PR_optzconfig_cr["zfp.rs"],
Helsinki_PR_optzconfig["zfp.rs"],
Helsinki_PR_optzconfig_cr["zfp.rs"],
["OptZConfig(ZFP)", r"$\epsilon_{abs} \cup \text{sgn}(PR)$", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sz3,
ERA5_PR_sz3_cr,
Belem_PR_sz3,
Belem_PR_sz3_cr,
Helsinki_PR_sz3,
Helsinki_PR_sz3_cr,
[r"SZ3($\epsilon_{abs}$)", "-", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sg_lossless["sz3.rs"],
ERA5_PR_sg_lossless_cr["sz3.rs"],
Belem_PR_sg_lossless["sz3.rs"],
Belem_PR_sg_lossless_cr["sz3.rs"],
Helsinki_PR_sg_lossless["sz3.rs"],
Helsinki_PR_sg_lossless_cr["sz3.rs"],
[
r"SZ3($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"lossless",
],
eb_abs,
(ERA5_PR_sz3, Belem_PR_sz3, Helsinki_PR_sz3),
),
table_precipitation(
ERA5_PR_sg["sz3.rs"],
ERA5_PR_sg_cr["sz3.rs"],
Belem_PR_sg["sz3.rs"],
Belem_PR_sg_cr["sz3.rs"],
Helsinki_PR_sg["sz3.rs"],
Helsinki_PR_sg_cr["sz3.rs"],
[
r"SZ3($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"one-shot",
],
eb_abs,
(ERA5_PR_sz3, Belem_PR_sz3, Helsinki_PR_sz3),
),
table_precipitation(
ERA5_PR_optzconfig["sz3.rs"],
ERA5_PR_optzconfig_cr["sz3.rs"],
Belem_PR_optzconfig["sz3.rs"],
Belem_PR_optzconfig_cr["sz3.rs"],
Helsinki_PR_optzconfig["sz3.rs"],
Helsinki_PR_optzconfig_cr["sz3.rs"],
["OptZConfig(SZ3)", r"$\epsilon_{abs} \cup \text{sgn}(PR)$", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sperr,
ERA5_PR_sperr_cr,
Belem_PR_sperr,
Belem_PR_sperr_cr,
Helsinki_PR_sperr,
Helsinki_PR_sperr_cr,
[r"SPERR($\epsilon_{abs}$)", "-", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_sg_lossless["sperr.rs"],
ERA5_PR_sg_lossless_cr["sperr.rs"],
Belem_PR_sg_lossless["sperr.rs"],
Belem_PR_sg_lossless_cr["sperr.rs"],
Helsinki_PR_sg_lossless["sperr.rs"],
Helsinki_PR_sg_lossless_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"lossless",
],
eb_abs,
(ERA5_PR_sperr, Belem_PR_sperr, Helsinki_PR_sperr),
),
table_precipitation(
ERA5_PR_sg["sperr.rs"],
ERA5_PR_sg_cr["sperr.rs"],
Belem_PR_sg["sperr.rs"],
Belem_PR_sg_cr["sperr.rs"],
Helsinki_PR_sg["sperr.rs"],
Helsinki_PR_sg_cr["sperr.rs"],
[
r"SPERR($\epsilon_{abs}$)",
r"$\epsilon_{abs} \cup \text{sgn}(PR) \cup \max(PR)$",
"one-shot",
],
eb_abs,
(ERA5_PR_sperr, Belem_PR_sperr, Helsinki_PR_sperr),
),
table_precipitation(
ERA5_PR_optzconfig["sperr.rs"],
ERA5_PR_optzconfig_cr["sperr.rs"],
Belem_PR_optzconfig["sperr.rs"],
Belem_PR_optzconfig_cr["sperr.rs"],
Helsinki_PR_optzconfig["sperr.rs"],
Helsinki_PR_optzconfig_cr["sperr.rs"],
["OptZConfig(SPERR)", r"$\epsilon_{abs} \cup \text{sgn}(PR)$", ""],
eb_abs,
None,
),
table_precipitation(
ERA5_PR_zstd,
ERA5_PR_zstd_cr,
Belem_PR_zstd,
Belem_PR_zstd_cr,
Helsinki_PR_zstd,
Helsinki_PR_zstd_cr,
["ZSTD(22)", "-", ""],
eb_abs,
None,
),
]
).set_index(["Compressor", "Safeguarded", "Corrections", "Station", "Source"])
Path("tables").joinpath("precipitation.tex").write_text(
pr_table.to_latex(escape=False)
.replace("%", r"\%")
.replace(
"\\cline{1-16} \\cline{2-16} \\cline{3-16} \\cline{4-16}\n\\bottomrule",
"\\bottomrule",
)
)
pr_table
| $L_{\infty}(\hat{PR})$ | $L_{2}(\hat{PR})$ | Integral | max(PR) | argmax(PR) | PR > 0 | FP + FN | PR < 0 | V | C | CR | |||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Compressor | Safeguarded | Corrections | Station | Source | |||||||||||
| Original | - | Helsinki | obs | 37.9mm | 4.44mm | 02.04 20h | 51h | 0h | 0h | $\times$ 1 | |||||
| ERA5 | 24.8mm | 3.00mm | 02.04 20h | 39h | 0h | 0h | |||||||||
| Belém | obs | 31.8mm | 22.20mm | 03.04 17h | 7h | 0h | 0h | $\times$ 1 | |||||||
| ERA5 | 26.0mm | 2.71mm | 03.04 18h | 71h | 0h | 0h | |||||||||
| ERA5 | avg | 7.1mm | 0.87mm | 0.65% | 0h | 0h | $\times$ 1 | ||||||||
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| ZSTD(22) | - | Helsinki | obs | 0.0 | 0.0 | +0mm | +0mm | +0h | +0h | 0h | 0h | 0 | $\times$ 2.74 | ||
| ERA5 | 0.0 | 0.0 | +0mm | +0mm | +0h | +0h | 0h | 0h | 0 | ||||||
| Belém | obs | 0.0 | 0.0 | +0mm | +0mm | +0h | +0h | 0h | 0h | 0 | $\times$ 8.11 | ||||
| ERA5 | 0.0 | 0.0 | +0mm | +0mm | +0h | +0h | 0h | 0h | 0 | ||||||
| ERA5 | avg | 0.0 | 0.0 | +0mm | +0mm | +0h | 0.65% | 0h | 0h | 0 | $\times$ 5.42 |
80 rows × 11 columns
Copied!
import json
with Path("observations").joinpath("precipitation.json").open("w") as f:
json.dump(observations, f)
import json
with Path("observations").joinpath("precipitation.json").open("w") as f:
json.dump(observations, f)
Copied!