diff --git a/docs/roadmap.rst b/docs/roadmap.rst index 3ccd8f6..1fdcdbc 100644 --- a/docs/roadmap.rst +++ b/docs/roadmap.rst @@ -88,8 +88,45 @@ Spatial models ~~~~~~~~~~~~~~~~ This category of methods includes models that have spatial dependence. -This far only the BCSD method (see Thrasher et al., 2012) for spatial disaggregation has -been implemented (for temperature and precipitation). +Thus far only the BCSD method (see Thrasher et al., 2012) for spatial disaggregation has +been implemented (for temperature and precipitation). Because the method involves regridding, +which depends on the xESMF package, a recipe is provided herein versus a full implementation +in scikit-downscale. + +.. code-block:: Python + + from skdownscale.spatial_models import SpatialDisaggregator, apply_weights + import xesmf as xe + import xarray as xr + + # da_climo_fine: xarray.DataArray (daily multi-decade climatology) + # da_obs : xarray.DataArray (daily obs data) + # da_model_train: xarray.DataArray (daily training data) + # ds_sd: xarray.Dataset (daily spatially disaggregated data) + + # regrid climo to model resolution + obs_to_mod_weights = "filename" + regridder_obs_to_mod = xe.Regridder(da_obs.isel(time=0), da_model_train.isel(time=0), + 'bilinear', filename=obs_to_mod_weights, reuse_weights=True) + climo_regrid = xr.map_blocks(apply_weights, regridder_obs_to_mod, + args=[da_climo_fine]) + climo_coarse = climo_regrid.compute() + + # fit the scaling factor model + sfc = SpatialDisaggregator.fit(da_model_train, climo_coarse, var_name='tasmax') + + # regrid scale factor + mod_to_obs_weights = "filename" + regridder_mod_to_obs = xe.Regridder(da_model_train.isel(time=0), + da_obs.isel(time=0), 'bilinear', + filename=mod_to_obs_weights, reuse_weights=True) + sff_regrid = xr.map_blocks(apply_weights, regridder_mod_to_obs, args=[sfc]) + sff = sff_regrid.compute() + + # predict using scale factor + ds_varname = 'scale_factor_fine' + sff_ds = sff.to_dataset(name=ds_varname) + ds_sd = SpatialDisaggregator.predict(sff_ds, da_climo_fine, var_name=ds_varname) Global models ~~~~~~~~~~~~~ @@ -119,7 +156,7 @@ Dependencies ------------ - Core: Xarray, Pandas, Dask, Scikit-learn, Numpy, Scipy -- Optional: Statsmodels, Keras, PyTorch, Tensorflow, etc. +- Optional: Statsmodels, Keras, PyTorch, Tensorflow, xESMF, etc. Related projects ---------------- diff --git a/skdownscale/spatial_models/__init__.py b/skdownscale/spatial_models/__init__.py new file mode 100644 index 0000000..3e2aa05 --- /dev/null +++ b/skdownscale/spatial_models/__init__.py @@ -0,0 +1,2 @@ +from .regridding import apply_weights +from .sd import SpatialDisaggregator diff --git a/skdownscale/spatial_models/regridding.py b/skdownscale/spatial_models/regridding.py new file mode 100644 index 0000000..0625889 --- /dev/null +++ b/skdownscale/spatial_models/regridding.py @@ -0,0 +1,7 @@ +import xesmf as xe + +def apply_weights(regridder, input_data): + regridder._grid_in = None + regridder._grid_out = None + result = regridder(input_data) + return result diff --git a/skdownscale/spatial_models/sd.py b/skdownscale/spatial_models/sd.py new file mode 100644 index 0000000..f947413 --- /dev/null +++ b/skdownscale/spatial_models/sd.py @@ -0,0 +1,124 @@ +import numpy as np +import pandas as pd +import xarray as xr + +class SpatialDisaggregator: + """ + Spatial disaggregation model class + Apply spatial disaggregation algorithm to an xarray Dataset with fit + and predict methods using NASA-NEX method for spatial disaggregation + (see Thrasher et al, 2012). + + Parameters + ---------- + var : str + specifies the variable being downscaled. Default is + temperature and other option is precipitation. + """ + + def __init__(self, var='temperature'): + self._var = var + + if var == "temperature": + pass + elif var == "precipitation": + pass + else: + raise NotImplementedError("functionality for spatial disaggregation" + " of %s has not yet been added" %var) + + def fit(self, ds_bc, climo_coarse, var_name, + lat_name='lat', lon_name='lon'): + """ + Fit the scaling factor used in spatial disaggregation + + Parameters + ----------- + ds_bc : xarray.Dataset + Daily bias corrected data at the model resolution + climo_coarse : xarray.DataArray or xarray.Dataset + Observed climatology that has been regridded (coarsened) + to the model resolution + var_name : str + Specifies the data variable name within ds_bc of the + daily bias corrected data + lat_name : str + Name of the latitude dimension of ds_bc and climo_coarse, + default is 'lat'. + lon_name : str + Name of the longitude dimension of ds_bc and climo_coarse, + default is 'lon'. + """ + + # check that climo has been regridded to model res + if not np.array_equal(ds_bc[lat_name], climo_coarse[lat_name]): + raise ValueError("climo latitude dimension does not match model res") + if not np.array_equal(ds_bc[lon_name], climo_coarse[lon_name]): + raise ValueError("climo longitude dimension does not match model res") + + scf = self._calculate_scaling_factor(ds_bc, climo_coarse, var_name, self._var) + + return scf + + def predict(self, scf, climo_fine, var_name, + lat_name='lat', lon_name='lon'): + """ + Predict (apply) the scaling factor to the observed climatology. + + Parameters + ----------- + scf : xarray.Dataset + Scale factor that has been regridded to the resolution of + the observed data. + climo_fine : xarray.DataArray or xarray.Dataset + Observed climatology at its native resolution + var_name : str + Specifies the data variable name within ds_bc of the + daily bias corrected data + lat_name : str + Name of the latitude dimension of ds_bc and climo_coarse, + default is 'lat'. + lon_name : str + Name of the longitude dimension of ds_bc and climo_coarse, + default is 'lon' + """ + + # check that scale factor has been regridded to obs res + if not np.array_equal(scf[lat_name], climo_fine[lat_name]): + raise ValueError("scale factor latitude dimension does not match obs res") + if not np.array_equal(scf[lon_name], climo_fine[lon_name]): + raise ValueError("scale factor longitude dimension does not match obs res") + + downscaled = self._apply_scaling_factor(scf, climo_fine, var_name, self._var) + + return downscaled + + def _calculate_scaling_factor(self, ds_bc, climo, var_name, var): + ''' + compute scaling factor + ''' + # Necessary workaround to xarray's check with zero dimensions + # https://github.com/pydata/xarray/issues/3575 + da = ds_bc[var_name] + if sum(da.shape) == 0: + return da + groupby_type = ds_bc.time.dt.dayofyear + gb = da.groupby(groupby_type) + + if var == 'temperature': + return gb - climo + elif var == 'precipitation': + return gb / climo + + def _apply_scaling_factor(self, scf, climo, var_name, var): + ''' + apply scaling factor + ''' + groupby_type = scf.time.dt.dayofyear + da = scf[var_name] + sff_daily = da.groupby(groupby_type) + + if var == 'temperature': + return sff_daily + climo + elif var == 'precipitation': + return sff_daily * climo