[Python] Scaling a GEV distribution for precipitation data by specifying a covariate-dependent...

Background: I am conducting a Climate Attribution study for which am applying the methods proposed in the protocol published by Philip et al. (2020). The data I am working with is daily precipitation data from a weather station over a roughly 100-year time period. The special focus lies on an heavy precipitation event that occured in 1998. The goal is to find out whether the likelihood of occurency and the magnitude of this event changed over the years. "The result of the fit to a statistical model provides an estimate of the return period of the event of interest and whether there is a trend outside the range of deviations expected by natural variability." Block maxima (per year) are computed from the precipitation data. This extreme value dataset is subsequently fitted to a Generalized extreme value (GEV) distribution and scaled it with an exponential function of the smoothed global mean temperature (GMST). The whole distribution is thus scaled with the location and scale parameters expressed in the following functions: µ = µ0 exp(αT'/µ0) and σ = σ0 exp(αT'/µ0) "The covariate-dependent function can be inverted and the distribution evaluated for a given year, e.g. a year in the past (with T' = T0) or the current year (T' = T1). This gives the probabilities for an event at least as extreme as the observed one in these 2 years, p0 and p1, or expressed as return periods τ0 = 1/p0 and τ1 = 1/p1."

Problem: My problem lies in the application of the exponential scaling functions. I have fitted the stationary GEV model and obtained location, scale and shape parameters from it. My understanding is, that the location and scale parameters for a non-stationary fit (with the GMST values as covariate) consist of functions that are solved for T0 and T1 respectively. I have not yet figured out how to realize this in python to calculate the probabilities p0 and p1.

I would be really grateful if someone could help me with this issue!

Datasets:

Dataset 1: daily precipitation values in mm resampled to display Three-day sums in mm only for spring months (#3,4,5) Unfortunately, I am not allowed to post my dataset, its values range from 0-45mm

Dataset 2: Yearly GMST (temperature) anomaly values in degrees Celsius Annual dataset can be downloaded here: https://www.ncei.noaa.gov/access/monitoring/global-temperature-anomalies/anomalies

What I tried:

• I computed the extremes of the dataset 1 with help of the module "pyextremes"

from pyextremes import EVA

# transform df into series
threeday_evaseries = pd.Series(threeday_springseries['threeday_prcp'])

model = EVA(threeday_evaseries)
model.get_extremes(method="BM", errors="coerce")

print(model.extremes.head())
extremes = model.extremes

• I imported the package climextremes (originally developed for R: climexRemes) to use its fit_gev function that allows for covariates

import climextremes


# dataframe conversion for R usage
import pandas as pd
import rpy2.robjects as ro
from rpy2.robjects.packages import importr
from rpy2.robjects import pandas2ri

from rpy2.robjects.conversion import localconverter

with localconverter(ro.default_converter + pandas2ri.converter):
r_from_df_extremes = ro.conversion.py2rpy(df_extremes.threeday_prcp)

r_from_df_extremes

with localconverter(ro.default_converter + pandas2ri.converter):
r_from_GMST = ro.conversion.py2rpy(GMST)

r_from_GMST


# stationary fit
result_stat = climextremes.fit_gev(r_from_df_extremes, getParams=True, bootSE=False)
result_stat

# Output:
# mle = maximum likelihood estimation
# {'mle_names': array(['location', 'scale', 'shape'], dtype='<U8'),
# 'mle': array([17.57752926, 6.09155533, 0.14391791]),
# 'se_mle_names': array(['location', 'scale', 'shape'], dtype='<U8'),
# 'se_mle': array([0.71883326, 0.56333485, 0.09658412]),
# 'nllh': array([340.05145515]),
# 'info': {'convergence': array([0], dtype=int32),
# 'counts': 'Compute Error',
# 'message': None,
# 'failure': array([0], dtype=int32)}}


# Parameters of stationary fit:
loc = 17.57752926
scale = 6.09155533
shape = 0.14391791
alpha = 0.9 --> assumed to be 0.9

• According to the exponential formula given above, I specified the following location and scale functions:

# Defining functions for non-stationary fit
locfun = loc**(alpha*r_from_GMST/loc)
scalefun = scale**(alpha*r_from_GMST/loc)

# Defining functions for non-stationary fit - year 1913
locfun_t0 = loc**(alpha*-0.35/loc)
scalefun_t0 = scale**(alpha*-0.35/loc)


For the year 1913, I inserted the GMST anomaly value which is -0.35

• For the non-stationary fit of T0 (1913), I applied the following and got an RRuntimeError

# Non-stationary fit
result_nonstat = climextremes.fit_gev(r_from_df_extremes, r_from_GMST, locationFun=locfun_t0, scaleFun=scalefun_t0, getParams=True, bootSE=False)
result_nonstat


File "C:\Users\X\anaconda3\envs\climRtest\lib\site-packages\rpy2\rinterface.py", line 810, in __call__
raise embedded.RRuntimeError(_rinterface._geterrmessage())

RRuntimeError: Error in parseParamInput(locationFun, names(x), .allowNoInt) :
parseParamInput: expecting integer-valued indices in locationFun.


What I want to achieve is this:

"[...] using a GEV distribution that scales with GMST and evaluated for the years 1900 and 2016. Reproduced from van der Wiel et al. (2017) under the licence agreement at https://creativecommons.org/licenses/by/4.0/ (last access: 31 October 2018); see original paper for more details."

I am open to suggestions for other packages or ways on how to arrive to the graphic shown above

Thank you in advance for any suggestion you might have

Continue reading...