-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Nov 12 14:02:28 2021
-
-@author: dboateng
-This routine contians all the utility classes and functions required for ESD functions
-"""
-
-import xarray as xr
-import pandas as pd
-import numpy as np
-import pickle
-import os
-from collections import OrderedDict
-import matplotlib.pyplot as plt
-
-
-
-[docs]class Dataset():
- def __init__(self, name, variables, domain_name):
- self.name = name
- self.variables = variables
- self.data = {}
- self.domain_name = domain_name
-
-
-[docs] def get(self, varname, select_domain= True, is_Dataset=False):
-
- try:
- data=self.data[varname]
-
- except KeyError:
-
- if is_Dataset == True:
- self.data[varname] = xr.open_dataset(self.variables[varname])
-
- #if self.data[varname].time[0].dt.is_month_start == False:
-
- self.data = self.data[varname]
-
- #convert date to month start
- self.data["time"] = self.data.time.values.astype("datetime64[M]")
-
- else:
- self.data[varname] = xr.open_dataarray(self.variables[varname])
-
-
- if self.data[varname].time[0].dt.is_month_start == False:
-
- #convert date to month start
- self.data["time"] = self.data.time.values.astype("datetime64[M]")
-
- data = self.data[varname]
-
-
-
- if select_domain == True:
-
- if self.domain_name == "NH":
-
- minlat, maxlat, minlon, maxlon = 10, 90, -90, 90
-
- elif self.domain_name == "Africa":
-
- minlat, maxlat, minlon, maxlon = -35, 35, -90, 90
-
-
- if hasattr(data, "longitude"):
- data = data.assign_coords({"longitude": (((data.longitude + 180) % 360) - 180)})
- data = data.where((data.latitude >= minlat) & (data.latitude <= maxlat), drop=True)
- data = data.where((data.longitude >= minlon) & (data.longitude <= maxlon), drop=True)
- else:
- data = data.assign_coords({"lon": (((data.lon + 180) % 360) - 180)})
- data = data.where((data.lat >= minlat) & (data.lat <= maxlat), drop=True)
- data = data.where((data.lon >= minlon) & (data.lon <= maxlon), drop=True)
-
- if hasattr(data, "level"):
- data = data.squeeze(dim="level")
- return data
-
- if hasattr(data, "level"):
- data = data.squeeze(dim="level")
- else:
- return data
-
-
-
-
-
-# function to estimate distance
-[docs]def haversine(lon1, lat1, lon2, lat2):
-# convert decimal degrees to radians
- lon1 = np.deg2rad(lon1)
- lon2 = np.deg2rad(lon2)
- lat1 = np.deg2rad(lat1)
- lat2 = np.deg2rad(lat2)
-
- # haversine formula
- dlon = lon2 - lon1
- dlat = lat2 - lat1
- a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
- c = 2 * np.arcsin(np.sqrt(a))
- r = 6371
- return c * r #km
-
-
-[docs]def extract_indices_around(dataset, lat, lon, radius):
- close_grids = lambda lat_, lon_: haversine(lat, lon, lat_, lon_) <= radius
-
-
- if hasattr(dataset, "longitude"):
- LON, LAT = np.meshgrid(dataset.longitude, dataset.latitude)
- else:
- LON, LAT = np.meshgrid(dataset.lon, dataset.lat)
-
- grids_index = np.where(close_grids(LAT, LON))
-
- return grids_index
-
-
-[docs]def map_to_xarray(X, datarray):
-
- coords = {}
- keys = []
-
- for k in datarray.coords.keys():
-
- if k != "time":
- coords[k] = datarray.coords.get(k)
- keys.append(k)
-
- if "expver" in keys:
- keys.remove("expver")
-
-
- if len(coords[keys[0]]) == np.shape(X)[0]:
- new = xr.DataArray(X, dims=keys, coords=coords)
- else:
- new = xr.DataArray(X, dims=keys[::-1], coords=coords)
-
- return new
-
-
-[docs]def store_pickle(stationname, varname, var, cachedir):
- filename = stationname.replace(' ', '_')
- fname = os.path.join(cachedir, filename + '_' + varname + '.pickle')
- with open(fname, 'wb') as f:
- pickle.dump(var, f)
-
-
-[docs]def load_pickle(stationname, varname, path):
- filename = stationname.replace(' ', '_')
- fname = os.path.join(path, filename + '_' + varname + '.pickle')
- with open(fname, 'rb') as f:
- return pickle.load(f)
-
-[docs]def store_csv(stationname, varname, var, cachedir):
- filename = stationname.replace(' ', '_')
- fname = os.path.join(cachedir, filename + '_' + varname + '.csv')
- var.to_csv(fname)
-
-
-[docs]def load_csv(stationname, varname, path):
- filename = stationname.replace(' ', '_')
- fname = os.path.join(path, filename + '_' + varname + '.csv')
- return pd.read_csv(fname, index_col=0, parse_dates=True)
-
-
-
-[docs]def load_all_stations(varname, path, stationnames):
- """
- This assumes that the stored quantity is a dictionary
-
- Returns a dictionary
- """
- values_dict = OrderedDict()
-
- for stationname in stationnames:
- values_dict[stationname] = load_pickle(stationname, varname, path)
-
- df = pd.DataFrame(values_dict).transpose()
-
- # get right order
- columns = list(values_dict[stationname].keys())
- df = df.loc[stationnames]
- return df[columns]
-
-
-
-
-
-
-
-
-
-
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Nov 7 17:28:48 2022
-
-@author: dboateng
-
-This module contains the regression routines. There are three layers for
-bootstrapped forward selection regression:
-
-- The ``BootstrappedRegression`` class is the outer layer. This implements the
- bootstrapping loop. This class has "regressor" member that implements the
- single regression step (i.e. a fit and a predict method). This can be the a
- ``ForwardSelection`` object, but can also be ``Lasso`` from sklearn or
- similar routines.
-- The ``ForwardSelection`` class is the next layer. This class implements a
- Forward Selection loop. This again has a regressor object that has to
- implement ``get_coefs``, ``set_coefs``, and ``average_coefs``. Additionally
- the regressor object has to implement ``fit_active``, ``fit``, and
- ``predict``.
- An example of such a regressor object is ``MultipleLSRegression``.
-
-"""
-
-from copy import copy
-import numpy as np
-import pandas as pd
-#import statsmodels.api as sm
-from sklearn.model_selection import check_cv
-from sklearn.linear_model import LinearRegression
-
-
-__all__ = ["BootstrappedRegression",
- "ForwardSelection",
- "BootstrappedForwardSelection",
- "MultipleLSRegression",
-]
-
-class MetaEstimator:
- """
- Meta estimators are classes that implement ``fit`` and ``predict``, but don't perform the
- regression themselves. They all have a member variable ``regressor`` that performs this task.
- This class contains getters and setters that are simply the ones of the regressor object.
- """
-
- def get_coefs(self):
- return self.regressor.get_coefs()
-
- def set_coefs(self, coefs):
- self.regressor.set_coefs(coefs)
- # set coef_ and intercept_ for ease of use
- if hasattr(self.regressor, 'coef_'):
- self.coef_ = self.regressor.coef_
- if hasattr(self.regressor, 'intercept_'):
- self.intercept_ = self.regressor.intercept_
- if hasattr(self.regressor, 'coefs_'):
- self.coefs_ = self.regressor.coefs_
- if hasattr(self.regressor, 'intercepts_'):
- self.intercepts_ = self.regressor.intercepts_
-
- def get_params(self):
- return self.regressor.get_params()
-
- def set_params(self, **params):
- self.regressor.set_params(**params)
-
-
-
-###############################################################################################
-# BootstrappedRegression and ForwardSelection
-###############################################################################################
-
-[docs]class BootstrappedRegression(MetaEstimator):
- """
- Performs a regression in a bootstrapping loop.
-
- This splits the data multiple times into training and test data and
- performs a regression for each split. In each loop
- the calculated parameters are stored. The final model uses the average of
- all predictors.
- If the model is a ``LinearModel`` from sklearn (i.e. it has the attributes
- ``coef_`` and ``intercept_``), the averaging routine does not have to be
- implemented. However, it can be implemented if something else than a
- arithmetic mean should be used (e.g. if only the average of robust
- predictors should be taken and everything else should be set to zero).
-
- Since this inherites from sklearn modules, it can to some extent be used
- interchangibly with other sklearn regressors.
-
-
- Parameters
- ----------
- regressor : regression object
- This should be an object similar to sklearn-like regressors that
- provides the methods ``fit(self, X_train, y_train, X_test, y_test)``
- and ``predict(self, X)``.
- This must also provide the methods ``get_coefs(self)``,
- ``set_coefs(self, coefs)``, and ``average_coefs(self, list_of_coefs)``.
- An example of this is ``ForwardSelection`` below.
- The regressor can also have a member variable ``additional_results``,
- which should be a dictionary of parameters that are calculated during
- fitting but not needed for predicting, for example metrics like the
- explained variance of predictors. In this case the regressor also needs
- the method ``average_additional_results(self, list_of_dicts)`` and
- ``set_additional_results(self, mean_additional_results)``.
- cv : integer or cross-validation generator (optional, default: None)
- This determines how the data are split:
-
- * If ``cv=None``, 3-fold cross-validation will be used.
- * If ``cv=n`` where ``n`` is an integer, n-fold cross-validation will be used.
- * If ``cv=some_object``, where ``some_object`` implements a
- ``some_object.split(X, y)`` method that returns indices for training
- and test set, this will be used. It is recommended to use
- ``YearlyBootstrapper()`` from ``stat_downscaling_tools.bootstrap``.
-
- Attributes
- ----------
- mean_coefs : type and shape depends on regressor, (only after fitting)
- Fitted coefficients (mean of all models where the coefficients
- were nonzero).
- cv_error : float (only after fitting)
- Mean of errors on test sets during bootstrapping loop.
-
- If the regressor object has the attributes ``intercept_`` and ``coef_``,
- these will also be set here.
- """
-
- def __init__(self, regressor, cv=None):
- self.regressor = regressor
- self.cv = cv
- self.has_additional_params = hasattr(self.regressor, 'additional_params')
-
-[docs] def fit(self, X, y):
- """
- Fits a model in a bootstrapping loop.
-
- Parameters
- ----------
- X : pd.DataFrame
- DataFrame of predictors
- y : pd.Series
- Series of predictand
- """
-
- #assign cross validation scheme
- cv = check_cv(self.cv, classifier=False)
-
- # number of splits in cross-validation:
- # when using cross-validaton procedures from sklearn, there is sometimes
- # no attribute cv.n_splits, therefore n_splits has to be get with the
- # built-in method cv.get_n_splits()
- if not hasattr(cv, 'n_splits'):
- n_splits = cv.get_n_splits(X)
- # otherwise n_splits is just taken directly:
- else:
- n_splits = cv.n_splits
-
- cv_error = np.zeros(n_splits)
- coefs = []
- if self.has_additional_results:
- additional_results = []
-
- for k, (train, test) in enumerate(cv.split(X, y)):
- # standardize
- X_train = X.values[train]
- X_test = X.values[test]
- y_train = y.values[train]
- y_test = y.values[test]
- # Regression
- self.regressor.fit(X_train, y_train, X_test, y_test)
- # get coefficients and error
- coefs.append(self.regressor.get_coefs())
- if self.has_additional_results:
- additional_results.append(copy(self.regressor.additional_results))
- cv_error[k] = np.sqrt(np.mean((y_test - self.regressor.predict(X_test))**2))
-
-
- # average the coefficients and the additional parameters
- self.mean_coefs = self.regressor.average_coefs(coefs)
- self.regressor.set_coefs(self.mean_coefs)
- self.cv_error = np.mean(cv_error)
- if self.has_additional_results:
- mean_add_results = self.regressor.average_additional_results(additional_results)
- self.regressor.set_additional_results(mean_add_results)
- # set coef_ and intercept_ for ease of use
- if hasattr(self.regressor, 'coef_'):
- self.coef_ = self.regressor.coef_
- if hasattr(self.regressor, 'intercept_'):
- self.intercept_ = self.regressor.intercept_
-
-
-
-[docs] def predict(self, X):
- """
- Predicts values from previously fitted coefficients.
-
- If the input X is a pandas DataFrame or Series, a Series is returned,
- otherwise only a numpy array.
-
- Parameters
- ----------
- X : pd.DataFrame
-
- Returns
- -------
- y : pd.Series
- """
- y = self.regressor.predict(X)
- if isinstance(X, pd.DataFrame) or isinstance(X, pd.Series):
- return pd.Series(data=y, index=X.index)
- else:
- return y
-
-
-
-
-[docs]class ForwardSelection(MetaEstimator):
- """
- Performs a forward selection regression.
-
- This stepwise selects the next most promising candidate predictor and adds
- it to the model if it is good enough. The method is outlined in
- "Statistical Analysis in Climate Research" (von Storch, 1999).
-
- Since this object is intended to be used in the BootstrappedRegression
- class, it implements all necessary methods.
-
- Parameters
- ----------
- regressor : regression object
- This should be an object similar to sklearn-like regressors that
- provides the methods fit and predict. Furthermore, it must also provide
- the methods ``get_coefs``, ``set_coefs``, ``average_coefs``, and ``fit_active``.
- An example of this is ``MultipleLSRegression`` below.
- min_explained_variance : float, optional (default: 0.02)
- If inclusion of the staged predictor doesn't improve the explained
- variance on the test set by at least this amount, stop the selection
- process.
-
- Attributes
- ----------
- explaned_variances : numpy array
-
-
-
- """
-
- def __init__(self, regressor, min_explained_variance=0.02):
- self.regressor = regressor
- self.min_explained_variance = min_explained_variance
- self.additional_results = {}
-
-
-[docs] def fit(self, X_train, y_train, X_test, y_test):
- """
- Cross-validated forward selection. This fits a regression model
- according to the following algorithm:
-
- 1) Start with yhat = mean(y), res = y - yhat, active = []
- 2) for each predictor in inactive set:
- - add to active set
- - perform regression
- - get error and uncertainty of error (standard deviation)
- - remove from active set
- 3) add predictor with lowest error on test set to active set
- 4) if improvement was not good enough, abort and use previous model.
-
- Parameters
- ----------
- X_train : numpy array of shape #samples x #predictors
- Array that holds the values of the predictors (columns) at
- different times (rows) for the training dataset.
- y_train : numpy array of length #samples
- Training predictand data
- X_test : numpy array of shape #samples x #predictors
- Test predictor data
- y_test : numpy array of length #samples
- Test predictand data
-
- Returns
- -------
- exp_var : numpy array of length #predictors
- explained variance of each predictor
- """
- # get some memory
- n_samples, n_predictors = X_train.shape
-
- X_train = np.array(X_train) # deep copy
- X_test = np.array(X_test) # deep copy
-
- # I store the index time series in a dictionary, so I can easily
- # remove the ones we already have and at the same time keep the number
- # of the index so I can set the correct coefficient.
- X_inactive = {i:X_train[:,i] for i in range(n_predictors)}
- active = []
-
- # initial model: no active predictor
- self.fit_active(X_train, y_train, active)
- #X_test_active = _get_active(X_test, active)
-
- residual_test = y_test - self.regressor.predict(X_test)
- SST = np.mean(residual_test**2)
- explained_variance = 0
- exp_var_predictors = np.zeros(n_predictors)
-
- error = np.zeros(n_predictors)
- old_coefs = self.regressor.get_coefs()
-
- for k in range(n_predictors):
- # perform regression with all predictors in inactive set
- inactive_mse_test = []
- inactive_mse_train = []
- inactive_coefs = []
- inactive = []
- for idx in X_inactive:
- inactive.append(idx)
- active.append(idx)
- self.fit_active(X_train, y_train, active)
- inactive_coefs.append(self.regressor.get_coefs())
- residual_test = y_test - self.regressor.predict(X_test) # turn the reset in _validate_data to False in sklean base
- residual_train = y_train - self.regressor.predict(X_train)
- inactive_mse_test.append(np.mean(residual_train**2))
- inactive_mse_train.append(np.mean(residual_train**2))
- active.pop()
- # find best predictor and add to active set/remove from inactive set
- imax = np.argmin(inactive_mse_train)
- idxmax = inactive[imax]
- SSE = inactive_mse_test[imax]
- new_explained_variance = 1 - SSE/SST
- delta_exp_var = new_explained_variance - explained_variance
- if delta_exp_var < self.min_explained_variance:
- # abort
- self.regressor.set_coefs(old_coefs)
- break
- else:
- # set the current best parameters as old parameters and start again
- old_coefs = inactive_coefs[imax]
- self.regressor.set_coefs(old_coefs)
- active.append(idxmax)
- del X_inactive[idxmax]
- exp_var_predictors[idxmax] = delta_exp_var
- explained_variance = new_explained_variance
-
- # done we only need to set the explained variance
- self.additional_results['explained variances'] = exp_var_predictors
-
-
-[docs] def fit_active(self, X, y, active):
- """
- Fits using only the columns of X whose index is in ``active``.
- """
- X_active = _get_active(X, active)
- self.regressor.fit(X_active, y)
- self.regressor.set_expand_coefs(active, X.shape[1])
-
-
-
-[docs] def predict_active(self, X, active):
- X_active = _get_active(X, active)
- return self.regressor.predict(X_active)
-
-[docs] def set_additional_results(self, add_results):
- self.additional_results = copy(add_results)
- self.explained_variances = add_results['explained variances']
-
-[docs] def average_additional_results(self, list_of_params):
- n_params = len(list_of_params)
- n_predictors = len(list_of_params[0]['explained variances'])
- exp_var = np.zeros((n_params, n_predictors))
- for i, p in enumerate(list_of_params):
- exp_var[i,:] = p['explained variances']
- add_results = {'explained variances':robust_average(exp_var)}
- return add_results
-
-
-
-###############################################################################################
-# Regressors for Forward Selection
-###############################################################################################
-
-
-class LinearCoefsHandlerMixin:
-
- def get_coefs(self):
- """
- Returns all fitted coefficients in the same order as ``set_coefs`` needs them.
- """
- coefs = np.zeros(len(self.coef_) + 1)
- coefs[0:-1] = self.coef_
- coefs[-1] = self.intercept_
- return coefs
-
- def set_coefs(self, coefs):
- self.coef_ = coefs[0:-1]
- self.intercept_ = coefs[-1]
-
- def average_coefs(self, list_of_coefs):
- """
- Calculates the average of robust predictors, i.e. of those that are
- nonzero in at least 50% of the cases.
- """
- coefs = np.asarray(list_of_coefs)
- return robust_average(coefs)
-
- def get_params(self):
- return {}
-
- def set_params(self, **params):
- pass
-
- def fit_active(self, X, y, active):
- """
- Fits using only the columns of X whose index is in ``active``.
- """
- X_active = _get_active(X, active)
- self.fit(X_active, y)
- self.set_expand_coefs(active, X.shape[1])
-
- def predict_active(self, X, active):
- """
- Predict using only the columns of X whose index is in ``active``.
- """
- X_active = _get_active(X, active)
- self.predict(X_active)
-
-
-[docs]class MultipleLSRegression(LinearCoefsHandlerMixin):
- """
- Implementation of multiple linear OLS regression to be used with
- ForwardSelection and BootstrappedRegression.
- The following methods are implemented:
-
- - ``fit``
- - ``predict``
- - ``get_coefs``
- - ``set_coefs``
- - ``average_coefs``
- - ``fit_active``
- """
-
- def __init__(self):
- self.lm = LinearRegression()
-
-[docs] def fit(self, X, y):
- if X.shape[1] == 0:
- # only intercept model
- self.intercept_ = np.mean(y)
- else:
- self.lm.fit(X, y)
-
-[docs] def predict(self, X):
- n, m = X.shape
- self.lm.coef_ = self.coef_
- self.lm.intercept_ = self.intercept_
- return self.lm.predict(X)
-
-[docs] def set_expand_coefs(self, active, n_predictors):
- """
- This will be called after ``fit``, since fit will often be called with only some of the
- predictors. This expands the current coefficients and expands them in a way such that
- ``predict`` can be called with all predictors.
- """
- # get a full coefficient vector back
- coefs = np.zeros(n_predictors)
- for i, idx in enumerate(active):
- coefs[idx] = self.lm.coef_[i]
- self.coef_ = coefs
-
- # get intercept_: if len(active) == 0, the intercept was already set
- if len(active) != 0:
- self.intercept_ = self.lm.intercept_
-
-
-# class GammaRegression(LinearCoefsHandlerMixin):
-# """
-# Implementation of generalized linear gamma regression.
-# """
-
-# def __init__(self, family=sm.families.Gamma()):
-# self.family = family
-
-# def fit(self, X, y):
-# n, m = X.shape
-# G = np.zeros((n, m+1))
-# G[:,0:-1] = X
-# G[:,-1] = np.ones(n)
-# self.glm = sm.GLM(y, G, family=self.family)
-# gamma_results = self.glm.fit()
-# self.coef_ = gamma_results.params[0:-1]
-# self.intercept_ = gamma_results.params[-1]
-
-# def predict(self, X):
-# n, m = X.shape
-# G = np.zeros((n, m+1))
-# G[:,0:-1] = X
-# G[:,-1] = np.ones(n)
-# params = np.zeros(len(self.coef_) + 1)
-# params[0:-1] = self.coef_
-# params[-1] = self.intercept_
-# return self.glm.predict(params, exog=G)
-
-
-# def set_expand_coefs(self, active, n_predictors):
-# # get a full coefficient vector back
-# coefs = np.zeros(n_predictors)
-# for i, idx in enumerate(active):
-# coefs[idx] = self.coef_[i]
-# self.coef_ = coefs
-
-
-
-###############################################################################################
-# Some functions
-###############################################################################################
-
-[docs]def _get_active(X, active):
- """
- Returns a new matrix X_active with only the columns of X whose index is in ``active``.
- """
- Xnew = np.empty((X.shape[0], len(active)))
- for i, idx in enumerate(active):
- Xnew[:,i] = X[:,idx]
- return Xnew
-
-def robust_average(coefs):
- """
- Takes the robust average of a coefficient matrix.
-
- Parameters
- ----------
- coefs : numpy 2d-array, n_coefs x n_predictors
-
- Returns
- -------
- mean_coefs : numpy array, length n_predictors
- """
- n_coefs, n_predictors = coefs.shape
- mean_coefs = np.zeros(n_predictors)
- for i in range(n_predictors):
- c = coefs[:,i] # one column
- if len(c[c != 0]) > 0.5*len(c):
- mean_coefs[i] = np.mean(c)
- return mean_coefs
-
-
-
-###############################################################################################
-# Easy to use classes
-###############################################################################################
-
-[docs]class BootstrappedForwardSelection(BootstrappedRegression):
- """
- This is an easy to use interface for BootstrappedRegression with ForwardSelection.
-
- Parameters
- ----------
- regressor : regression object
- This should be an object similar to sklearn-like regressors that
- provides the methods fit and predict. Furthermore, it must also provide
- the methods ``get_coefs``, ``set_coefs``, ``average_coefs``, and ``fit_active``.
- An example of this is ``MultipleLSRegression`` below.
- min_explained_variance : float, optional (default: 0.02)
- If inclusion of the staged predictor doesn't improve the explained
- variance on the test set by at least this amount, stop the selection
- process.
- cv : integer or cross-validation generator (optional, default: None)
- This determines how the data are split:
-
- * If ``cv=None``, 3-fold cross-validation will be used.
- * If ``cv=n`` where ``n`` is an integer, n-fold cross-validation will be used.
- * If ``cv=some_object``, where ``some_object`` implements a
- ``some_object.split(X, y)`` method that returns indices for training
- and test set, this will be used. It is recommended to use
- ``YearlyBootstrapper()`` from ``stat_downscaling_tools.bootstrap``.
- """
-
- def __init__(self, regressor, min_explained_variance=0.02, cv=None):
- self.regressor = ForwardSelection(regressor, min_explained_variance)
- self.cv = cv
- self.has_additional_results = True
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Nov 12 14:02:45 2021
-
-@author: dboateng
-"""
-
-# importing modules
-from abc import ABC, abstractmethod
-import os
-import sys
-import pickle
-import numpy as np
-import pandas as pd
-
-def _check_data_available(data, daterange):
- for d in daterange:
- data.loc[d]
-
-
-[docs]class Predictor(ABC):
-
- def __init__(self, name, longname=None, cachedir=None):
-
- self.name = name
- if longname is not None:
- self.longname = longname
- else:
- self.longname = name
-
- if cachedir is not None:
- self.cachedir = cachedir
- else:
- self.cachedir = '.predictors'
-
- self.data = {}
- self.params = {}
- self.patterns = {}
-
-
-[docs] def save(self):
-
- if not os.path.isdir(self.cachedir):
- try:
- os.makedirs(self.cachedir)
- except:
- print("There might be problem making the directory" + self.cachedir +
- "which is required to store predictors", file=sys.stderr)
- raise
-
-
- filename_to_store = os.path.join(self.cachedir, self.longname + ".pickle")
-
- with open(filename_to_store, "wb") as f:
- predictordata = {"data":self.data, "params":self.params, "patterns":self.patterns}
-
- # serialize the the predictor data with dump()
- pickle.dump(predictordata, f)
-
-[docs] def load(self):
- filename_to_store = os.path.join(self.cachedir, self.longname + ".pickle")
-
- if not os.path.exists(filename_to_store):
- raise FileNotFoundError("Predictor data may not be available in serialize form")
-
-
- with open(filename_to_store, "rb") as f:
-
- predictordata = pickle.load(f)
- self.data = predictordata["data"]
- self.params = predictordata["params"]
- self.patterns = predictordata["patterns"]
-
-
-[docs] def get(self, daterange, dataset, fit, regenerate=False, patterns_from=None, params_from=None):
-
- if patterns_from is None:
- patterns_from = dataset.name
- if params_from is None:
- params_from = dataset.name
-
- data_key = "data=" + dataset.name + "_patterns=" +params_from + "_params=" + params_from
-
- if not self.data and not self.params and not regenerate:
- try:
- self.load()
- except FileNotFoundError:
- pass
- if not self.data:
- regenerate = True
-
- if not regenerate:
- try:
- _check_data_available(self.data[data_key], daterange)
- except KeyError:
- regenerate=True
-
- if regenerate:
- print("Regenerating predictor data for", self.name, "using dataset", dataset.name,
- "with loading patterns and params from", patterns_from, "and", params_from)
-
-
- if dataset.name not in self.params:
- self.params[dataset.name] = {}
- data = self._generate(daterange, dataset, fit, patterns_from, params_from)
-
- if data_key in self.data:
- self.data[data_key] = self.data[data_key].combine_first(data)
- else:
- self.data[data_key] = data
-
- self.save()
-
- data = self.data[data_key].loc[daterange]
-
- try:
- _check_data_available(data, daterange)
- except KeyError:
- print("Predictor data for", self.name, "could not be generated for all required timesteps",
- file=sys.stderr)
- raise
- return data
-
-
-
-
- @abstractmethod
- def _generate(self, daterange, dataset, fit, patterns_from, params_from):
- ...
-
-
-
-
-[docs] def plot(self, daterange, dataset, fit, regenerate=False, patterns_from=None, params_from=None,
- **plot_kwargs):
- data = self.get(daterange, dataset, fit, regenerate, patterns_from, params_from)
- handle = data.plot(**plot_kwargs)
- return handle
-
-
-
-
-
-
-
-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Nov 12 14:03:09 2021
-
-@author: dboateng
-"""
-# importing modules
-
-import numpy as np
-import xarray as xr
-import pandas as pd
-
-
-
-
-try:
- from .Predictor_Base import Predictor
- from .ESD_utils import extract_indices_around
-except:
- from Predictor_Base import Predictor
- from ESD_utils import extract_indices_around
-
-
-[docs]class RegionalAverage(Predictor):
-
- def __init__(self, name, lat, lon, standardizer_constructor=None,
- radius=250, **kwargs):
-
- self.lon = lon
- self.lat = lat
- self.varname = name
- self.radius = radius
- self.standardizer_constructor = standardizer_constructor
-
- super().__init__(name, name + '_{:2.2f}N_{:2.2f}E'.format(self.lat, self.lon),
- **kwargs)
-
-
- def _generate(self, daterange, dataset, fit, patterns_from, params_from):
-
- params = self.params[params_from]
-
- da = dataset.get(self.varname)
-
- da = da.sel(time=daterange)
-
- #da = da[self.varname]
-
- if hasattr(da, "longitude"):
- da = da.rename({"longitude":"lon", "latitude":"lat"})
-
- if "indices" not in params or fit is True:
- params["indices"] = extract_indices_around(da, self.lat, self.lon, self.radius)
- values = da.isel(lat=params["indices"][0], lon=params["indices"][1])
-
- data = values.mean(dim=("lat", "lon")).to_series().astype(np.double)
-
- if self.standardizer_constructor is not None:
- if fit:
- params["standardizer"] = self.standardizer_constructor()
- params["standardizer"].fit(data)
- data = params["standardizer"].transform(data)
-
- return pd.Series(data)
-
-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Sun Nov 21 00:55:37 2021
-
-@author: dboateng
-"""
-
-import pickle
-import numpy as np
-import pandas as pd
-import os
-
-
-
-try:
- from .Predictor_Generator import *
- from .standardizer import MonthlyStandardizer
- from .predictand import PredictandTimeseries
- from .teleconnections import NAO, SCAN, EA, EAWR
-except:
- from Predictor_Generator import *
- from standardizer import MonthlyStandardizer
- from predictand import PredictandTimeseries
- from teleconnections import NAO, SCAN, EA, EAWR
-
-
-[docs]class StationOperator():
-
- def __init__(self, data, name, lat, lon, elevation):
-
- self.variables = {}
- for varname in data:
- self.variables[varname] = PredictandTimeseries(data[varname])
- self.name = name
- self.lat = lat
- self.lon = lon
- self.elevation = elevation
- print(name, lat, lon, elevation)
-
-[docs] def get_var(self, variable, daterange, anomalies=True):
-
- y = self.variables[variable].get(daterange, anomalies=anomalies)
-
- return y
-
-[docs] def set_predictors(self, variable, predictors, cachedir, radius=250, detrending=False, scaling=False,
- standardizer=None):
- predictor_list = []
-
- for name in predictors:
- if name == "NAO":
- predictor_list.append(NAO(cachedir=cachedir))
- elif name == "EA":
- predictor_list.append(EA(cachedir=cachedir))
-
- elif name == "SCAN":
- predictor_list.append(SCAN(cachedir=cachedir))
-
- elif name == "EAWR":
- predictor_list.append(EAWR(cachedir=cachedir))
- else:
-
- if standardizer == None:
- predictor_list.append(RegionalAverage(name, self.lat, self.lon, radius=radius, cachedir=cachedir,
- standardizer_constructor=lambda:
- MonthlyStandardizer(detrending=detrending, scaling=scaling)))
- else:
- predictor_list.append(RegionalAverage(name, self.lat, self.lon, radius=radius, cachedir=cachedir,
- standardizer_constructor=lambda:
- standardizer))
-
- self.variables[variable].set_predictors(predictor_list)
-
-[docs] def set_transform(self, variable, transform):
-
- self.variables[variable].set_transform(transform)
-
-[docs] def set_standardizer(self, variable, standardizer):
-
- self.variables[variable].set_standardizer(standardizer)
-
-
-
-[docs] def set_model(self, variable, method, ensemble_learning=False, estimators=None, cv=10, final_estimator_name=None,
- daterange =None, predictor_dataset=None, fit_predictors=True,
- scoring=["r2", "neg_root_mean_squared_error"], **predictor_kwargs):
-
-
- self.variables[variable].set_model(method, ensemble_learning=ensemble_learning, estimators=estimators, cv=cv, final_estimator_name=final_estimator_name,
- daterange =daterange , predictor_dataset=predictor_dataset,
- fit_predictors=fit_predictors, scoring=scoring, **predictor_kwargs)
-
-
-
- def _get_predictor_data(self,variable, daterange , dataset, fit_predictors=True, **predictor_kwargs):
-
- return self.variables[variable]._get_predictor_data(daterange , dataset, fit_predictors=fit_predictors, **predictor_kwargs)
-
-
-
-[docs] def predictor_correlation(self, variable, daterange, predictor_dataset, fit_predictors=True, fit_predictand=True,
- method="pearson", **predictor_kwargs):
-
- return self.variables[variable].predictor_correlation(daterange, predictor_dataset, fit_predictors=fit_predictors,
- fit_predictand=fit_predictand,
- method= method, **predictor_kwargs)
-
-
-
-[docs] def fit_predictor(self, variable, name, daterange, predictor_dataset):
-
- self.variables[variable].fit_predictor(name, daterange, predictor_dataset)
-
-
-
-[docs] def fit(self, variable, daterange , predictor_dataset, fit_predictors=True , predictor_selector=True, selector_method="Recursive",
-
- selector_regressor="Ridge", num_predictors=None, selector_direction=None, cal_relative_importance=False,
- fit_predictand=True, impute=False, impute_method=None,
- impute_order=None,**predictor_kwargs):
-
-
-
- return self.variables[variable].fit(daterange , predictor_dataset, fit_predictors=fit_predictors , predictor_selector=predictor_selector,
- selector_method=selector_method,
- selector_regressor= selector_regressor,
- num_predictors=num_predictors,
- selector_direction= selector_direction,
- cal_relative_importance = cal_relative_importance,
- fit_predictand = fit_predictand,
- impute=impute, impute_method=impute_method,
- impute_order= impute_order,
- **predictor_kwargs)
-
-
-[docs] def predict(self, variable, daterange , predictor_dataset, fit_predictand=True, fit_predictors=True, **predictor_kwargs):
-
- return self.variables[variable].predict(daterange , predictor_dataset, fit_predictand=fit_predictand,
- fit_predictors=fit_predictors,
- **predictor_kwargs)
-
-
-[docs] def cross_validate_and_predict(self, variable, daterange , predictor_dataset, fit_predictand=True,
- return_cv_scores=False, **predictor_kwargs):
-
- return self.variables[variable].cross_validate_and_predict(daterange , predictor_dataset,
- fit_predictand=fit_predictand,
- return_cv_scores=return_cv_scores,
- **predictor_kwargs)
-
-[docs] def evaluate(self, variable, daterange, predictor_dataset, fit_predictand=True, **predictor_kwargs):
-
- return self.variables[variable].evaluate(daterange, predictor_dataset, fit_predictand=fit_predictand, **predictor_kwargs)
-
-
-[docs] def ensemble_transform(self, variable, daterange, predictor_dataset, **predictor_kwargs):
-
- return self.variables[variable].ensemble_transform(daterange, predictor_dataset, **predictor_kwargs)
-
-[docs] def relative_predictor_importance(self, variable):
-
- return self.variables[variable].relative_predictor_importance()
-
-
-
-
-[docs] def tree_based_feature_importance(self, variable, daterange, predictor_dataset, fit_predictand=True,
- plot=False, **predictor_kwargs):
-
- return self.variables[variable].tree_based_feature_importance(daterange, predictor_dataset, fit_predictand=fit_predictand,
- plot=plot, **predictor_kwargs)
-
-
-[docs] def tree_based_feature_permutation_importance(self, variable, daterange, predictor_dataset, fit_predictand=True,
- plot=False, **predictor_kwargs):
-
- return self.variables[variable].tree_based_feature_permutation_importance(daterange, predictor_dataset, fit_predictand=fit_predictand,
- plot=plot, **predictor_kwargs)
-
-
-
-[docs] def climate_score(self, variable, fit_period, score_period, predictor_dataset,
- **predictor_kwargs):
- """
- Calculate the climate score of a fitted model for the given variable.
-
- Parameters
- ----------
- variable : string
- Variable name. "Temperature" or "Precipitation"
- fit_period : pd.DatetimeIndex
- Range of data that should will be used for creating the reference prediction.
- score_period : pd.DatetimeIndex
- Range of data for that the prediction score is evaluated
- predictor_dataset : stat_downscaling_tools.Dataset
- The dataset that should be used to calculate the predictors
- predictor_kwargs : keyword arguments
- These arguments are passed to the predictor's get function
-
- Returns
- -------
- cscore : double
- Climate score (similar to rho squared). 1 for perfect fit, 0 for no
- skill, negative for even worse skill than mean prediction.
- """
- return self.variables[variable].climate_score(fit_period, score_period, predictor_dataset,
- **predictor_kwargs)
-
-
-[docs] def get_explained_variance(self, variable):
- """
- If the model is fitted and has the attribute ``explained_variance``,
- returns it, otherwise returns an array of zeros.
- """
- return self.variables[variable].explained_variance_predictors
-
-
-[docs] def save(self, directory=None, fname=None):
- """
- Saves the weatherstation object to a file (pickle).
-
- Parameters
- ----------
- directory : str, optional (default : None)
- Directory name where the pickle-file should be stored. Defaults to
- the current directory.
- fname : str, optional (default: None)
- Filename of the file where the station should be stored. Defaults
- to ``self.name.replace(' ', '_') + '.pickle'``.
- """
- if directory is None:
- directory = './'
- if fname is None:
- fname = self.name.replace(' ', '_') + '.pickle'
- filename = os.path.join(directory, fname)
-
- with open(filename, 'wb') as f:
- pickle.dump(self, f)
-
-
-[docs]def load_station(fname):
- """
- Loads a pickled station from the given file
- """
- with open(fname, 'rb') as f:
- so = pickle.load(f)
- return so
-
-
-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Nov 12 14:01:43 2021
-
-
-This routine handles the preprocessing of data downloaded directly from DWD. The default time series is monthly, others frequency must be pass to the function
-1. Extracting only stations with required number of years
-2. Writing additional information into files (eg. station name, lat, lon and elevation), since its downloaded into a separate file using station codes
-3. All utils function to read stations into pyESD Station operator class
-
-Note: This routine is specifically designed for data downloded from DWD (otherwise please contact daniel.boateng@uni-tuebingen.de for assistance on other datasets)
-
-@author: dboateng
-"""
-
-# importing packages
-import os
-import pandas as pd
-import numpy as np
-
-#from StationOperator import StationOperator
-# #local packages
-try:
- from .StationOperator import StationOperator
-except:
- from StationOperator import StationOperator
-
-
-
-
-[docs]def read_weatherstationnames(path_to_data):
- """
- This function reads all the station names in the data directory
-
- Parameters
- ----------
- path_to_data : TYPE: str
- DESCRIPTION. The directory path to where all the station data are stored
-
- Returns
- -------
- namedict : TYPE: dict
- DESCRIPTION.
-
- """
-
- nr, name = np.loadtxt(os.path.join(path_to_data, 'stationnames.csv'),
- delimiter=',', skiprows=1, usecols=(0,1),
- dtype=str, unpack=True)
-
- nr = [int(i) for i in nr]
- namedict = dict(zip(nr, name))
- return namedict
-
-
-[docs]def read_station_csv(filename, varname, return_all=False):
- """
-
-
- Parameters
- ----------
- filename : TYPE: str
- DESCRIPTION. Name of the station in path
- varname : TYPE: str
- DESCRIPTION. The name of the varibale to downscale (eg. Precipitation, Temperature)
-
- Raises
- ------
- ValueError
- DESCRIPTION.
-
- Returns
- -------
- ws : TYPE
- DESCRIPTION.
-
- """
-
- # reading headers info with readline
- with open(filename, "r") as f:
- name = f.readline().split(',')[1].replace("\n", "")
- lat = float(f.readline().split(',')[1].replace("\n",""))
- lon = float(f.readline().split(',')[1].replace("\n",""))
- elev = float(f.readline().split(',')[1].replace("\n",""))
-
-
-
- data = pd.read_csv(filename, sep=',', skiprows=6, usecols=[0,1,],
- parse_dates=[0], index_col=0, names=['Time',
- varname])
-
- data = data.dropna()
-
- if varname == "Precipitation":
- pr = data[varname]
- pr[pr == -9999] = np.nan
- assert not np.any(pr < -1e-2)
- assert not np.any(pr > 2000)
-
- data = {varname:pr}
-
- elif varname == "Temperature":
- t = data[varname]
- t[t == -8888] = np.nan
- assert not np.any(t < -50)
- assert not np.any(t > 80)
-
- data = {varname:t}
- else:
- raise ValueError("The model does not recognize the variable name")
-
- if return_all == False:
-
- so = StationOperator(data, name, lat, lon, elev)
-
- return so
-
- else:
- return data, lat, lon
-
-
-
-
-[docs]def read_weatherstations(path_to_data):
- """
- Read all the station data in a directory.
-
- Parameters
- ----------
- path_to_data : TYPE: STR
- DESCRIPTION. relative or absolute path to the station folder
-
- Returns
- -------
- stations : TYPE: DICT
- DESCRIPTION. Dictionary containing all the datasets
-
- """
- namedict = read_weatherstationnames(path_to_data)
- stations = {}
- for i in namedict:
- filename = namedict[i].replace(' ', '_') + '.csv'
- print("Reading", filename)
- ws = read_station_csv(os.path.join(path_to_data, filename))
- stations[i] = ws
- return stations
-
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Mar 14 11:02:35 2022
-
-@author: dboateng
-"""
-
-# importing models
-import sys
-from sklearn.ensemble import StackingRegressor, VotingRegressor
-from sklearn.model_selection import cross_val_score, cross_validate, cross_val_predict
-
-try:
- from models import Regressors
-except:
- from .models import Regressors
-
-
-[docs]class EnsembleRegressor():
-
- def __init__(self, estimators, final_estimator_name=None, cv=10, n_jobs=-1, passthrough=False, method="Stacking",
- scoring=None):
- self.estimators = estimators
- self.final_estimator_name = final_estimator_name
- self.cv = cv
- self.n_jobs = n_jobs
- self.passthrough = passthrough
- self.method = method
- self.scoring = scoring
-
-
- if self.final_estimator_name == None:
- regressor = Regressors(method="RandomForest", cv=10)
- regressor.set_model()
- print("....using random forest as final estimator......")
- self.final_estimator = regressor.estimator
- else:
-
- regressor = Regressors(method=final_estimator_name, cv=10)
- regressor.set_model()
- print("-------using " + final_estimator_name + " as final estimator ......")
- self.final_estimator = regressor.estimator
-
-
-
- if self.method == "Stacking":
- self.ensemble = StackingRegressor(estimators=self.estimators, final_estimator=self.final_estimator, cv=self.cv,
- passthrough=self.passthrough, n_jobs=self.n_jobs)
- elif self.method == "Voting":
- self.ensemble = VotingRegressor(estimators=self.estimators, n_jobs=self.n_jobs)
-
- else:
- raise ValueError("The ensembles are Stacking Generalization or Voting, check the name of the method")
-
-
- if scoring == None:
- self.scoring= ["r2", "neg_root_mean_squared_error"]
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-[docs] def predict_average(self, X):
- y_preds = self.transform(X)
- y_avg = y_preds.mean(axis=1)
-
- return y_avg
-
-[docs] def cross_val_score(self, X, y):
- return cross_val_score(self.ensemble, X, y, cv=self.cv, scoring=self.scoring)
-
-[docs] def cross_validate(self, X, y):
- return cross_validate(self.ensemble, X, y, scoring=self.scoring,
- n_jobs=2, verbose=0, cv=self.cv)
-
-[docs] def cross_val_predict(self, X, y):
- return cross_val_predict(self.ensemble, X, y, n_jobs=2, verbose=0)
-
-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Jan 3 17:18:14 2022
-
-@author: dboateng
-"""
-
-# importing modules
-import pandas as pd
-import numpy as np
-import matplotlib.pyplot as plt
-from sklearn.feature_selection import RFECV, SelectFromModel, SequentialFeatureSelector
-from sklearn.linear_model import Lasso, LassoCV, Ridge, BayesianRidge, ARDRegression
-from sklearn.inspection import permutation_importance
-from sklearn.ensemble import RandomForestRegressor
-from sklearn.model_selection import TimeSeriesSplit, LeaveOneOut, LeaveOneGroupOut
-from sklearn.svm import SVR
-from sklearn.ensemble import ExtraTreesRegressor, ExtraTreesClassifier
-
-
-
-# class for recursive feature elimination
-[docs]class RecursiveFeatureElimination():
-
- def __init__(self, regressor_name="ARD"):
- self.regressor_name = regressor_name
- self.cv = TimeSeriesSplit()
- self.n_jobs = -1
- self.min_features = 5
- self.scoring = "r2"
-
-
- if self.regressor_name == "ARD":
- self.estimator = ARDRegression()
- elif self.regressor_name == "BayesianRidge":
- self.estimator = BayesianRidge()
- elif self.regressor_name == "lasso":
- self.estimator = Lasso()
- elif self.regressor_name == "lassocv":
- self.estimator = LassoCV()
- elif self.regressor_name == "Ridge":
- self.estimator = Ridge()
- elif self.regressor_name == "RandomForest":
- self.estimator = RandomForestRegressor()
- else:
- raise ValueError("Check the regressor if implemented")
-
-
-[docs] def fit(self, X, y):
-
- self.regressor = RFECV(estimator=self.estimator, scoring=self.scoring, cv=self.cv, n_jobs= self.n_jobs,
- min_features_to_select=self.min_features).fit(X, y)
-
-
-[docs] def print_selected_features(self, X):
- num_features = self.regressor.n_features_
- select_names = X.columns[self.regressor.support_]
-
- self.select_names = select_names
-
- print("{0} : optimal number of predictors and selected variables are {1}".format(num_features, select_names))
-
-
-
-[docs] def cv_test_score(self):
- cv_score = self.regressor.cv_results_["mean_test_score"].mean()
- return cv_score
-
-
-
-
-
-[docs]class TreeBasedSelection():
- def __init__(self, regressor_name="RandomForest"):
- self.regressor_name = regressor_name
- self.n_jobs = -1
- self.bootstrap = True
- self.criterion = "squared_error"
- self.scoring = "r2"
-
- if self.regressor_name == "RandomForest":
- self.estimator = RandomForestRegressor(n_jobs=self.n_jobs, criterion=self.criterion, bootstrap=self.bootstrap,
- n_estimators=200)
- elif self.regressor_name == "ExtraTree":
- self.estimator = ExtraTreesRegressor(n_jobs=self.n_jobs, criterion=self.criterion, bootstrap=self.bootstrap)
- else:
- raise ValueError("Tree regressor estimator is not defined properly")
-
-
-[docs] def fit(self, X, y):
-
- self.regressor = SelectFromModel(estimator=self.estimator, prefit=False).fit(X,y)
-
-
-
-[docs] def feature_importance(self, X,y, plot=False):
- self.estimator.fit(X,y)
- importance = self.estimator.feature_importances_
- feature_names = X.columns
- forest_importances = pd.Series(importance, index=feature_names)
-
- if plot == True:
- std = np.std([tree.feature_importances_ for tree in self.estimator.estimators_], axis=0)
- fig,ax = plt.subplots()
- forest_importances.plot.bar(yerr=std, ax=ax)
- ax.set_title("Feature importances using tree regressor")
- fig.tight_layout()
- plt.show()
-
- return forest_importances
-
-[docs] def permutation_importance_(self, X,y, plot=False):
- self.estimator.fit(X,y)
- importance = permutation_importance(estimator=self.estimator, X=X, y=y, scoring=self.scoring,
- n_repeats=10, n_jobs=self.n_jobs)
- sorted_idx = importance.importance_mean.argsort()
- if plot == True:
- fig,ax = plt.subplots()
- ax.boxplot(importance.importances[sorted_idx].T, vert=False, labels=X.columns[sorted_idx])
- ax.set_title("Permutation Importances (On test data)")
- fig.tight_layout()
- plt.show()
-
- return sorted_idx
-
-
-[docs] def print_selected_features(self, X):
- select_names = X.columns[self.regressor.get_support(indices=True)]
- num_features = len(select_names)
- print("{0} : optimal number of predictors and selected variables are {1}".format(num_features, select_names))
-
- self.select_names = select_names
-
-
-[docs]class SequentialFeatureSelection():
- def __init__(self, regressor_name = "Ridge", n_features= 10, direction="forward"):
- self.regressor_name = regressor_name
- self.n_features =n_features
- self.scoring = "r2"
- self.direction = direction
-
- if self.regressor_name == "ARD":
- self.estimator = ARDRegression()
- elif self.regressor_name == "BayesianRidge":
- self.estimator = BayesianRidge()
- elif self.regressor_name == "lasso":
- self.estimator = Lasso()
- elif self.regressor_name == "lassocv":
- self.estimator = LassoCV()
- elif self.regressor_name == "Ridge":
- self.estimator = Ridge()
- else:
- raise ValueError("Check the regressor if implemented")
-
-[docs] def fit(self, X,y):
- self.regressor = SequentialFeatureSelector(estimator=self.estimator, n_features_to_select=self.n_features, scoring=self.scoring,
- direction=self.direction).fit(X,y)
-
-
-
-
-
-[docs] def print_selected_features(self, X):
- select_names = X.columns[self.regressor.get_support(indices=True)]
- num_features = len(select_names)
- print("{0} : optimal number of predictors and selected variables are {1}".format(num_features, select_names))
-
-
- self.select_names = select_names
-
-
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Mar 16 11:34:25 2022
-
-@author: dboateng
-"""
-
-#importing models
-
-import numpy as np
-
-from sklearn.metrics import mean_squared_error, accuracy_score, r2_score, explained_variance_score, max_error
-from sklearn.metrics import mean_absolute_error, mean_squared_log_error
-
-
-[docs]class Evaluate():
-
- def __init__(self, y_true, y_pred):
- self.y_true = y_true
- self.y_pred = y_pred
-
-[docs] def RMSE(self):
- error = (self.y_true - self.y_pred).dropna()
- score = np.sqrt(np.mean(error **2))
- print("RMSE: {:2f}".format(score))
-
- return score
-
-[docs] def NSE(self):
- score = (1-(np.sum((self.y_pred - self.y_true)**2))/np.sum((self.y_true - np.mean(self.y_true))**2))
-
- print("Nach-Sutcliffe Efficiency(NSE): {:2f}".format(score))
-
- return score
-
-[docs] def MSE(self):
- score = mean_squared_error(self.y_true, self.y_pred, squared=False)
- print("Mean Squared Error): {:.2f}".format(score))
-
- return score
-
-
-[docs] def MAE(self):
- score = mean_absolute_error(self.y_true, self.y_pred)
- print("Mean Absolute Error): {:.2f}".format(score))
-
- return score
-
-
-
-[docs] def R2_score(self):
- score = r2_score(self.y_true, self.y_pred)
- print("R² (Coefficient of determinaiton): {:.2f}".format(score))
-
- return score
-
-[docs] def explained_variance(self):
-
- score = explained_variance_score(self.y_true, self.y_pred)
- print("Explained Variance: {:.2f}".format(score))
-
- return score
-
-
-
-[docs] def max_error(self):
-
- score = max_error(self.y_true, self.y_pred)
- print("Maximum error: {:2f}".format(score))
-
- return score
-
-[docs] def adjusted_r2(self):
-
- r2 = r2_score(self.y_true, self.y_pred)
-
- adj_r2 = (1 - (1- r2) * ((self.y_true.shape[0]- 1) /
- (self.y_true.shape[0] - self.y_true[1] -1)))
-
- print("Adjusted R²: {:.2f}".format(adj_r2))
-
- return adj_r2
-
-
-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Jan 25 16:00:11 2022
-
-@author: dboateng
-
-"""
-
-# imporitng modules
-import numpy as np
-import pandas as pd
-import sys
-
-
-#from sklearn
-from sklearn.linear_model import LassoCV, LassoLarsCV
-from sklearn.model_selection import GridSearchCV, RandomizedSearchCV
-from sklearn.linear_model import ARDRegression, BayesianRidge
-from sklearn.linear_model import GammaRegressor, PoissonRegressor
-from sklearn.linear_model import SGDRegressor, RidgeCV
-from sklearn.neural_network import MLPRegressor
-from sklearn.svm import SVR
-from sklearn.ensemble import RandomForestRegressor, ExtraTreesRegressor, BaggingRegressor, AdaBoostRegressor
-from sklearn.ensemble import HistGradientBoostingRegressor, GradientBoostingRegressor
-from sklearn.model_selection import cross_val_score, cross_validate, cross_val_predict
-
-from xgboost import XGBRegressor
-
-from skopt import BayesSearchCV
-from skopt.space import Real, Categorical, Integer
-
-
-try:
- from .splitter import MonthlyBooststrapper
-
-except:
- from splitter import MonthlyBooststrapper
-
-
-
-
-
-
-[docs]class HyperparameterOptimize(MetaAttributes):
-
- def __init__(self, method, param_grid, regressor, scoring="r2", cv=10):
- self.method = method
- self.param_grid = param_grid
- self.scoring = scoring
- self.cv = cv
- self.regressor = regressor
-
- if self.method == "GridSearchCV":
- self.hyper = GridSearchCV(estimator= self.regressor, param_grid=self.param_grid,
- scoring=self.scoring, cv=self.cv)
-
- elif self.method == "RandomizedSearchCV":
- self.hyper = RandomizedSearchCV(estimator= self.regressor, param_distributions=self.param_grid,
- scoring=self.scoring, cv=self.cv)
- elif self.method == "BayesSearchCV":
- self.hyper = BayesSearchCV(estimator=self.regressor, search_spaces=self.param_grid, scoring=self.scoring,
- cv=self.cv,)
- else:
- raise ValueError("The defined hyperparameter search can't be recognized")
-
-[docs] def fit(self, X, y):
- self.hyper = self.hyper.fit(X,y)
- self.estimator = self.hyper.best_estimator_
- return self.estimator
-
-
-
-
-
-
-
-
-
-
-
-[docs] def cross_validate(self, X, y):
- return cross_validate(self.hyper, X, y, scoring=["r2", "neg_root_mean_squared_error"],
- n_jobs=2, verbose=0, cv=self.cv)
-
-
-[docs] def cross_val_predict(self, X, y):
- return cross_val_predict(self.estimator, X, y, n_jobs=2, verbose=0)
-
-
-
-
-[docs]class Regressors(MetaAttributes):
-
- def __init__(self, method, cv=None, hyper_method=None, scoring=None):
- self.method = method
- self.cv = cv
- self.hyper_method = hyper_method
- self.scoring = scoring
-
-
- if self.cv == None:
-
- print(".....Using monthly bootstrapper as default splitter....")
- self.cv = MonthlyBooststrapper(n_splits=500, block_size=12)
-
- self.selection = "random"
-
- if hyper_method == None:
- self.hyper_method = "GridSearchCV"
-
-
- if scoring == None:
- self.scoring= ["r2", "neg_root_mean_squared_error"]
-
-
-[docs] def set_model(self):
-
- # regression with variable selction
- if self.method == "LassoCV":
- self.estimator = LassoCV(cv=self.cv, selection=self.selection)
-
- elif self.method == "LassoLarsCV":
- self.estimator = LassoLarsCV(cv=self.cv, normalize=False)
-
- elif self.method == "RidgeCV":
- self.estimator = RidgeCV(cv=self.cv, scoring="r2", alphas=[1e-3, 1e-2, 1e-1, 1, 10])
-
-
- #Bayesian regression algorithms
- elif self.method == "ARD": # Automatic Relevance Determination regression
- self.estimator = ARDRegression(n_iter=500)
-
- elif self.method == "BayesianRidge":
- self.estimator = BayesianRidge(n_iter=500)
-
- #Generalized Linear Models
- elif self.method == "Gamma":
- self.estimator = GammaRegressor()
-
- elif self.method == "Poisson":
- self.estimator = PoissonRegressor()
-
- # Neural Networks models (Perceptron)
- elif self.method == "MLP":
- regressor= MLPRegressor(random_state=42, max_iter=1000, early_stopping=True, batch_size=50,
- n_iter_no_change=20)
- param_grid = {"hidden_layer_sizes": [200,300], "alpha": [0.0001, 1.5, 2,5, 10],
- "learning_rate": ["adaptive"], "solver": ["adam"]}
- self.hyper = HyperparameterOptimize(method=self.hyper_method , param_grid= param_grid, regressor=regressor)
-
- #Support Vector Machines (it very expensive to do hyperparamiter search)
- elif self.method == "SVR":
- regressor = SVR()
- if self.hyper_method == "BayesSearchCV":
- param_grid = {"C":Real(0.1, 10, prior="uniform"),
- "gamma": Real(0.1,10, prior="uniform"),
- "kernel":Categorical(["linear", "poly", "rbf", "sigmoid"]),
- }
- else:
- param_grid = {"C":[0.1, 1, 10], "gamma":[0.1,1, 10],
- "kernel":["linear", "poly", "rbf", "sigmoid"]}
-
- self.hyper = HyperparameterOptimize(method=self.hyper_method, param_grid= param_grid, regressor=regressor)
-
- # Ensemble tree based algorithms
- elif self.method == "RandomForest":
- self.estimator = RandomForestRegressor(n_estimators=100, random_state=42, n_jobs=-1)
-
- elif self.method == "ExtraTree":
- self.estimator = ExtraTreesRegressor(n_estimators=100, random_state=42, n_jobs=-1)
-
- # Bagging
- elif self.method == "Bagging":
- self.estimator = BaggingRegressor(n_estimators=100, random_state=42, n_jobs=-1)
-
- #Boosting
- elif self.method == "AdaBoost":
- self.estimator = AdaBoostRegressor(n_estimators=100, loss= "linear", random_state=42)
-
- elif self.method == "HistGradientBoost":
- self.estimator = HistGradientBoostingRegressor(loss="squared_error", max_iter=200, validation_fraction=0.1,
- random_state=42)
-
- elif self.method == "GradientBoost":
- self.estimator = GradientBoostingRegressor(loss="squared_error", n_estimators=200, validation_fraction=0.1,
- random_state=42, criterion="friedman_mse")
-
- elif self.method == "XGBoost":
- self.estimator = XGBRegressor(n_estimators=100, n_jobs=-1, random_state=42)
-
-
- # gradient descent regression
- elif self.method == "SGD":
- self.estimator = SGDRegressor(loss = "squared_error", max_iter=2000, early_stopping=True,
- random_state=42, validation_fraction=0.1, learning_rate="invscaling",
- )
- else:
- raise ValueError("method name not recognized")
-
-[docs] def fit(self, X,y):
- if self.method == "MLP" or self.method=="SVR":
- self.estimator = self.hyper.fit(X,y)
-
- return self.estimator
- else:
-
- self.estimator = self.estimator.fit(X,y)
- return self.estimator
-
-
-
-
-
-
-[docs] def cross_val_score(self, X, y):
- return cross_val_score(self.estimator, X, y, cv=self.cv, scoring=self.scoring)
-
-[docs] def cross_validate(self, X, y):
- return cross_validate(self.estimator, X, y, scoring=self.scoring,
- n_jobs=2, verbose=0, cv=self.cv)
-
-[docs] def cross_val_predict(self, X, y):
- return cross_val_predict(self.estimator, X, y, n_jobs=2, verbose=0)
-
-
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Mar 16 11:34:25 2022
-
-@author: dboateng
-"""
-
-import matplotlib.pyplot as plt
-import pandas as pd
-import numpy as np
-import seaborn as sns
-import os
-import matplotlib as mpl
-import seaborn as sns
-from matplotlib.dates import YearLocator
-import matplotlib.dates as mdates
-from cycler import cycler
-
-
-try:
- from plot_utils import *
-
-except:
- from .plot_utils import *
-
-[docs]def plot_monthly_mean(means, stds, color, ylabel=None, ax=None,
- fig_path=None, fig_name=None, lolims=False):
-
- if ax is None:
- fig,ax = plt.subplots(1,1, sharex=False, figsize=(20, 15))
-
- plot = means.plot(kind="bar", yerr=stds, rot=0, ax=ax, fontsize=20, capsize=4,
- width=0.8, color=color, edgecolor=black,
- error_kw=dict(ecolor='black',elinewidth=0.5, lolims=lolims))
-
- for ch in plot.get_children():
- if str(ch).startswith("Line2D"):
- ch.set_marker("_")
- ch.set_markersize(10)
- break
-
- if ylabel is not None:
- ax.set_ylabel(ylabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_yticklabels([])
-
- plt.tight_layout()
- if fig_path is not None:
- plt.savefig(os.path.join(fig_path, fig_name), bbox_inches="tight", format= "svg")
-
-[docs]def correlation_heatmap(data, cmap, ax=None, vmax=None, vmin=None, center=0, cbar_ax=None,
- add_cbar=True, title=None, label= "Correlation Coefficinet", fig_path=None, fig_name=None,
- xlabel=None, ylabel=None, fig=None):
-
- if ax is None:
- fig,ax = plt.subplots(1,1, sharex=False, figsize=(15, 13))
-
-
- if add_cbar == True:
- if cbar_ax is None:
- cbar_ax = [0.90, 0.4, 0.02, 0.25]
-
-
- cbar_ax = fig.add_axes(cbar_ax)
- cbar_ax.get_xaxis().set_visible(False)
- cbar_ax.yaxis.set_ticks_position('right')
- cbar_ax.set_yticklabels([])
- cbar_ax.tick_params(size=0)
-
- sns.set(font_scale=1.2)
- if all(parameter is not None for parameter in [vmin, vmax]):
-
- sns.heatmap(ax=ax, data=data, cmap=cmap, vmax=vmax, vmin=vmin, center=center, cbar=add_cbar,
- square=True, cbar_ax = cbar_ax, cbar_kws={"label": label, "shrink":0.5,
- "drawedges": False},
- linewidth=0.5, linecolor="black")
- else:
-
- sns.heatmap(ax=ax, data=data, cmap=cmap, robust=True, cbar=add_cbar,
- square=True, cbar_ax = cbar_ax, cbar_kws={"label": label, "shrink":0.5,
- "drawedges": False},
- linewidth=0.5, linecolor="black")
-
- if xlabel is not None:
- ax.set_xlabel(xlabel, fontsize=18)
- ax.set_ylabel(ylabel, fontsize=18)
-
- plt.tight_layout()
-
- plt.subplots_adjust(left=0.15, right=0.88, top=0.97, bottom=0.05)
-
- if fig_path is not None:
- plt.savefig(os.path.join(fig_path, fig_name), bbox_inches="tight", format= "svg")
-
-
-[docs]def barplot(methods, stationnames, path_to_data, ax=None, xlabel=None, ylabel=None,
- varname="test_r2", varname_std="test_r2_std", filename="validation_score_", legend=True,
- fig_path=None, fig_name=None, show_error=False, width=0.5, rot=0):
-
- if ax is None:
- fig,ax = plt.subplots(1,1, sharex=False, figsize=(18, 15))
-
- df, df_std = barplot_data(methods, stationnames, path_to_data, varname=varname, varname_std=varname_std,
- filename=filename, use_id=True)
-
- colors = [selector_method_colors[m] for m in methods]
- mpl.rcParams["axes.prop_cycle"] = cycler("color", colors)
-
- if show_error == True:
-
- df.plot(kind="bar", yerr=df_std, rot=rot, ax=ax, legend = legend, fontsize=20, capsize=4,
- width=width, edgecolor=black)
- else:
-
- df.plot(kind="bar", rot=rot, ax=ax, legend = legend, fontsize=20, width=width, edgecolor=black)
-
-
-
- if ylabel is not None:
- ax.set_ylabel(ylabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_yticklabels([])
-
- if xlabel is not None:
- ax.set_xlabel(xlabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_xticklabels([])
-
-
-
- if legend ==True:
- ax.legend(loc="upper right", bbox_to_anchor=(1.15, 1), borderaxespad=0., frameon=True, fontsize=20)
- plt.tight_layout()
- plt.subplots_adjust(left=0.05, right=0.95, top=0.97, bottom=0.05)
-
- if fig_path is not None:
- plt.savefig(os.path.join(fig_path, fig_name), bbox_inches="tight", format= "svg")
-
-
-[docs]def boxplot(regressors, stationnames, path_to_data, ax=None, xlabel=None, ylabel=None,
- varname="test_r2", filename="validation_score_",
- fig_path=None, fig_name=None, colors=None, patch_artist=False, rot=45):
-
- if ax is None:
- fig,ax = plt.subplots(1,1, sharex=False, figsize=(20, 15))
-
-
- scores = boxplot_data(regressors, stationnames, path_to_data, filename=filename,
- varname=varname)
-
- color = { "boxes": black,
- "whiskers": black,
- "medians": red,
- "caps": black,
- }
-
-
- boxplot = scores.plot(kind= "box", rot=rot, ax=ax, fontsize=20, color= color, sym="+b", grid=False,
- widths=0.9, notch=False, patch_artist=patch_artist, return_type="dict")
-
-
- if colors is not None:
-
- for patch, color in zip(boxplot["boxes"], colors):
- patch.set_facecolor(color)
-
-
- if ylabel is not None:
- ax.set_ylabel(ylabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_yticklabels([])
-
- if xlabel is not None:
- ax.set_xlabel(xlabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_xticklabels([])
-
- plt.tight_layout()
- plt.subplots_adjust(left=0.05, right=0.95, top=0.97, bottom=0.05)
-
- if fig_path is not None:
- plt.savefig(os.path.join(fig_path, fig_name), bbox_inches="tight", format= "svg")
-
-
-[docs]def heatmaps(data, cmap, label=None, title=None, vmax=None, vmin=None, center=None, ax=None,
- cbar=True, cbar_ax=None, xlabel=None):
-
- if ax is None:
- fig,ax = plt.subplots(1,1, figsize=(20,15))
- plt.subplots_adjust(left=0.02, right=1-0.02, top=0.94, bottom=0.45, hspace=0.25)
-
- if all(parameter is not None for parameter in [vmax, vmin, center]):
-
- if cbar == False:
- sns.heatmap(data=data, ax=ax, cmap=cmap, vmax=vmax, vmin=vmin, center=center,
- square=True, cbar=cbar, linewidth=0.3, linecolor="black")
- else:
-
- if cbar_ax is not None:
- sns.heatmap(data=data, ax=ax, cmap=cmap, vmax=vmax, vmin=vmin, center=center,
- square=True, cbar=cbar, cbar_kws={"label":label,
- "shrink":.80,},
- linewidth=0.3, linecolor="black", cbar_ax=cbar_ax)
- else:
- sns.heatmap(data=data, ax=ax, cmap=cmap, vmax=vmax, vmin=vmin, center=center,
- square=True, cbar=cbar, cbar_kws={"label":label,
- "shrink":.80,},
- linewidth=0.3, linecolor="black")
- else:
- if cbar == False:
- sns.heatmap(data=data, ax=ax, cmap=cmap, square=True, cbar=cbar, linewidth=0.3,
- linecolor="black")
- else:
-
- if cbar_ax is not None:
- sns.heatmap(data=data, ax=ax, cmap=cmap, square=True, cbar=cbar,
- cbar_kws={"label":label,"shrink":.80,},
- linewidth=0.3, linecolor="black", cbar_ax=cbar_ax)
- else:
- sns.heatmap(data=data, ax=ax, cmap=cmap, square=True, cbar=cbar,
- cbar_kws={"label":label,"shrink":.80,},
- linewidth=0.3, linecolor="black")
-
- if title is not None:
- ax.set_title(title, fontsize=20, fontweight="bold", loc="left")
-
- if xlabel is not None:
- ax.set_xlabel(xlabel, fontweight="bold", fontsize=20)
-
- else:
- ax.set_xticklabels([])
-
-
-
-
-
-[docs]def scatterplot(station_num, stationnames, path_to_data, filename, ax=None,
- obs_train_name="obs 1958-2010",
- obs_test_name="obs 2011-2020",
- val_predict_name="ERA5 1958-2010",
- test_predict_name="ERA5 2011-2020",
- method = "Stacking", ylabel=None, xlabel=None,
- fig_path=None, fig_name=None,
- ):
-
- if ax is None:
- fig,ax = plt.subplots(1,1, figsize=(20,15))
- plt.subplots_adjust(left=0.02, right=1-0.02, top=0.94, bottom=0.45, hspace=0.25)
-
- station_info = prediction_example_data(station_num, stationnames, path_to_data, filename,
- obs_test_name=obs_test_name, obs_train_name=obs_train_name,
- val_predict_name=val_predict_name, test_predict_name=test_predict_name,
- method=method)
-
- obs_train = station_info["obs_train"]
- obs_test = station_info["obs_test"]
- ypred_train = station_info["ypred_train"]
- ypred_test = station_info["ypred_test"]
- obs = station_info["obs"]
-
- from scipy import stats
-
- regression_stats = stats.linregress(obs_test, ypred_test)
-
- regression_slope = regression_stats.slope * obs + regression_stats.intercept
-
- r2 = regression_stats.rvalue
-
- ax.scatter(obs_train, ypred_train, alpha=0.3, c=black, s=100, label=val_predict_name)
-
- ax.scatter(obs_test, ypred_test, alpha=0.3, c=red, s=100, label=test_predict_name)
-
- ax.plot(obs, regression_slope, color=red, label="R = {:.2f}".format(r2))
-
- ax.legend(loc= "upper left", fontsize=20)
-
- # Plot design
-
- ax.spines["top"].set_visible(False)
- ax.spines["right"].set_visible(False)
- ax.get_xaxis().tick_bottom()
- ax.get_yaxis().tick_left()
- ax.spines["left"].set_position(("outward", 20))
- ax.spines["bottom"].set_position(("outward", 20))
-
-
-
- if ylabel is not None:
- ax.set_ylabel(ylabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_yticklabels([])
-
- if xlabel is not None:
- ax.set_xlabel(xlabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_xticklabels([])
-
- plt.tight_layout()
- plt.subplots_adjust(left=0.05, right=0.95, top=0.97, bottom=0.05)
-
- if fig_path is not None:
- plt.savefig(os.path.join(fig_path, fig_name), bbox_inches="tight", format= "svg")
-
-
-
-
-[docs]def lineplot(station_num, stationnames, path_to_data, filename, ax=None, fig=None,
- obs_train_name="obs 1958-2010",
- obs_test_name="obs 2011-2020",
- val_predict_name="ERA5 1958-2010",
- test_predict_name="ERA5 2011-2020",
- method = "Stacking", ylabel=None, xlabel=None,
- fig_path=None, fig_name=None,
- ):
-
-
-
- if ax is None:
- fig, ax = plt.subplots(1, 1, figsize= (20, 15), sharex=True)
-
- plt.subplots_adjust(left=0.12, right=1-0.01, top=0.98, bottom=0.06, hspace=0.01)
-
- station_info = prediction_example_data(station_num, stationnames, path_to_data, filename,
- obs_test_name=obs_test_name, obs_train_name=obs_train_name,
- val_predict_name=val_predict_name, test_predict_name=test_predict_name,
- method=method)
-
-
- ypred_train = station_info["ypred_train"].rolling(3, min_periods=1, win_type="hann",
- center=True).mean()
- ypred_test = station_info["ypred_test"].rolling(3, min_periods=1, win_type="hann",
- center=True).mean()
- obs = station_info["obs"].rolling(3, min_periods=1, win_type="hann",
- center=True).mean()
-
-
-
- ax.plot(obs, linestyle="-", color=green, label="Obs")
- ax.plot(ypred_train, linestyle="-.", color= blue, label = val_predict_name)
- ax.plot(ypred_test, linestyle="--", color=red, label=test_predict_name)
-
-
- ax.xaxis.set_major_locator(YearLocator(5))
- ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
- ax.axhline(y=0, linestyle="--", color=grey, linewidth=2)
- ax.legend(bbox_to_anchor=(0.01, 1.02, 1., 0.102), loc=3, ncol=3, borderaxespad=0., frameon = True,
- fontsize=20)
-
-
- if ylabel is not None:
- ax.set_ylabel(ylabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_yticklabels([])
-
- if xlabel is not None:
- ax.set_xlabel(xlabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_xticklabels([])
-
-
- plt.tight_layout()
-
- if fig_path is not None:
- plt.savefig(os.path.join(fig_path, fig_name), bbox_inches="tight", format= "svg")
-
-
-[docs]def plot_time_series(stationnames, path_to_data, filename, id_name,
- daterange, color, label, ymax=None, ymin=None, ax=None, ylabel=None, xlabel=None,
- fig_path=None, fig_name=None, method="Stacking", window=12):
-
- if ax is None:
- fig, ax = plt.subplots(1, 1, figsize= (20, 15), sharex=True)
-
- df = extract_time_series(stationnames, path_to_data, filename, id_name,
- method, daterange,)
-
-
- df = df.rolling(window, min_periods=1, win_type="hann", center=True).mean()
-
- ax.plot(df["mean"], "--", color=color, label=label)
-
-
- #try with max and min to notice the difference with 5 years window
-
- ax.fill_between(df.index, df["mean"] - df["std"], df["mean"] + df["std"], color=color,
- alpha=0.2,)
-
- ax.xaxis.set_major_locator(YearLocator(10))
- ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
- ax.axhline(y=0, linestyle="--", color=grey, linewidth=2)
-
- ax.legend(frameon=True, fontsize=12, loc="lower left")
- if ymax is not None:
-
- ax.set_ylim([ymin, ymax])
-
-
- if ylabel is not None:
- ax.set_ylabel(ylabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_yticklabels([])
-
- # if xlabel is not None:
- # ax.set_xlabel(xlabel, fontweight="bold", fontsize=20)
- # ax.grid(True)
- # else:
- # ax.grid(True)
- # ax.set_xticklabels([])
-
-
- plt.tight_layout()
-
- if fig_path is not None:
- plt.savefig(os.path.join(fig_path, fig_name), bbox_inches="tight", format= "svg")
-
-
-
-[docs]def plot_projection_comparison(stationnames, path_to_data,
- filename, id_name, method, stationloc_dir,
- daterange, datasets, variable, dataset_varname,
- ax=None, xlabel=None, ylabel=None, legend=True,
- figpath=None, figname=None, width=0.5, title=None):
-
- df = extract_comparison_data_means(stationnames, path_to_data, filename, id_name,
- method, stationloc_dir, daterange, datasets,
- variable, dataset_varname)
-
- models_col_names = ["ESD", "MPIESM", "CESM5", "HadGEM2"]
-
-
- if ax is None:
- fig,ax = plt.subplots(1,1, sharex=False, figsize=(18, 15))
-
- colors = [Models_colors[c] for c in models_col_names]
- mpl.rcParams["axes.prop_cycle"] = cycler("color", colors)
-
-
- df.plot(kind="bar", rot=0, ax=ax, legend=legend, fontsize=20, width=width)
-
-
- if ylabel is not None:
- ax.set_ylabel(ylabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_yticklabels([])
-
- if xlabel is not None:
- ax.set_xlabel(xlabel, fontweight="bold", fontsize=20)
- ax.grid(True, linestyle="--", color=gridline_color)
- else:
- ax.grid(True, linestyle="--", color=gridline_color)
- ax.set_xticklabels([])
-
- if title is not None:
- ax.set_title(title, fontsize=20, fontweight="bold", loc="left")
-
-
-
- if legend ==True:
- ax.legend(loc="upper right", bbox_to_anchor=(1.15, 1), borderaxespad=0., frameon=True, fontsize=20)
- plt.tight_layout()
- plt.subplots_adjust(left=0.05, right=0.95, top=0.97, bottom=0.05)
-
-
-
-
-
-
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Apr 11 09:03:49 2022
-
-@author: dboateng
-"""
-
-import matplotlib.pyplot as plt
-import matplotlib as mpl
-import matplotlib.patches as mpatches
-import pandas as pd
-import numpy as np
-from cycler import cycler
-from matplotlib import rc
-
-try:
- from ESD_utils import load_all_stations, load_csv, load_pickle
-except:
- from .ESD_utils import load_all_stations, load_csv, load_pickle
-
-# A4 paper size: 210 mm X 297 mm
-
-cm = 0.3937 # 1 cm in inch for plot size
-pt = 1/72.27 # pt in inch from latex geometry package
-textwidth = 345*pt
-big_width = textwidth + 2*3*cm
-
-
-
-# colors
-gridline_color = "#d5dbdb"
-indianred = "#ec7063"
-seablue = "#5dade2"
-orange = 'orangered'
-lightblue = 'teal'
-brown = 'sienna'
-red = 'red'
-blue = 'blue'
-green = '#55a868'
-purple = '#8172b2'
-lightbrown = '#ccb974'
-pink = 'fuchsia'
-lightgreen = 'lightgreen'
-skyblue = "skyblue"
-tomato = "tomato"
-gold = "gold"
-magenta = "magenta"
-black = "black"
-grey = "grey"
-
-#specific colors
-
-recursive_color = "#D35400"
-treebased_color = "#1abc9c"
-sequential_color = "#1C2833"
-
-esd_color = "#212F3D"
-mpi_color = "#52BE80"
-cesm_color = "#CB4335"
-had_color = "#D4AC0D"
-cordex_color = "#8E44AD"
-
-
-RdBu_r = plt.cm.RdBu_r
-RdBu = plt.cm.RdBu
-BrBG = plt.cm.BrBG
-seismic = plt.cm.seismic
-
-selector_method_colors = {
-
- "Recursive": recursive_color,
- "TreeBased": treebased_color,
- "Sequential": sequential_color,
- }
-
-
-Models_colors = {
- "ESD" : esd_color,
- "MPIESM": mpi_color,
- "CESM5": cesm_color,
- "HadGEM2": had_color,
- "CORDEX": cordex_color}
-
-
-
-
-
-# colors for individual predictors
-
-
-predictor_colors = {
- 'r250':green, 'r500':lightgreen, 'r700':green, 'r850':green, 'r1000':green,
- 'z250':pink, 'z500':green, 'z700':green, 'z850':green, 'z1000':green,
- 't250':tomato, 't500':green, 't700':green, 't850':purple, 't1000':green,
- 'v250':green, 'v500':green, 'v700':green, 'v850':green, 'v1000':green,
- 'u250':green, 'u500':green, 'u700':green, 'u850':green, 'u1000':green,
- 'dt250':green, 'dt500':green, 'dt700':green, 'dt850':lightblue, 'dt1000':green,
- 'NAO':magenta,
- 'EA':brown,
- 'EAWR':brown,
- 'SCAN':brown,
- 't2m':gold,
- 'msl':blue,
- "u10": skyblue
- }
-
-[docs]def apply_style(fontsize=20, style=None, linewidth=2):
- """
-
-
- Parameters
- ----------
- fontsize : TYPE, optional
- DESCRIPTION. The default is 10.
- style : TYPE, optional
- DESCRIPTION. The default is "bmh". ["seaborn", "fivethirtyeight",]
-
- Returns
- -------
- None.
-
- """
- if style is not None:
- plt.style.use(style)
-
- rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
- mpl.rc('text', usetex=True)
- mpl.rc('xtick', labelsize=fontsize)
- mpl.rc('ytick', labelsize=fontsize)
- mpl.rc('legend', fontsize=fontsize)
- mpl.rc('axes', labelsize=fontsize)
- mpl.rc('lines', linewidth=linewidth)
- mpl.rc("font", weight="bold")
-
-
-
-
-[docs]def barplot_data(methods, stationnames, path_to_data, varname="test_r2", varname_std="test_r2_std",
- filename="validation_score_", use_id=False):
-
- df = pd.DataFrame(index=stationnames, columns=methods,)
- df_std = pd.DataFrame(index=stationnames, columns=methods,)
-
- for method in methods:
- scores = load_all_stations(filename + method, path_to_data, stationnames)
- df[method] = scores[varname]
-
- if varname_std is not None:
- df_std[method] = scores[varname_std]
-
- if use_id == True:
-
- df.reset_index(drop=True, inplace=True)
- df.index += 1
- df_std.reset_index(drop=True, inplace=True)
- df_std.index += 1
-
- return df, df_std
-
-
-[docs]def correlation_data(stationnames, path_to_data, filename, predictors, use_id=False):
-
- df = pd.DataFrame(index=stationnames, columns=predictors)
-
-
- for i,idx in enumerate(stationnames):
-
- df.iloc[i] = load_csv(idx, filename, path_to_data)
-
- df = df.astype(float)
-
- if use_id == True:
- df.reset_index(drop=True, inplace=True)
- df.index += 1
-
- return df
-
-[docs]def count_predictors(methods, stationnames, path_to_data, filename, predictors,
- ):
-
- df_count = pd.DataFrame(index=methods, columns=predictors)
-
-
- for method in methods:
-
- df = pd.DataFrame()
- for i, idx in enumerate(stationnames):
-
- selected = load_pickle(idx, filename + method, path_to_data)
-
- df = pd.concat([df, selected], axis=0)
-
- for predictor in predictors:
-
- if predictor in df:
- df_count.loc[method, predictor] = len(df[predictor].dropna())
-
- df_count = df_count.astype(float)
-
- return df_count
-
-[docs]def boxplot_data(regressors, stationnames, path_to_data, filename="validation_score_",
- varname="test_r2"):
-
- df = pd.DataFrame(index=stationnames, columns=regressors)
-
- for regressor in regressors:
- scores = load_all_stations(filename + regressor, path_to_data, stationnames)
-
- df[regressor] = scores[varname]
-
- df.reset_index(drop=True, inplace=True)
- df.index += 1
-
- return df
-
-
-
-[docs]def resample_seasonally(data, daterange):
- df = data[daterange].resample("Q-NOV").mean()
- winter = df[df.index.quarter == 1].mean()
- spring = df[df.index.quarter == 2].mean()
- summer = df[df.index.quarter == 3].mean()
- autumn = df[df.index.quarter == 4].mean()
-
- return winter, spring, summer, autumn
-
-[docs]def resample_monthly(data, daterange):
-
- df = data.resample("MS").mean()
- df = df[daterange]
-
- Jan = df[df.index.month == 1].mean()
- Feb = df[df.index.month == 2].mean()
- Mar = df[df.index.month == 3].mean()
- Apr = df[df.index.month == 4].mean()
- May = df[df.index.month == 5].mean()
- Jun = df[df.index.month == 6].mean()
- Jul = df[df.index.month == 7].mean()
- Aug = df[df.index.month == 8].mean()
- Sep = df[df.index.month == 9].mean()
- Oct = df[df.index.month == 10].mean()
- Nov = df[df.index.month == 11].mean()
- Dec = df[df.index.month == 12].mean()
-
-
- month_means = [Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec]
-
- return month_means
-
-
-[docs]def seasonal_mean(stationnames, path_to_data, filename, daterange, id_name,
- method, use_id=False):
-
- columns = ["DJF", "MAM", "JJA", "SON", "Annum"]
-
- df_stations = pd.DataFrame(index=stationnames, columns=columns)
-
- filename = filename + method
-
- for i,stationname in enumerate(stationnames):
-
- df = load_csv(stationname, filename, path_to_data)
- obs = df[id_name]
- winter, spring, summer, autumn = resample_seasonally(obs, daterange)
-
- obs_mean = obs[daterange].mean()
-
- means = [ winter, spring, summer, autumn, obs_mean]
-
- for j,season in enumerate(columns):
- df_stations.loc[stationname][season] = means[j]
-
-
- if use_id == True:
- df_stations.reset_index(drop=True, inplace=True)
- df_stations.index +=1
-
- df_stations = df_stations.T.astype(float)
-
- return df_stations
-
-[docs]def monthly_mean(stationnames, path_to_data, filename, daterange, id_name,
- method, use_id=False):
-
- import calendar
- month_names = [calendar.month_abbr[im+1] for im in np.arange(12)]
- df_stations = pd.DataFrame(index=stationnames, columns=month_names)
- filename = filename + method
-
- for i,stationname in enumerate(stationnames):
-
- df = load_csv(stationname, filename, path_to_data)
- obs = df[id_name]
- month_means = resample_monthly(obs, daterange)
-
- for j,month in enumerate(month_names):
- df_stations.loc[stationname][month] = month_means[j]
-
- if use_id == True:
- df_stations.reset_index(drop=True, inplace=True)
- df_stations.index +=1
-
-
- df_stations = df_stations.T.astype(float)
-
- return df_stations
-
-
-[docs]def prediction_example_data(station_num, stationnames, path_to_data, filename,
- obs_train_name="obs 1958-2010",
- obs_test_name="obs 2011-2020",
- val_predict_name="ERA5 1958-2010",
- test_predict_name="ERA5 2011-2020",
- method = "Stacking",
- use_cv_all=False):
-
- stationname = stationnames[station_num]
- print("extracting information for the station: ", stationname)
-
- filename = filename + method
- df = load_csv(stationname, filename, path_to_data)
-
- obs_train = df[obs_train_name].dropna()
- obs_test = df[obs_test_name].dropna()
- ypred_validation = df[val_predict_name][~np.isnan(df[obs_train_name])]
- ypred_test = df[test_predict_name][~np.isnan(df[obs_test_name])]
- obs_full = df["obs anomalies"].dropna()
-
-
- # if use_cv_all == True:
-
- # scores_df = load_pickle(stationname, varname="CV_scores_" + method, path=path_to_data)
-
-
- # test_score = scores_df["test_r2"].loc[stationname].max()
-
- # else:
- # test_score_filename = "test_score_" + method
-
- # test_score = load_pickle(stationname, varname=test_score_filename,
- # path=path_to_data)
-
-
- # load score
-
- # validation_score_filename = "validation_score_" + method
-
-
- # validation_score = load_pickle(stationname, varname=validation_score_filename ,
- # path=path_to_data)
-
-
- station_info = {"obs_train": obs_train,
- "obs_test": obs_test,
- "ypred_train": ypred_validation,
- "ypred_test": ypred_test,
- "obs": obs_full}
-
- return station_info
-
-
-
-[docs]def extract_time_series(stationnames, path_to_data, filename, id_name, method,
- daterange):
-
- df = pd.DataFrame(columns=stationnames)
-
- filename = filename + method
-
- for i,stationname in enumerate(stationnames):
-
- stn_data = load_csv(stationname, filename, path_to_data)
-
- pred_data = stn_data[id_name][daterange].dropna()
-
- df[stationname] = pred_data
-
- # calculation of statistics
-
- df["mean"] = df.mean(axis=1)
- df["std"] = df.std(axis=1)
- df["max"] = df.max(axis=1)
- df["min"] = df.min(axis=1)
-
- return df
-
-
-[docs]def extract_comparison_data_means(stationnames, path_to_data,
- filename, id_name, method, stationloc_dir,
- daterange, datasets, variable, dataset_varname,
- ):
-
- models_col_names = ["ESD", "MPIESM", "CESM5", "HadGEM2", "CORDEX"]
-
- df = pd.DataFrame(index=stationnames, columns=models_col_names)
-
- df_stn_loc = pd.read_csv(stationloc_dir, index_col=False,
- usecols=["Latitude", "Longitude"])
-
- filename = filename + method
- for i,stationname in enumerate(stationnames):
- stn_data = load_csv(stationname, filename, path_to_data)
-
- df["ESD"].loc[stationname] = stn_data[id_name][daterange].dropna().mean()
-
- lon = df_stn_loc.iloc[i]["Longitude"]
- lat = df_stn_loc.iloc[i]["Latitude"]
-
- for j in range(len(datasets)):
-
- data = datasets[j].get(dataset_varname, is_Dataset=True)
- if hasattr(data, "rlat"):
- df_proj = data.sel(rlat=lat, rlon=lon, method= "nearest").to_series()
-
- else:
- data = data.sortby("lon")
- df_proj = data.sel(lat=lat, lon=lon, method= "nearest").to_series()
-
- df_proj_mean = df_proj[daterange].mean()
-
- if variable == "Temperature":
- df[models_col_names[j+1]].loc[stationname] = df_proj_mean - 273.15
- elif variable == "Precipitation":
- df[models_col_names[j+1]].loc[stationname] = df_proj_mean *60*60*24*30
-
-
- df.reset_index(drop=True, inplace=True)
- df.index +=1
- df = df.astype(float)
-
- return df
-
-
-
-
-
-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Sun Nov 21 00:55:22 2021
-
-@author: dboateng
-"""
-import sys
-import numpy as np
-import pandas as pd
-import seaborn as sns
-import matplotlib.pyplot as plt
-from copy import copy
-from collections import OrderedDict
-
-
-try:
- from standardizer import MonthlyStandardizer, NoStandardizer
- from feature_selection import RecursiveFeatureElimination, TreeBasedSelection, SequentialFeatureSelection
- from models import Regressors
- from ensemble_models import EnsembleRegressor
- from metrics import Evaluate
- from MLR_model import BootstrappedForwardSelection, MultipleLSRegression
-
-except:
- from .standardizer import MonthlyStandardizer, NoStandardizer
- from .feature_selection import RecursiveFeatureElimination, TreeBasedSelection, SequentialFeatureSelection
- from .models import Regressors
- from .ensemble_models import EnsembleRegressor
- from .metrics import Evaluate
- from .MLR_model import BootstrappedForwardSelection, MultipleLSRegression
-
-
-
-
-
-[docs]class PredictandTimeseries():
-
- def __init__(self, data, transform=None, standardizer=None):
-
- self.data = data
- self.transform = None
- self.standardizer = None
-
- if transform is not None:
- self.set_transform(transform)
- self.set_standardizer(standardizer)
-
-
-[docs] def get(self, daterange=None, anomalies=True):
- if anomalies:
- if daterange is not None:
- return self.data_st.loc[daterange]
- else:
- return copy(self.data_st)
- else:
- if daterange is not None:
- return self.data.loc[daterange]
- else:
- return copy(self.data)
-
-[docs] def set_transform(self, transform):
- if transform == "log":
- self.transform = np.log
- self.backtransform = np.exp
- else:
- raise NotImplementedError("Yet to implement additional transformers to log")
- self.data_st = self.transform(self.data)
-
-[docs] def set_standardizer(self, standardizer):
- if standardizer is not None:
- self.standardizer = standardizer
- else:
- self.standardizer = NoStandardizer()
-
- if self.transform is not None:
- self.data_st = self.transform(self.data)
- self.data_st = self.standardizer.fit_transform(self.data_st)
- else:
- self.data_st = self.standardizer.fit_transform(self.data)
-
-[docs] def set_model(self, method, ensemble_learning=False, estimators=None, cv=10, final_estimator_name=None,
- daterange=None, predictor_dataset=None, fit_predictors=True,
- scoring=["r2", "neg_root_mean_squared_error"], MLR_learning=False, **predictor_kwargs):
-
- self.cv = cv
- self.scoring = scoring
-
- if ensemble_learning == True:
- if method not in ["Stacking", "Voting"]:
- raise ValueError("ensemble method used either stacking or voting")
-
-
- if estimators is None:
- raise ValueError("...estimators list must be provided for ensemble models")
-
- regressors = []
- for i in range(len(estimators)) :
- regressor = Regressors(method=estimators[i], cv=self.cv, scoring=self.scoring)
-
- regressor.set_model()
-
- if estimators[i] in ["MLP", "SVR"]:
-
- X = self._get_predictor_data(daterange=daterange, dataset=predictor_dataset,
- fit_predictors = fit_predictors,
- **predictor_kwargs)
-
- y = self.get(daterange = daterange, anomalies=fit_predictors)
-
- X = X.loc[~np.isnan(y)]
-
- y = y.dropna()
-
- regressor.fit(X.values, y.values)
-
- regressors.append((estimators[i], regressor.estimator))
-
-
-
-
-
-
- self.model = EnsembleRegressor(estimators=regressors, cv=self.cv, method=method,
- final_estimator_name=final_estimator_name,
- scoring=self.scoring)
-
- elif MLR_learning == True:
- if method == "OLSForward":
- self.model = BootstrappedForwardSelection(MultipleLSRegression(), cv=self.cv)
-
- else:
-
- self.model = Regressors(method=method, cv=self.cv, scoring=self.scoring)
- self.model.set_model()
-
-
-[docs] def set_predictors(self, predictors):
- self.predictors = OrderedDict()
- for p in predictors:
- self.predictors[p.name] = p
-
- def _get_predictor_data(self, daterange, dataset, fit_predictors=True, **predictor_kwargs):
- Xs = []
-
- for p in self.predictors:
- Xs.append(self.predictors[p].get(daterange, dataset, fit=fit_predictors))
-
- return pd.concat(Xs, axis=1)
-
-[docs] def predictor_correlation(self, daterange, predictor_dataset, fit_predictors=True, fit_predictand=True,
- method="pearson", **predictor_kwargs):
-
- X = self._get_predictor_data(daterange, predictor_dataset, fit_predictors=fit_predictors,
- **predictor_kwargs)
-
- y = self.get(daterange, anomalies=fit_predictand)
-
- # dropna values
-
- X = X.loc[~np.isnan(y)]
-
- y = y.dropna()
-
-
- corr = X.corrwith(other=y, axis=0, drop=True, method=method)
-
- corr = corr.to_frame()
-
- return corr.T # change the code with scipy stats in order to estimate the significance
-
-
-[docs] def fit_predictor(self, name, daterange, predictor_dataset):
-
- if type(name) == list:
-
- for n in name:
- self.predictors[n].fit(daterange, predictor_dataset)
- else:
- self.predictors[name].fit(daterange, predictor_dataset)
-
-
-
-
-[docs] def fit(self, daterange, predictor_dataset, fit_predictors=True , predictor_selector=True, selector_method="Recursive",
- selector_regressor="Ridge", num_predictors=None, selector_direction=None, cal_relative_importance=False,
- fit_predictand=True, impute=False, impute_method=None, impute_order=None,
- **predictor_kwargs):
-
- # checking attributes required before fitting
-
- if not hasattr(self, "model"):
- raise ValueError("...set model before fitting...")
-
- if not hasattr(self, "predictors"):
- raise ValueError("-----define predictor set first with set_predictors method....")
-
-
- X = self._get_predictor_data(daterange, predictor_dataset, fit_predictors=fit_predictors, **predictor_kwargs)
-
- y = self.get(daterange, anomalies=fit_predictand)
-
- # dropna values
-
-
- if impute == False:
-
- X = X.loc[~np.isnan(y)]
-
- y = y.dropna()
- else:
- if impute_method is None:
- raise ValueError("Enter the imputation method, i.e is either linear or spline ..")
-
-
- else:
- y = y.fillna(y.interpolate(method = impute_method,
- order= impute_order))
-
-
- if predictor_selector ==True:
-
- if selector_method == "Recursive":
- self.selector = RecursiveFeatureElimination(regressor_name=selector_regressor)
-
- elif selector_method == "TreeBased":
- self.selector = TreeBasedSelection(regressor_name=selector_regressor)
-
- elif selector_method == "Sequential":
- if num_predictors == None and selector_direction == None:
- self.selector = SequentialFeatureSelection(regressor_name=selector_regressor)
- else:
- self.selector = SequentialFeatureSelection(regressor_name=selector_regressor,
- n_features=num_predictors, direction=selector_direction)
-
- else:
- raise ValueError("....selector method not recognized .....")
-
-
- X = X.loc[~np.isnan(y)] # because the imputation won't fill series of nan for more years (it just interpolate)
-
- y = y.dropna()
-
-
- self.selector.fit(X, y)
- self.selector.print_selected_features(X)
-
-
-
- X_selected = self.selector.transform(X)
-
- self.model.fit(X_selected, y)
-
- else:
- self.model.fit(X, y)
-
- # explained variance for OLS model
- if hasattr(self.model, "regressor"):
- if hasattr(self.model.regressor, "explained_variances"):
- self.explained_variance_predictors = self.model.regressor.explained_variances
-
-
-
- if cal_relative_importance == True:
-
- if not hasattr(self.model, "coef"):
- raise ValueError("The estimator should have coef_attributes..or must be fitted before this method....")
-
- else:
- coef_ = self.model.coef()
-
- if predictor_selector == True:
-
- score = self.model.score(X_selected,y)
- residual = np.sqrt(1 - score)
- normalized_coef_ = coef_ * np.std(X_selected, axis=0)
- else:
- score = self.model.score(X,y)
- residual = np.sqrt(1 - score)
- normalized_coef_ = coef_ * np.std(X, axis=0)
-
- total = residual + np.sum(np.abs(normalized_coef_))
-
- self.predictor_relative_contribution = normalized_coef_ / total
-
-
-
-
-
-[docs] def predict(self, daterange, predictor_dataset, fit_predictors=True, fit_predictand=True,
- **predictor_kwargs):
-
- X = self._get_predictor_data(daterange, predictor_dataset, fit_predictors, **predictor_kwargs)
-
- if not hasattr(self, "selector"):
-
- yhat = pd.Series(data=self.model.predict(X), index=daterange)
-
- else:
- X_selected = self.selector.transform(X)
-
- yhat = pd.Series(data=self.model.predict(X_selected), index=daterange)
-
- if fit_predictand == False:
- if self.standardizer is not None:
- yhat = self.standardizer.inverse_transform(yhat)
-
- if self.transform is not None:
- yhat = self.backtransform(yhat)
-
- return yhat
-
-
-[docs] def cross_validate_and_predict(self, daterange, predictor_dataset, fit_predictand=True, return_cv_scores=False,
- **predictor_kwargs):
-
- X = self._get_predictor_data(daterange, predictor_dataset, **predictor_kwargs)
-
- y = self.get(daterange, anomalies=fit_predictand)
-
- X = X.loc[~np.isnan(y)]
-
- y = y.dropna()
-
- if hasattr(self, "selector"):
- X_selected = self.selector.transform(X)
-
- val_score = self.model.cross_validate(X_selected, y)
- fit_score = self.model.score(X_selected, y)
- y_pred = self.model.cross_val_predict(X_selected, y)
-
-
- else:
- fit_score = self.model.score(X, y)
- val_score = self.model.cross_validate(X, y)
- y_pred = self.model.cross_val_predict(X, y)
-
-
-
- # write code to work for every scoring variable
- if "r2" in self.scoring:
-
- scores = {"test_r2": np.mean(val_score["test_r2"]),
- "test_r2_std": np.std(val_score["test_r2"]),
- "train_r2": fit_score,
- "test_rmse": -np.mean(val_score["test_neg_root_mean_squared_error"]),
- "test_rmse_std": np.std(val_score["test_neg_root_mean_squared_error"]),
- "test_mae": -np.mean(val_score["test_neg_mean_absolute_error"]),
- "test_mae_std": np.std(val_score["test_neg_mean_absolute_error"]),
- }
-
- else:
-
- scores = {"train_r2": fit_score,
- "test_rmse": -np.mean(val_score["test_neg_root_mean_squared_error"]),
- "test_rmse_std": np.std(val_score["test_neg_root_mean_squared_error"]),
- "test_mae": -np.mean(val_score["test_neg_mean_absolute_error"]),
- "test_mae_std": np.std(val_score["test_neg_mean_absolute_error"]),
- "test_mae_precent": -np.mean(val_score["test_neg_mean_absolute_percentage_error"]),
- "test_mae_precent_std": np.std(val_score["test_neg_mean_absolute_percentage_error"]),
- }
-
-
-
- y_pred = pd.DataFrame({"obs": y,
- "prediction" : y_pred})
-
- if return_cv_scores == True:
- return scores, y_pred, val_score
-
- else:
- return scores, y_pred
-
-
-[docs] def evaluate(self, daterange, predictor_dataset, fit_predictand=True, **predictor_kwargs):
-
- y_true = self.get(daterange, anomalies=fit_predictand)
-
- y_pred = self.predict(daterange, predictor_dataset, anomalies=fit_predictand, **predictor_kwargs)
-
- y_pred = y_pred.loc[~np.isnan(y_true)]
-
- y_true = y_true.dropna()
-
-
- self.evaluate = Evaluate(y_true, y_pred)
-
- rmse = self.evaluate.RMSE()
- nse = self.evaluate.NSE()
- mse = self.evaluate.MSE()
- mae = self.evaluate.MAE()
- exp_var = self.evaluate.explained_variance()
- r2 = self.evaluate.R2_score()
- max_error = self.evaluate.max_error()
- adj_r2 = self.evaluate.adjusted_r2()
-
- scores = {"RMSE": rmse,
- "MSE": mse,
- "NSE": nse,
- "MAE": mae,
- "explained_variance": exp_var,
- "r2": r2,
- "max_error": max_error,
- "adj_r2": adj_r2}
-
- return scores
-
-[docs] def ensemble_transform(self, daterange, predictor_dataset, **predictor_kwargs):
-
- X = self._get_predictor_data(daterange, predictor_dataset, **predictor_kwargs)
-
- if not hasattr(self, "selector"):
-
- y_preds = self.model.transform(X)
-
- else:
- X_selected = self.selector.transform(X)
-
- y_preds = self.model.transform(X_selected)
-
- return y_preds
-
-
-[docs] def relative_predictor_importance(self):
-
- if not hasattr(self,"predictor_relative_contribution"):
-
- raise ValueError("The relative varince must be calculated during fit")
-
- return self.predictor_relative_contribution
-
-
-[docs] def selected_names(self):
-
- if not hasattr(self, "selector"):
- raise ValueError("Predictor selection must be defined when fitting the model")
-
- names = self.selector.select_names
-
- names = names.to_frame()
-
- return names.T
-
-
-
-[docs] def tree_based_feature_importance(self, daterange, predictor_dataset, fit_predictand=True, plot=False, **predictor_kwargs):
-
- if not hasattr(self, "selector"):
- raise ValueError("Predictor selection must be defined when fitting the model")
-
- # if not hasattr(self.selector, "feature_importance"):
- # raise TypeError("the feature selector must be treebased")
-
- X = self._get_predictor_data(daterange, predictor_dataset, **predictor_kwargs)
-
- y = self.get(daterange, anomalies=fit_predictand)
-
- X = X.loc[~np.isnan(y)]
-
- y = y.dropna()
-
-
- return self.selector.feature_importance(X,y, plot=plot)
-
-
-[docs] def tree_based_feature_permutation_importance(self, daterange, predictor_dataset, fit_predictand=True,
- plot=False, **predictor_kwargs):
-
- if not hasattr(self, "selector"):
- raise ValueError("Predictor selection must be defined when fitting the model")
-
- # if not hasattr(self.selector.estimator, "feature_importances_"):
- # raise TypeError("the feature selector must be treebased")
-
- X = self._get_predictor_data(daterange, predictor_dataset, **predictor_kwargs)
-
- y = self.get(daterange, anomalies=fit_predictand)
-
-
- X = X.loc[~np.isnan(y)]
-
- y = y.dropna()
-
- return self.selector.feature_importance(X,y, plot=plot)
-
-
-[docs] def climate_score(self, fit_period, score_period, predictor_dataset, **predictor_kwargs):
- """
- How much better the prediction for the given period is then the
- annual mean.
-
- Parameters
- ----------
- fit_period : pd.DatetimeIndex
- Range of data that should will be used for creating the reference prediction.
- score_period : pd.DatetimeIndex
- Range of data for that the prediction score is evaluated
- predictor_dataset : stat_downscaling_tools.Dataset
- The dataset that should be used to calculate the predictors
- predictor_kwargs : keyword arguments
- These arguments are passed to the predictor's get function
-
- Returns
- -------
- cscore : double
- Climate score (similar to rho squared). 1 for perfect fit, 0 for no
- skill, negative for even worse skill than mean prediction.
- """
- if isinstance(self.standardizer, MonthlyStandardizer):
- # the reference prediction is the seasonal cycle
- y_fit_period = self.get(fit_period, anomalies=False)
- new_standardizer = copy(self.standardizer)
- new_standardizer.fit(y_fit_period)
- zero_prediction = pd.Series(np.zeros(len(y_fit_period)), index=fit_period)
- mean_prediction = new_standardizer.inverse_transform(zero_prediction)
- climate_mse = np.mean((y_fit_period - mean_prediction).dropna()**2)
- y_predict_period = self.get(score_period, anomalies=False)
- # Anomalies need to be predicted and the anomalies converted using
- # the new standardizer, otherwise the prediction has more info than the reference
- yhat_anomalies = self.predict(score_period, predictor_dataset, anomalies=True, **predictor_kwargs)
- yhat = new_standardizer.inverse_transform(yhat_anomalies)
- prediction_mse = np.mean((y_predict_period - yhat).dropna()**2)
- else:
- # don't use this together with a transform
- y_fit = self.get(fit_period, anomalies=True).dropna()
- y_predict = self.get(score_period, anomalies=True)
- yhat = self.predict(score_period, predictor_dataset, anomalies=True, **predictor_kwargs)
- error = (y_predict - yhat).dropna()
- prediction_mse = np.mean(error**2)
- climate_mse = np.mean((np.mean(y_fit) - y_fit)**2)
- if climate_mse == 0:
- return 0.0 # identifier for such a case (avoid dividing by 0)
- else:
- return 1 - prediction_mse/climate_mse
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Jan 25 16:52:13 2022
-
-@author: dboateng
-"""
-
-# importing models
-from sklearn.model_selection import KFold, LeaveOneOut, LeaveOneGroupOut, RepeatedKFold, TimeSeriesSplit
-import numpy as np
-
-[docs]class Splitter():
- # try more on how to use the customized splitter with the model fitting
- def __init__(self, method, shuffle=False, n_splits=5):
- self.method = method
- self.shuffle = shuffle
- self.n_splits = n_splits
- self.random_state = None
-
- if self.method == "Kfold":
- self.estimator = KFold(n_splits=self.n_splits, shuffle=self.shuffle, random_state=self.random_state)
- elif self.method == "LeaveOneOut":
- self.estimator = LeaveOneOut()
- elif self.method == "LeaveOneGroupOut":
- self.estimator = LeaveOneGroupOut()
- elif self.method == "RepeatedKFold":
- self.estimator = RepeatedKFold()
- elif self.method == "TimeSeriesSplit":
- self.estimator = TimeSeriesSplit()
-
- else:
- raise ValueError("Invalid splitter method might have been defined")
-
-
-[docs] def get_n_splits(self, X=None, y=None, groups=None):
- return self.method.get_n_splits(X, y, groups)
-
-
-
-[docs]class MonthlyBooststrapper():
- def __init__(self, n_splits=500, test_size=0.1, block_size=12):
- self.n_splits = n_splits
- self.test_size = test_size
- self.block_size = int(block_size)
-
-[docs] def split(self, X, y, groups=None):
- """
- num_blocks * block_size = test_size*num_samples
- --> n_blocks = test_size/block_size*num_samples
-
- Parameters
- ----------
- X : TYPE
- DESCRIPTION.
- y : TYPE
- DESCRIPTION.
- groups : TYPE, optional
- DESCRIPTION. The default is None.
-
- Returns
- -------
- None.
-
- """
- num_samples = len(y)
- num_blocks = round(self.test_size/self.block_size * num_samples)
-
- for i in range(self.n_splits):
- test_mask = np.zeros(num_samples, dtype=np.bool)
-
- for k in range(num_blocks):
- train_mask = np.zeros(num_samples, dtype=np.bool)
-
- for j in range(num_samples - self.block_size):
- train_mask[j] = not test_mask[j] and not test_mask[j+self.block_size]
-
- train = np.where(train_mask)[0]
- rand = np.random.choice(train)
- test_mask[rand:rand + 12] = True
- train= np.where(~test_mask)
- test = np.where(test_mask)
- yield train, test
-
-[docs]class YearlyBootstrapper:
- """
- Splits data in training and test set by picking complete years. You can use it like this::
-
- X = ...
- y = ...
- yb = YearlyBootstrapper(10)
-
- for i, (train, test) in enumerate(yb.split(X, y)):
- X_train, y_train = X.iloc[train], y.iloc[train]
- X_test, y_test = X.iloc[test], y.iloc[test]
- ...
-
-
- Parameters
- ----------
- n_splits : int (optional, default: 500)
- number of splits
- test_size : float (optional, default: 1/3)
- Ratio of test years.
- min_month_per_year : int (optional, default: 9)
- minimum number of months that must be available in a year to use this
- year in the test set.
- """
-
- def __init__(self, n_splits=500, test_size=1/3, min_month_per_year=9):
- self.n_splits = n_splits
- self.test_size = test_size
- self.min_month_per_year = min_month_per_year
-
-[docs] def split(self, X, y, groups=None):
- """
- Returns ``n_splits`` pairs of indices to training and test set.
-
- Parameters
- ----------
- X : pd.DataFrame
- y : pd.Series
- groups : dummy
-
- X and y should both have the same DatetimeIndex as index.
-
- Returns
- -------
- train : array of ints
- Array of indices of training data
- test : array of ints
- Array of indices of test data
- """
- if np.any(X.index.values != y.index.values):
- raise ValueError("X and y must have the same index")
-
- years = X.index.values.astype('datetime64[Y]').astype(int)
- existing_years, counts = np.unique(years, return_counts=True)
- # we only use years with at least self.min_month_per_year for the test
- # set.
- existing_years = existing_years[counts >= self.min_month_per_year]
- N = len(existing_years)
- size = int(self.test_size*N)
-
- for i in range(self.n_splits):
- test_years = np.random.choice(existing_years, size=size,
- replace=False)
- test_mask = np.zeros(len(years), dtype=np.bool)
- for k in test_years:
- test_mask = test_mask | (years == k)
- train = np.where(~test_mask)
- test = np.where(test_mask)
- yield train, test
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Sun Nov 21 00:55:02 2021
-
-@author: dboateng
-"""
-
-from sklearn.base import TransformerMixin, BaseEstimator
-from sklearn.preprocessing import StandardScaler, RobustScaler, Normalizer
-from sklearn.preprocessing import QuantileTransformer, PowerTransformer
-from sklearn.decomposition import PCA, KernelPCA, IncrementalPCA, SparsePCA
-from copy import deepcopy
-import numpy as np
-import pandas as pd
-import scipy.stats as st
-
-
-[docs]def add_seasonal_cycle(t, anomalies, mean):
-
- """
- Adds a seasonal cycle such that
-
- X = anomalies + mean_seasonal_cycle
-
- Parameters
- ----------
- t : numpy array of ints
- time in number of months
- res : numpy array
- Array of standardized values
- mean : array of shape 12 x #columns(res)
- Mean values for each month and each column in res
-
- Returns
- -------
- X : unstandardized values
- """
- N = len(t)
- M = anomalies.shape[1]
- for i in range(N):
- for j in range(M):
- anomalies[i,j] += mean[t[i] % 12,j]
- return anomalies
-
-
-[docs]def remove_seasonal_cycle(t, X, mean):
-
- """
- Inverse operation to add_seasonal_cycle
- """
- N = len(t)
- M = X.shape[1]
- for i in range(N):
- for j in range(M):
- X[i,j] -= mean[t[i] % 12,j]
- return X
-
-
-[docs]def get_annual_mean_cycle(t, X):
- np.seterr("raise")
-
- N = len(t)
- M = X.shape[1]
- mean = np.zeros((12, M))
- counts = np.zeros((12, M), dtype=np.int)
- for i in range(N):
- for j in range(M):
- if not np.isnan(X[i,j]):
- mean[t[i] % 12,j] += X[i,j]
- counts[t[i] % 12,j] += 1
- for i in range(12):
- for j in range(M):
- mean[i,j] = mean[i,j]/counts[i,j]
- return mean
-
-
-[docs]def get_mean_prediction(t, mean):
- N = len(t)
- mean_prediction = np.zeros(N)
- for i in range(N):
- mean_prediction[i] = mean[t[i] % 12, 0]
- return mean_prediction
-
-[docs]class NoStandardizer(BaseEstimator, TransformerMixin):
- """
- This is just a dummy standardizer that does nothing.
- """
-
-
-
-
-
-
-
-
-[docs]class MonthlyStandardizer(BaseEstimator, TransformerMixin):
- """
- Standardizes monthly data that has a seasonal cycle and possibly a linear trend.
-
- Since the seasonal cycle might affect the trend estimation, the
- seasonal cycle is removed first (by subtracting the mean annual cycle)
- and the trend is estimated by linear regression. Afterwards the data is
- scaled to variance 1.
-
- Parameters
- ----------
- detrending : bool, optional (default: False)
- Whether to remove a linear trend
- """
-
- def __init__(self, detrending=False, scaling=False):
- self.detrending = detrending
- self.scaling = scaling
-
-[docs] def fit(self, X, y=None):
- """
- Fits the standardizer to the provided data, i.e. calculates annual mean
- cycle and trends.
-
- Parameters
- ----------
- X : pd.DataFrame or pd.Series
- DataFrame or Series which holds the data
- y : dummy (optional, default: None)
- Not used
-
- Returns
- -------
- self
- """
- # be careful: X might contain NaNs
-
- t = X.index.values.astype('datetime64[M]').astype(int)
-
-
- # if np.issubdtype(X.index.values.dtype, np.datetime64) == False:
- # t = X.index.values.astype('datetime64[M]').astype(int)
-
- # else:
- # t = X.index.values.astype(int)
-
- values = np.array(X.values) # deep copy
- # X might be a Series or a DataFrame
- if values.ndim == 1:
- values = values[:,np.newaxis]
- N, M = values.shape
-
- # get seasonal mean
- self.mean = get_annual_mean_cycle(t, values)
-
-
- # get slope
- if self.detrending:
- # remove seasonal cycle
- values = remove_seasonal_cycle(t, values, self.mean)
-
- self.slopes = np.zeros(M)
- self.intercepts = np.zeros(M)
- for col in range(M):
- # remove NaNs
- nans = np.isnan(values[:,col])
- x = values[~nans,col]
- time = t[~nans]
- slope, intercept, r, p, stderr = st.linregress(time, x)
- self.slopes[col] = slope
- self.intercepts[col] = intercept
- # remove trend
- values -= self.intercepts + np.outer(t, self.slopes)
-
- if self.scaling:
- self.std = np.nanstd(values, axis=0)
-
- return self
-
-[docs] def transform(self, X, y=None):
- """
- Standardizes the values based on the previously calculated parameters
- Parameters
- ----------
- X : pd.DataFrame or pd.Series
- DataFrame or Series which holds the data
- y : dummy (optional, default: None)
- Not used
-
- Returns
- -------
- X_transformed : pd.DataFrame or pd.Series
- Transformed data
- """
- t = X.index.values.astype('datetime64[M]').astype(int)
- values = np.array(X.values) # deep copy
- # X might be a Series or a DataFrame
- if values.ndim == 1:
- values = values[:,np.newaxis]
-
- values = remove_seasonal_cycle(t, values, self.mean)
- if self.detrending:
- values -= self.intercepts + np.outer(t, self.slopes) # remove trend
-
- if self.scaling:
- values /= self.std
-
- if X.values.ndim == 1:
- return pd.Series(data=values[:,0], index=X.index, name=X.name)
- else:
- return pd.DataFrame(data=values, index=X.index, columns=X.columns)
-
-
-[docs] def inverse_transform(self, X, y=None):
- """
- De-standardizes the values based on the previously calculated parameters
- Parameters
- ----------
- X : pd.DataFrame or pd.Series
- DataFrame or Series which holds the standardized data
- y : dummy (optional, default: None)
- Not used
-
- Returns
- -------
- X_unstandardized : pd.DataFrame or pd.Series
- Unstandardized data
- """
- t = X.index.values.astype('datetime64[M]').astype(int)
- values = np.array(X.values)
- # X might be a Series or a DataFrame
- if values.ndim == 1:
- values = values[:,np.newaxis]
-
- if self.scaling:
- values *= self.std
-
- # add the trend
- if self.detrending:
- values += self.intercepts + np.outer(t, self.slopes)
-
- # add seasonal cycle
- values = add_seasonal_cycle(t, values, self.mean)
-
-
- if X.values.ndim == 1:
- return pd.Series(data=values[:,0], index=X.index, name=X.name)
- else:
- return pd.DataFrame(data=values, index=X.index, columns=X.columns)
-
-
-[docs]class StandardScaling(BaseEstimator, TransformerMixin):
-
- def __init__(self, method=None, with_std=True, with_mean=True, unit_variance=False,
- norm="l2"):
- self.method = method
- self.with_std = with_std
- self.with_mean = with_mean
- self.unit_variance = unit_variance
- self.norm = norm
-
-
- if self.method == None or "standardscaler":
- self.scaler = StandardScaler(with_mean=self.with_mean,
- with_std=self.with_std)
- elif self.method == "robustscaler":
- self.scaler = RobustScaler(with_centering=self.with_mean, with_scaling=self.with_std)
-
- elif self.method == "normalize":
- self.scaler = Normalizer(norm=self.norm)
-
- elif self.method == "powertransformer":
- self.scaler = PowerTransformer(method="yeo-johnson", standardize=True)
-
- elif self.method == "quantiletransformer":
- self.scaler = QuantileTransformer(n_quantiles=1000, output_distribution="unifom")
-
- else:
- raise ValueError("The standardizer do not recognize the defined method")
-
-
-[docs] def fit(self, X):
- X_values = X.values
- if X_values.ndim == 1:
- X_values = X_values[:,np.newaxis]
-
- self.scaler.fit(X_values)
-
- return self
-
- # def fit_transform(self, X, y=None):
-
- # values = self.scaler.fit_transform(X=X, y=y)
-
- # if X.values.ndim == 1:
- # return pd.Series(data=values[:,0], index=X.index, name=X.name)
- # else:
- # return pd.DataFrame(data=values, index=X.index, columns=X.columns)
-
-
-[docs] def inverse_transform(self, X):
-
- X_values = X.values
- if X_values.ndim == 1:
- X_values = X_values[:,np.newaxis]
-
- values = self.scaler.inverse_transform(X_values)
-
- if X.values.ndim == 1:
- return pd.Series(data=values[:,0], index=X.index, name=X.name)
- else:
- return pd.DataFrame(data=values, index=X.index, columns=X.columns)
-
-
-[docs] def transform(self, X):
- X_values = X.values
- if X_values.ndim == 1:
- X_values = X_values[:,np.newaxis]
-
- values = self.scaler.transform(X_values)
-
- if X.values.ndim == 1:
- return pd.Series(data=values[:,0], index=X.index, name=X.name)
- else:
- return pd.DataFrame(data=values, index=X.index, columns=X.columns)
-
-
-
-[docs]class PCAScaling(TransformerMixin, BaseEstimator):
-
- def __init__(self, n_components=None, kernel="linear", method=None):
- self.n_components = n_components
- self.kernel = kernel
- self.method = method
-
-
- if self.method == "PCA" or None:
- self.scaler = PCA(n_components=self.n_components)
-
- elif self.method == "IncrementalPCA":
- self.scaler = IncrementalPCA(n_components=self.n_components)
-
- elif self.method == "KernelPCA":
- self.scaler = KernelPCA(n_components=self.n_components, kernel=self.kernel)
-
- elif self.method == "SparsePCA":
- self.scaler = SparsePCA(n_components=self.n_components)
-
- else:
- raise ValueError("The method passed to the PCAscaling object is not defined")
-
-
-[docs] def fit(self, X):
-
- X_values = X.values
- if X_values.ndim == 1:
- X_values = X_values[:,np.newaxis]
-
- self.scaler.fit(X_values)
-
-
-[docs] def fit_transform(self, X):
-
- X_values = X.values
- if X_values.ndim == 1:
- X_values = X_values[:,np.newaxis]
-
-
- values = self.scaler.fit_transform(X_values)
-
- if X.values.ndim == 1:
- return pd.Series(data=values[:,0], index=X.index)
- else:
- return pd.DataFrame(data=values, index=X.index)
-
-
-[docs] def inverse_transform(self, X):
-
- X_values = X.values
- if X_values.ndim == 1:
- X_values = X_values[:,np.newaxis]
-
- values = self.scaler.inverse_transform(X_values)
-
- if X.values.ndim == 1:
- return pd.Series(data=values[:,0], index=X.index)
- else:
- return pd.DataFrame(data=values, index=X.index)
-
-[docs] def transform(self, X):
- X_values = X.values
-
- if X_values.ndim == 1:
- X_values = X_values[:,np.newaxis]
-
- values = self.scaler.transform(X_values)
-
- if X.values.ndim == 1:
- return pd.Series(data=values[:,0], index=X.index)
- else:
- return pd.DataFrame(data=values, index=X.index)
-
-
-
-
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Mar 14 16:58:59 2022
-
-@author: dboateng
-"""
-
-# import modules
-import os
-import sys
-import xarray as xr
-import numpy as np
-import pandas as pd
-from eofs.xarray import Eof
-from sklearn import decomposition
-import cartopy.crs as ccrs
-import cartopy.feature as cfeature
-from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
-from cartopy.util import add_cyclic_point
-
-
-try:
- from .Predictor_Base import Predictor
- from .ESD_utils import map_to_xarray
-except:
- from Predictor_Base import Predictor
- from ESD_utils import map_to_xarray
-
-
-[docs]def eof_analysis(data, neofs, method="eof_package", apply_equal_wtgs=True,
- pcscaling=1):
-
- if hasattr(data, "longitude"):
- data = data.rename({"longitude":"lon", "latitude":"lat"})
-
- if apply_equal_wtgs == True:
- wtgs = np.sqrt(np.abs(np.cos(data.lat * np.pi / 180)))
- data = data * wtgs
-
- if True in data.isnull():
- data = data.dropna(dim="time")
-
- if method == "eof_package":
-
- print("The EoF Package implementation of EOF analysis is used! for the teleconnections")
-
- solver = Eof(data)
- eofs_cov = solver.eofsAsCovariance(neofs=neofs, pcscaling=pcscaling)
-
- eofs_cov = eofs_cov.sortby(eofs_cov.lon)
-
- pcs = solver.pcs(pcscaling=pcscaling, npcs=neofs)
-
- elif method == "sklearn_package":
-
- print("The sklearn implementation of EOF analysis is used! for the teleconnections")
-
- time_mean = data.mean(dim="time")
- pcs = np.empty(neofs)
-
- PCA = decomposition.PCA(n_components=neofs)
-
- X = np.reshape(data.values, (len(data.time), time_mean.size))
-
- PCA.fit(X)
-
- shape = [neofs] + list(np.shape(time_mean))
-
- np_eofs = np.reshape(PCA.components_, shape)
- pcs = PCA.explained_variance_
-
- ls_eofs= []
- for i in range(neofs):
- ls_eofs.append(map_to_xarray(X=np_eofs[i, ...], datarray=data))
-
- eofs_cov = xr.concat(ls_eofs, pd.Index(range(1, neofs+1), name="eof_number"))
-
- eofs_cov = eofs_cov.sortby(eofs_cov.lon)
-
-
- else:
- raise ValueError("Define the method or defaut of EOF package is used")
-
-
-
- if apply_equal_wtgs == True:
- return eofs_cov, pcs, wtgs
-
- else:
- return eofs_cov, pcs
-
-
-[docs]def extract_region(data, datarange, varname, minlat, maxlat, minlon, maxlon):
-
- data = data.get(varname).sel(time=datarange)
-
- if hasattr(data, "longitude"):
- data = data.rename({"longitude":"lon", "latitude":"lat"})
-
- #data = data.assign_coords({"lon": (((data.lon + 180) % 360) - 180)})
-
- data = data.where((data.lat >=minlat) & (data.lat <=maxlat), drop=True)
- data = data.where((data.lon >=minlon) & (data.lon <= maxlon), drop=True)
-
- return data
-
-
-
-
-[docs]class NAO(Predictor):
- def __init__(self, **kwargs):
- super().__init__(name="NAO", longname= "North Atlantic Oscilation", **kwargs)
-
- def _generate(self, datarange, dataset, fit, patterns_from, params_from):
-
- if hasattr(dataset, "expver"):
- data = dataset.drop("expver")
-
- params = self.params[params_from]
-
- data = extract_region(dataset, datarange, varname="msl", minlat=20, maxlat=80,
- minlon=-80, maxlon=60)
-
- # removing monthly cycle
- group = data.groupby("time.month")
-
- if fit:
- params["monthly_means"] = group.mean(dim="time")
-
- anomalies = group.apply(
- lambda x: x - params["monthly_means"].sel(month=_get_month(x[0].time.values))
- )
- params["monthly_means"] = group.mean(dim="time")
-
-
- #anomalies = group - params["monthly_means"]
- anomalies = anomalies.drop("month")
-
- if fit:
- params["std_field"] = anomalies.std(dim="time")
-
- anomalies /= params["std_field"]
-
- area_weighted = anomalies * np.sqrt(np.abs(np.cos(anomalies.lat*np.pi/180)))
-
-
- method_name = "eof_package"
-
- if fit:
-
- self.patterns[dataset.name] = {}
- eofs, pcs, wtgs = eof_analysis(data=anomalies, neofs=1, method=method_name, apply_equal_wtgs=True,
- pcscaling=0)
-
-
- if hasattr(eofs, "eof_number"):
- self.patterns[dataset.name]["eof"] = eofs.sel(eof_number=1)
-
- else:
- self.patterns[dataset.name]["eof"] = eofs.sel(mode=0)
-
- if method_name == "sklearn_package":
- nao = (self.patterns[patterns_from]["eof"] * area_weighted).sum(dim=("lat", "lon"))
-
- else:
-
- nao = pcs.sel(mode=0)
-
- nao.name = "NAO"
-
-
- if fit:
- self.patterns[dataset.name]["std"] = nao.std()
-
- nao /= self.patterns[patterns_from]["std"]
-
- nao_series = nao.to_series()
-
-
- return nao_series
-
-
-
-
-
-[docs]class EA(Predictor):
- def __init__(self, **kwargs):
- super().__init__(name="EA", longname= "East Atlantic Oscilation", **kwargs)
-
- def _generate(self, datarange, dataset, fit, patterns_from, params_from):
-
- if hasattr(dataset, "expver"):
- data = dataset.drop("expver")
-
- params = self.params[params_from]
-
- data = extract_region(dataset, datarange, varname="msl", minlat=20, maxlat=80,
- minlon=-80, maxlon=60)
-
- # removing monthly cycle
- group = data.groupby("time.month")
-
- if fit:
- params["monthly_means"] = group.mean(dim="time")
-
- anomalies = group.apply(
- lambda x: x - params["monthly_means"].sel(month=_get_month(x[0].time.values))
- )
- params["monthly_means"] = group.mean(dim="time")
-
-
- #anomalies = group - params["monthly_means"]
- anomalies = anomalies.drop("month")
-
- if fit:
- params["std_field"] = anomalies.std(dim="time")
-
- anomalies /= params["std_field"]
-
- area_weighted = anomalies * np.sqrt(np.abs(np.cos(anomalies.lat*np.pi/180)))
-
- method_name = "eof_package"
-
- if fit:
-
- self.patterns[dataset.name] = {}
- eofs, pcs, wtgs = eof_analysis(data=anomalies, neofs=2, method=method_name, apply_equal_wtgs=True,
- pcscaling=0)
-
-
- if hasattr(eofs, "eof_number"):
- self.patterns[dataset.name]["eof"] = eofs.sel(eof_number=2)
-
- else:
- self.patterns[dataset.name]["eof"] = eofs.sel(mode=1)
-
- if method_name == "sklearn_package":
- ea = (self.patterns[patterns_from]["eof"] * area_weighted).sum(dim=("lat", "lon"))
-
- else:
-
- ea = pcs.sel(mode=1)
-
- ea.name = "EA"
-
-
- if fit:
- self.patterns[dataset.name]["std"] = ea.std()
-
- ea /= self.patterns[patterns_from]["std"]
-
- ea_series = ea.to_series()
-
-
- return ea_series
-
-
-
-
-[docs]class SCAN(Predictor):
- def __init__(self, **kwargs):
- super().__init__(name="SCAN", longname= "Scandinavian Oscilation", **kwargs)
- def _generate(self, datarange, dataset, fit, patterns_from, params_from):
-
-
- if hasattr(dataset, "expver"):
- data = dataset.drop("expver")
-
- params = self.params[params_from]
-
- data = extract_region(dataset, datarange, varname="msl", minlat=20, maxlat=80,
- minlon=-80, maxlon=60)
-
- # removing monthly cycle
- group = data.groupby("time.month")
-
- if fit:
- params["monthly_means"] = group.mean(dim="time")
-
- anomalies = group.apply(
- lambda x: x - params["monthly_means"].sel(month=_get_month(x[0].time.values))
- )
- params["monthly_means"] = group.mean(dim="time")
-
-
- #anomalies = group - params["monthly_means"]
- anomalies = anomalies.drop("month")
-
- if fit:
- params["std_field"] = anomalies.std(dim="time")
-
- anomalies /= params["std_field"]
-
- area_weighted = anomalies * np.sqrt(np.abs(np.cos(anomalies.lat*np.pi/180)))
-
- method_name = "sklearn_package"
-
- if fit:
-
- self.patterns[dataset.name] = {}
- eofs, pcs, wtgs = eof_analysis(data=anomalies, neofs=3, method=method_name, apply_equal_wtgs=True,
- pcscaling=0)
-
-
- if hasattr(eofs, "eof_number"):
- self.patterns[dataset.name]["eof"] = eofs.sel(eof_number=3)
-
- else:
- self.patterns[dataset.name]["eof"] = eofs.sel(mode=2)
-
- if method_name == "sklearn_package":
- scan = (self.patterns[patterns_from]["eof"] * area_weighted).sum(dim=("lat", "lon"))
-
- else:
-
- scan = pcs.sel(mode=2)
-
- scan.name = "SCAN"
-
-
- if fit:
- self.patterns[dataset.name]["std"] = scan.std()
-
- scan /= self.patterns[patterns_from]["std"]
-
- scan_series = scan.to_series()
-
-
- return scan_series
-
-
-
-
-[docs]class EAWR(Predictor):
- def __init__(self, **kwargs):
- super().__init__(name="EAWR", longname= "East Atlantic_West Russian Oscilation", **kwargs)
- def _generate(self, datarange, dataset, fit, patterns_from, params_from):
-
- if hasattr(dataset, "expver"):
- data = dataset.drop("expver")
-
- params = self.params[params_from]
-
- data = extract_region(dataset, datarange, varname="msl", minlat=20, maxlat=80,
- minlon=-80, maxlon=60)
-
- # removing monthly cycle
- group = data.groupby("time.month")
-
-
- if fit:
- params["monthly_means"] = group.mean(dim="time")
-
- anomalies = group.apply(
- lambda x: x - params["monthly_means"].sel(month=_get_month(x[0].time.values))
- )
- params["monthly_means"] = group.mean(dim="time")
-
-
- #anomalies = group - params["monthly_means"]
- anomalies = anomalies.drop("month")
-
- if fit:
- params["std_field"] = anomalies.std(dim="time")
-
- anomalies /= params["std_field"]
-
- area_weighted = anomalies * np.sqrt(np.abs(np.cos(anomalies.lat*np.pi/180)))
-
- method_name = "sklearn_package"
-
- if fit:
-
- self.patterns[dataset.name] = {}
- eofs, pcs, wtgs = eof_analysis(data=anomalies, neofs=4, method=method_name, apply_equal_wtgs=True,
- pcscaling=0)
-
-
- if hasattr(eofs, "eof_number"):
- self.patterns[dataset.name]["eof"] = eofs.sel(eof_number=4)
-
- else:
- self.patterns[dataset.name]["eof"] = eofs.sel(mode=3)
-
- if method_name == "sklearn_package":
- ea_wr = (self.patterns[patterns_from]["eof"] * area_weighted).sum(dim=("lat", "lon"))
-
- else:
-
- ea_wr = pcs.sel(mode=3)
-
- ea_wr.name = "EAWR"
-
-
- if fit:
- self.patterns[dataset.name]["std"] = ea_wr.std()
-
- ea_wr /= self.patterns[patterns_from]["std"]
-
- ea_wr_series = ea_wr.to_series()
-
-
- return ea_wr_series
-
-
-
-
-[docs]def _get_month(npdatetime64):
-
- """
- Returns the month for a given npdatetime64 object, 1 for January, 2 for
- February, ...
- """
- month = npdatetime64.astype('datetime64[M]').astype(int) % 12 + 1
-
- return month
-' + - '' + - Documentation.gettext("Hide Search Matches") + - "
" - ) - ); - }, - - /** - * helper function to hide the search marks again - */ - hideSearchWords: () => { - document - .querySelectorAll("#searchbox .highlight-link") - .forEach((el) => el.remove()); - document - .querySelectorAll("span.highlighted") - .forEach((el) => el.classList.remove("highlighted")); - const url = new URL(window.location); - url.searchParams.delete("highlight"); - window.history.replaceState({}, "", url); - }, - - /** - * helper function to focus on search bar - */ - focusSearchBar: () => { - document.querySelectorAll("input[name=q]")[0]?.focus(); - }, - - /** - * Initialise the domain index toggle buttons - */ - initDomainIndexTable: () => { - const toggler = (el) => { - const idNumber = el.id.substr(7); - const toggledRows = document.querySelectorAll(`tr.cg-${idNumber}`); - if (el.src.substr(-9) === "minus.png") { - el.src = `${el.src.substr(0, el.src.length - 9)}plus.png`; - toggledRows.forEach((el) => (el.style.display = "none")); - } else { - el.src = `${el.src.substr(0, el.src.length - 8)}minus.png`; - toggledRows.forEach((el) => (el.style.display = "")); - } - }; - - const togglerElements = document.querySelectorAll("img.toggler"); - togglerElements.forEach((el) => - el.addEventListener("click", (event) => toggler(event.currentTarget)) - ); - togglerElements.forEach((el) => (el.style.display = "")); - if (DOCUMENTATION_OPTIONS.COLLAPSE_INDEX) togglerElements.forEach(toggler); - }, - - initOnKeyListeners: () => { - // only install a listener if it is really needed - if ( - !DOCUMENTATION_OPTIONS.NAVIGATION_WITH_KEYS && - !DOCUMENTATION_OPTIONS.ENABLE_SEARCH_SHORTCUTS - ) - return; - - const blacklistedElements = new Set([ - "TEXTAREA", - "INPUT", - "SELECT", - "BUTTON", - ]); - document.addEventListener("keydown", (event) => { - if (blacklistedElements.has(document.activeElement.tagName)) return; // bail for input elements - if (event.altKey || event.ctrlKey || event.metaKey) return; // bail with special keys - - if (!event.shiftKey) { - switch (event.key) { - case "ArrowLeft": - if (!DOCUMENTATION_OPTIONS.NAVIGATION_WITH_KEYS) break; - - const prevLink = document.querySelector('link[rel="prev"]'); - if (prevLink && prevLink.href) { - window.location.href = prevLink.href; - event.preventDefault(); - } - break; - case "ArrowRight": - if (!DOCUMENTATION_OPTIONS.NAVIGATION_WITH_KEYS) break; - - const nextLink = document.querySelector('link[rel="next"]'); - if (nextLink && nextLink.href) { - window.location.href = nextLink.href; - event.preventDefault(); - } - break; - case "Escape": - if (!DOCUMENTATION_OPTIONS.ENABLE_SEARCH_SHORTCUTS) break; - Documentation.hideSearchWords(); - event.preventDefault(); - } - } - - // some keyboard layouts may need Shift to get / - switch (event.key) { - case "/": - if (!DOCUMENTATION_OPTIONS.ENABLE_SEARCH_SHORTCUTS) break; - Documentation.focusSearchBar(); - event.preventDefault(); - } - }); - }, -}; - -// quick alias for translations -const _ = Documentation.gettext; - -_ready(Documentation.init); diff --git a/_static/documentation_options.js b/_static/documentation_options.js deleted file mode 100755 index 077fd7d..0000000 --- a/_static/documentation_options.js +++ /dev/null @@ -1,14 +0,0 @@ -var DOCUMENTATION_OPTIONS = { - URL_ROOT: document.getElementById("documentation_options").getAttribute('data-url_root'), - VERSION: '0.0.1', - LANGUAGE: 'en', - COLLAPSE_INDEX: false, - BUILDER: 'html', - FILE_SUFFIX: '.html', - LINK_SUFFIX: '.html', - HAS_SOURCE: true, - SOURCELINK_SUFFIX: '.txt', - NAVIGATION_WITH_KEYS: false, - SHOW_SEARCH_SUMMARY: true, - ENABLE_SEARCH_SHORTCUTS: false, -}; \ No newline at end of file diff --git a/_static/file.png b/_static/file.png deleted file mode 100755 index a858a41..0000000 Binary files a/_static/file.png and /dev/null differ diff --git a/_static/jquery-3.5.1.js b/_static/jquery-3.5.1.js deleted file mode 100755 index 5093733..0000000 --- a/_static/jquery-3.5.1.js +++ /dev/null @@ -1,10872 +0,0 @@ -/*! - * jQuery JavaScript Library v3.5.1 - * https://jquery.com/ - * - * Includes Sizzle.js - * https://sizzlejs.com/ - * - * Copyright JS Foundation and other contributors - * Released under the MIT license - * https://jquery.org/license - * - * Date: 2020-05-04T22:49Z - */ -( function( global, factory ) { - - "use strict"; - - if ( typeof module === "object" && typeof module.exports === "object" ) { - - // For CommonJS and CommonJS-like environments where a proper `window` - // is present, execute the factory and get jQuery. - // For environments that do not have a `window` with a `document` - // (such as Node.js), expose a factory as module.exports. - // This accentuates the need for the creation of a real `window`. - // e.g. var jQuery = require("jquery")(window); - // See ticket #14549 for more info. - module.exports = global.document ? - factory( global, true ) : - function( w ) { - if ( !w.document ) { - throw new Error( "jQuery requires a window with a document" ); - } - return factory( w ); - }; - } else { - factory( global ); - } - -// Pass this if window is not defined yet -} )( typeof window !== "undefined" ? window : this, function( window, noGlobal ) { - -// Edge <= 12 - 13+, Firefox <=18 - 45+, IE 10 - 11, Safari 5.1 - 9+, iOS 6 - 9.1 -// throw exceptions when non-strict code (e.g., ASP.NET 4.5) accesses strict mode -// arguments.callee.caller (trac-13335). But as of jQuery 3.0 (2016), strict mode should be common -// enough that all such attempts are guarded in a try block. -"use strict"; - -var arr = []; - -var getProto = Object.getPrototypeOf; - -var slice = arr.slice; - -var flat = arr.flat ? function( array ) { - return arr.flat.call( array ); -} : function( array ) { - return arr.concat.apply( [], array ); -}; - - -var push = arr.push; - -var indexOf = arr.indexOf; - -var class2type = {}; - -var toString = class2type.toString; - -var hasOwn = class2type.hasOwnProperty; - -var fnToString = hasOwn.toString; - -var ObjectFunctionString = fnToString.call( Object ); - -var support = {}; - -var isFunction = function isFunction( obj ) { - - // Support: Chrome <=57, Firefox <=52 - // In some browsers, typeof returns "function" for HTML