Repository: carbonplan/cmip6-downscaling (language: python)

Issue:
Refactor/fullrun



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml
index 7f681820..17373e12 100644
--- a/.github/workflows/main.yaml
+++ b/.github/workflows/main.yaml
@@ -12,7 +12,7 @@ jobs:
     runs-on: ubuntu-latest
     timeout-minutes: 20
     container:
-      image: carbonplan/cmip6-downscaling-prefect:2022.02.05
+      image: carbonplan/cmip6-downscaling-prefect:latest
       options: --user root
     steps:
       - uses: actions/checkout@v2
diff --git a/.isort.cfg b/.isort.cfg
new file mode 100644
index 00000000..5d7bf33d
--- /dev/null
+++ b/.isort.cfg
@@ -0,0 +1,2 @@
+[tool.isort]
+profile = "black"
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 62080811..26da970d 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -16,7 +16,7 @@ repos:
       args: ["--autofix", "--indent=2", "--no-sort-keys"]
 
 -   repo: https://github.com/psf/black
-    rev: 21.12b0
+    rev: 22.3.0
     hooks:
     - id: black
       args: ["--line-length", "100", "--skip-string-normalization"]
@@ -34,6 +34,7 @@ repos:
     rev: 5.10.1
     hooks:
     - id: isort
+      args: ["--profile", "black"]
 
 -   repo: https://github.com/deathbeds/prenotebook
     rev: f5bdb72a400f1a56fe88109936c83aa12cc349fa
diff --git a/cmip6_downscaling/__init__.py b/cmip6_downscaling/__init__.py
index 0e70844a..a74613de 100644
--- a/cmip6_downscaling/__init__.py
+++ b/cmip6_downscaling/__init__.py
@@ -1,7 +1,14 @@
 from dask.utils import SerializableLock
 from donfig import Config
+from pkg_resources import DistributionNotFound, get_distribution
 
-from cmip6_downscaling.config import _defaults
+from .config import _defaults
+
+try:
+    version = get_distribution(__name__).version
+except DistributionNotFound:  # pragma: no cover
+    version = '0.0.0'  # pragma: no cover
+__version__ = version
 
 config = Config("cmip6_downscaling", defaults=[_defaults])
 config.config_lock = SerializableLock()
diff --git a/cmip6_downscaling/analysis/analyses_template.ipynb b/cmip6_downscaling/analysis/analyses_template.ipynb
index 0eac9338..9ba66628 100644
--- a/cmip6_downscaling/analysis/analyses_template.ipynb
+++ b/cmip6_downscaling/analysis/analyses_template.ipynb
@@ -33,13 +33,15 @@
     "import pytest\n",
     "from cmip6_downscaling.analysis import analysis, metrics\n",
     "from cmip6_downscaling.analysis.analysis import qaqc_checks, grab_top_city_data, load_top_cities, get_seasonal, change_ds\n",
-    "from cmip6_downscaling.analysis.plot import plot_cdfs, plot_values_and_difference, plot_seasonal, plot_each_step_bcsd\n",
+    "from cmip6_downscaling.analysis.plot import plot_cdfs, plot_values_and_difference, plot_seasonal#, plot_each_step_bcsd\n",
     "from cmip6_downscaling.analysis.qaqc import make_qaqc_ds\n",
-    "from cmip6_downscaling.workflows.utils import BBox\n",
+    "from cmip6_downscaling.methods.common.containers import BBox\n",
     "import matplotlib.pyplot as plt\n",
     "import cartopy.crs as ccrs\n",
-    "from cmip6_downscaling.workflows.paths import make_bcsd_output_path\n",
     "from carbonplan import styles\n",
+    "from cmip6_downscaling import config\n",
+    "from upath import UPath\n",
+    "\n",
     "\n",
     "styles.mpl.set_theme(style='carbonplan_light')\n"
    ]
@@ -78,22 +80,21 @@
    "outputs": [],
    "source": [
     "# Pass as a parameter to be used by papermill: the runid which will be used to create the template for the raw and summary files\n",
-    "gcm_identifier = \"MIROC6/ssp370/tasmax/40.0_50.0_80.0_90.0/1981_2010/1981_2099/\"  # \"MIROC6-ssp370-1990-1990-2090-2090-55-60-80-85-tasmax\"\n",
-    "obs_identifier = \"ERA5/tasmax/40.0_50.0_80.0_90.0/1981_2010/\"\n",
-    "gcm_grid_spec = \"128x256_gridsize_14_14_llcorner_-88_-180\"\n",
-    "result_dir = \"az://flow-outputs/testing_results\"\n",
-    "intermediate_dir = \"az://flow-outputs/testing_intermediates\"\n",
-    "var = \"tasmax\"\n",
-    "gcm = \"MIROC6\"\n",
+    "method = \"bcsd\"\n",
+    "obs = \"ERA5\"\n",
+    "model = \"MIROC6\"\n",
+    "member = \"r1i1p1f1\"\n",
+    "grid_label = \"gn\"\n",
+    "table_id = \"day\"\n",
     "scenario = \"ssp370\"\n",
-    "train_period_start = \"1981\"\n",
-    "train_period_end = \"2010\"\n",
-    "predict_period_start = \"1981\"\n",
-    "predict_period_end = \"2099\"\n",
-    "latmin = \"40\"\n",
-    "latmax = \"50\"\n",
-    "lonmin = \"80\"\n",
-    "lonmax = \"90\""
+    "variable = \"tasmax\"\n",
+    "latmin = \"-90\"\n",
+    "latmax = \"90\"\n",
+    "lonmin = \"-180\"\n",
+    "lonmax = \"180\"\n",
+    "train_dates = [\"1990\", \"1990\"]\n",
+    "predict_dates = [\"2020\", \"2020\"]\n",
+    "run_id = \"bcsd_ERA5_MIROC6_ssp370_tasmax_-90_90_-180_180_1990_1990_2020_2020\""
    ]
   },
   {
@@ -103,7 +104,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "\"/\".join([result_dir, make_bcsd_output_path(gcm_identifier)])"
+    "intermediate_dir = UPath(config.get(\"storage.intermediate.uri\"))\n",
+    "results_dir = UPath(config.get(\"storage.results.uri\"))\n",
+    "daily_path = results_dir / \"daily\" / run_id\n",
+    "print(intermediate_dir)\n",
+    "print(results_dir)\n",
+    "print(daily_path)"
    ]
   },
   {
@@ -133,7 +139,7 @@
    },
    "outputs": [],
    "source": [
-    "ds = xr.open_zarr(\"/\".join([result_dir, make_bcsd_output_path(gcm_identifier)]))"
+    "ds = xr.open_zarr(daily_path)"
    ]
   },
   {
@@ -157,8 +163,8 @@
     "    lonmin=float(lonmin),\n",
     "    lonmax=float(lonmax),\n",
     ")\n",
-    "train_period = slice(train_period_start, train_period_end)\n",
-    "future_period = slice(predict_period_start, predict_period_end)"
+    "train_period = slice(train_dates[0], train_dates[1])\n",
+    "future_period = slice(predict_dates[0], predict_dates[1])"
    ]
   },
   {
@@ -181,9 +187,9 @@
    "source": [
     "obs_ds = open_era5(\n",
     "    time_period=train_period,\n",
-    "    variables=[var],\n",
+    "    variables=variable,\n",
     ").sel(lat=lat_slice, lon=lon_slice)\n",
-    "if var == \"pr\":\n",
+    "if variable == \"pr\":\n",
     "    obs_ds *= (\n",
     "        86400  # can remove this once ERA5 is rescaled in open_era5 function\n",
     "    )"
@@ -205,12 +211,12 @@
    "outputs": [],
    "source": [
     "gcm_ds = get_gcm(\n",
-    "    gcm,\n",
-    "    scenario,\n",
-    "    variables=[var],\n",
-    "    train_period=train_period,\n",
-    "    predict_period=future_period,\n",
-    "    bbox=bbox,\n",
+    "    scenario=scenario,\n",
+    "    member_id=member,\n",
+    "    table_id=table_id,\n",
+    "    variable=variable,\n",
+    "    source_id=model,\n",
+    "    grid_label=grid_label,\n",
     ")  # .sel(lat=lat_slice, lon=lon_slice)"
    ]
   },
@@ -811,9 +817,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "notebook-kernel",
+   "display_name": "Python [conda env:notebook] *",
    "language": "python",
-   "name": "notebook-kernel"
+   "name": "conda-env-notebook-py"
   },
   "language_info": {
    "codemirror_mode": {
diff --git a/cmip6_downscaling/analysis/analysis.py b/cmip6_downscaling/analysis/analysis.py
index 44150612..a23dd9fa 100644
--- a/cmip6_downscaling/analysis/analysis.py
+++ b/cmip6_downscaling/analysis/analysis.py
@@ -1,23 +1,12 @@
-import os
 from typing import Tuple
 
 import cartopy.crs as ccrs
-import fsspec
 import matplotlib.pyplot as plt
 import pandas as pd
-import papermill as pm
 import xarray as xr
-from azure.storage.blob import BlobServiceClient, ContentSettings
-from prefect import task
-
-from cmip6_downscaling.workflows.paths import get_notebook_paths
-from cmip6_downscaling.workflows.utils import BBox
 
 from .qaqc import make_qaqc_ds
 
-connection_string = os.environ.get("AZURE_STORAGE_CONNECTION_STRING")
-fs = fsspec.filesystem('az', connection_string=connection_string)
-
 
 def qaqc_checks(ds: xr.Dataset) -> Tuple[xr.Dataset, xr.Dataset]:
     '''
@@ -41,125 +30,6 @@ def qaqc_checks(ds: xr.Dataset) -> Tuple[xr.Dataset, xr.Dataset]:
     return annual_qaqc_ts, qaqc_maps
 
 
-def monthly_summary(ds: xr.Dataset, **kwargs) -> xr.Dataset:
-    """Creates an monthly summary dataset.
-
-    Parameters
-    ----------
-    ds : xr.Dataset
-        Downscaled dataset at daily timestep
-
-    Returns
-    ----------
-    xr.Dataset
-        Downscaled dataset at monthly timestep
-    """
-    out_ds = xr.Dataset()
-    for var in ds:
-        if var in ['tasmax', 'tasmin']:
-            out_ds[var] = ds[var].resample(time='1MS').mean(dim='time')
-        elif var in ['pr']:
-            out_ds[var] = ds[var].resample(time='1MS').sum(dim='time')
-        else:
-            print(f'{var} not implemented')
-
-    return out_ds
-
-
-def annual_summary(ds: xr.Dataset, **kwargs) -> xr.Dataset:
-    """Creates an annual summary dataset.
-
-    Parameters
-    ----------
-    ds : xr.Dataset
-        Downscaled dataset at daily timestep
-
-    Returns
-    ----------
-    xr.Dataset
-        Downscaled dataset at annual timestep
-    """
-    out_ds = xr.Dataset()
-    for var in ds:
-        if var in ['tasmax', 'tasmin']:
-            out_ds[var] = ds[var].resample(time='YS').mean()
-        elif var in ['pr']:
-            out_ds[var] = ds[var].resample(time='YS').sum()
-        else:
-            print(f'{var} not implemented')
-
-    return out_ds
-
-
-@task(log_stdout=True, tags=['dask-resource:TASKSLOTS=1'])
-def run_analyses(
-    parameters: dict, bbox: BBox, train_period: slice, predict_period: slice, web_blob: str
-) -> str:
-    """Prefect task to run the analyses on results from a downscaling run.
-
-    Parameters
-    ----------
-    parameters : dict
-        Parameters provided from the downscaling run to uniquely identify
-        the output dataset (e.g. gcm, training period, obs dataset, etc.)
-    bbox : BBox
-        The bounding box over which you are running.
-    train_period : slice
-        Training period (e.g. slice('1981', '2010'))
-    predict_period : slice
-        Prediction period (e.g. slice('1981', '2099')). Currently this assumes
-        that the downscaled dataset has the historical and future timeseries all
-        in one file.
-    web_blob : str
-        The location in the web bucket where you want these files to go. Likely
-        specified in a yaml file in your local setup.
-
-    Returns
-    -------
-    str
-        The local location of an executed notebook path.
-    """
-    gcm_identifier = parameters['gcm_identifier']
-    template_path, executed_notebook_path, executed_html_path = get_notebook_paths(
-        gcm_identifier.replace('/', '_')
-    )
-
-    parameters['train_period_start'] = train_period.start
-    parameters['train_period_end'] = train_period.stop
-    parameters['predict_period_start'] = predict_period.start
-    parameters['predict_period_end'] = predict_period.stop
-    parameters['latmax'] = bbox.lat_slice.start
-    parameters['latmin'] = bbox.lat_slice.stop
-    parameters['lonmin'] = bbox.lon_slice.start
-    parameters['lonmax'] = bbox.lon_slice.stop
-    pm.execute_notebook(template_path, executed_notebook_path, parameters=parameters)
-    # convert from ipynb to html
-    os.system(f"jupyter nbconvert {executed_notebook_path} --to html")
-    connection_string = os.getenv('AZURE_STORAGE_CONNECTION_STRING', None)
-    if connection_string is not None:
-        # if you have a connection_string, copy the html to azure, if not just return
-        # because it is already in your local machine
-        blob_service_client = BlobServiceClient.from_connection_string(connection_string)
-        # TODO: fix b/c the run_id has slashes now!!!
-        blob_name = web_blob + gcm_identifier + 'analyses.html'
-        blob_client = blob_service_client.get_blob_client(container='$web', blob=blob_name)
-        # clean up before writing
-        try:
-            blob_client.delete_blob()
-        except:
-            pass
-
-        #  need to specify html content type so that it will render and not download
-        with open(executed_html_path, "rb") as data:
-            blob_client.upload_blob(
-                data, content_settings=ContentSettings(content_type='text/html')
-            )
-        print(
-            f'**** Your notebook is hosted here! *****\nhttps://cmip6downscaling.z6.web.core.windows.net/{blob_name}'
-        )
-    return executed_notebook_path
-
-
 def load_top_cities(
     num_cities: int = 100, plot: bool = False, add_additional_cities: bool = True
 ) -> pd.DataFrame:
diff --git a/cmip6_downscaling/analysis/plot.py b/cmip6_downscaling/analysis/plot.py
index a7585070..fcefd76c 100644
--- a/cmip6_downscaling/analysis/plot.py
+++ b/cmip6_downscaling/analysis/plot.py
@@ -6,17 +6,6 @@
 import xarray as xr
 from carbonplan import styles
 
-from cmip6_downscaling.workflows.paths import (
-    make_bcsd_output_path,
-    make_bias_corrected_path,
-    make_coarse_obs_path_full_space,
-    make_coarse_obs_path_full_time,
-    make_gcm_predict_path,
-    make_rechunked_gcm_path,
-    make_return_obs_path,
-    make_spatial_anomalies_path,
-)
-
 styles.mpl.set_theme(style='carbonplan_light')
 
 
@@ -247,67 +236,3 @@ def plot_seasonal(ds1: xr.Dataset, ds2: xr.Dataset) -> mpl.figure.Figure:
     plt.tight_layout()
     plt.close()
     return fig
-
-
-def plot_each_step_bcsd(
-    gcm_identifier: str,
-    obs_identifier: str,
-    gcm_grid_spec: str,
-    result_dir: str,
-    intermediate_dir: str,
-    train_period: slice,
-    var: str,
-) -> mpl.figure.Figure:
-    """Plot the training period mean of each intermediary and output file in the
-     bcsd process. For the spatial anomalies it just plots the first month.
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        unique identifier for a run
-    obs_identifier : str
-        unique identifier for the obs used in training
-    result_dir : str
-        Location of the final results
-    intermediate_dir : str
-        Location of intermediate files
-    train_period : slice
-        Period used for training
-    var : str
-        Variable of interest
-
-    Returns
-    -------
-    mpl.figure.Figure
-        Figure
-    """
-    steps = [
-        make_return_obs_path(obs_identifier),
-        make_coarse_obs_path_full_time(obs_identifier=obs_identifier, gcm_grid_spec=gcm_grid_spec),
-        make_coarse_obs_path_full_space(obs_identifier=obs_identifier, gcm_grid_spec=gcm_grid_spec),
-        make_spatial_anomalies_path(obs_identifier=obs_identifier, gcm_grid_spec=gcm_grid_spec),
-        make_rechunked_gcm_path(gcm_identifier),
-        make_gcm_predict_path(gcm_identifier),
-        make_bias_corrected_path(gcm_identifier),
-        make_bcsd_output_path(gcm_identifier),
-    ]
-
-    fig, axarr = plt.subplots(ncols=len(steps), figsize=(20, 3))
-    for i, path in enumerate(steps):
-        prefix = path.split('/')[0]
-        print(prefix)
-
-        if prefix == 'bcsd_output':
-            data_location = result_dir
-        else:
-            data_location = intermediate_dir
-        print('/'.join([data_location, path]))
-        ds = xr.open_zarr('/'.join([data_location, path]))
-        if prefix == 'spatial_anomalies':
-            ds[var].isel(month=0).plot(ax=axarr[i])
-        else:
-            ds[var].sel(time=train_period).mean(dim='time').plot(ax=axarr[i])
-        axarr[i].set_title(prefix)
-    plt.tight_layout()
-    plt.close()
-    return fig
diff --git a/cmip6_downscaling/config.py b/cmip6_downscaling/config.py
index 53c49558..b58f3231 100644
--- a/cmip6_downscaling/config.py
+++ b/cmip6_downscaling/config.py
@@ -1,5 +1,6 @@
 """Config file used by donfig"""
 _defaults = {
+    'chunk_dims': {'full_space': ('time',), 'full_time': ('lat', 'lon')},
     'storage': {
         'top_level': {
             'uri': 'az://',
@@ -19,6 +20,10 @@
             'uri': 'az://flow-outputs/temporary',
             'storage_options': {"connection_string": "$AZURE_STORAGE_CONNECTION_STRING"},
         },
+        'scratch': {
+            'uri': 'az://scratch',
+            'storage_options': {"connection_string": "$AZURE_STORAGE_CONNECTION_STRING"},
+        },
         'xpersist_store_name': '/xpersist_metadata_store',
         'xpersist_overwrite': {'mode': 'w'},
         'web_results': {
@@ -62,7 +67,7 @@
             'storage_options': {"account_name": "cmip6downscaling"},
         },
     },
-    'run_options': {'runtime': None, 'cleanup_flag': True},
+    'run_options': {'runtime': None, 'cleanup_flag': True, 'use_cache': True},
     "runtime": {
         "cloud": {
             "storage_prefix": "az://",
diff --git a/cmip6_downscaling/data/cmip.py b/cmip6_downscaling/data/cmip.py
index 4537bee1..e96c8067 100644
--- a/cmip6_downscaling/data/cmip.py
+++ b/cmip6_downscaling/data/cmip.py
@@ -1,42 +1,68 @@
-from typing import List, Optional, Union
-
 import dask
-import numpy as np
 import xarray as xr
-import zarr
-
-from cmip6_downscaling import config
-from cmip6_downscaling.workflows.paths import make_rechunked_gcm_path
-from cmip6_downscaling.workflows.utils import (
-    lon_to_180,
-    rechunk_zarr_array_with_caching,
-    subset_dataset,
-)
 
 from . import cat
+from .utils import lon_to_180
+
 
+def postprocess(ds: xr.Dataset) -> xr.Dataset:
+    """Post process input experiment
 
-def maybe_drop_band_vars(ds):
+    - Drops band variables (if present)
+    - Drops height variable (if present)
+    - Squeezes length 1 dimensions (if present)
+    - Standardizes longitude convention to [-180, 180]
+    - Reorders latitudes to [-90, 90]
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        Input dataset
+
+    Returns
+    -------
+    ds : xr.Dataset
+        Post processed dataset
+    """
+
+    # drop band variables
     if ('lat_bnds' in ds) or ('lat_bnds' in ds.coords):
         ds = ds.drop('lat_bnds')
     if ('lon_bnds' in ds) or ('lon_bnds' in ds.coords):
         ds = ds.drop('lon_bnds')
     if ('time_bnds' in ds) or ('time_bnds' in ds.coords):
         ds = ds.drop('time_bnds')
+
+    # drop height variable
+    if 'height' in ds:
+        ds = ds.drop('height')
+
+    # squeeze length 1 dimensions
+    ds = ds.squeeze(drop=True)
+
+    # standardize longitude convention
+    ds = lon_to_180(ds)
+
+    # Reorders latitudes to [-90, 90]
+    if ds.lat[0] > ds.lat[-1]:
+        ds = ds.reindex({"lat": ds.lat[::-1]})
+
+    # Shifts time from Noon (12:00) start to Midnight (00:00) start to match with Obs
+    ds['time'] = ds['time'].resample(time='1D').first()
     return ds
 
 
 def load_cmip(
-    activity_ids: str = "CMIP",
-    experiment_ids: str = "historical",
-    member_ids: str = "r1i1p1f1",
-    source_ids: str = "MIROC6",
-    table_ids: str = "day",
-    grid_labels: str = "gn",
-    variable_ids: List[str] = ["tasmax"],
-    return_type: str = 'zarr',
+    activity_ids: str = None,
+    experiment_ids: str = None,
+    member_ids: str = None,
+    source_ids: str = None,
+    table_ids: str = None,
+    grid_labels: str = None,
+    variable_ids: list[str] = None,
 ) -> xr.Dataset:
     """Loads CMIP6 GCM dataset based on input criteria.
+
     Parameters
     ----------
     activity_ids : list, optional
@@ -53,213 +79,93 @@ def load_cmip(
         grid_labels in CMIP6 catalog, by default ["gn"]
     variable_ids : list, optional
         variable_ids in CMIP6 catalog, by default ['tasmax']
+
     Returns
     -------
-    ds : xr.Dataset or zarr group
+    ds : xr.Dataset
         Dataset or zarr group with CMIP data
     """
     with dask.config.set(**{'array.slicing.split_large_chunks': False}):
 
-        if isinstance(variable_ids, str):
-            variable_ids = [variable_ids]
-
         col = cat.cmip6()
+        col_subset = col.search(
+            activity_id=activity_ids,
+            experiment_id=experiment_ids,
+            member_id=member_ids,
+            source_id=source_ids,
+            table_id=table_ids,
+            grid_label=grid_labels,
+            variable_id=variable_ids,
+        )
+        keys = list(col_subset.keys())
+        if len(keys) != 1:
+            raise ValueError(f'intake-esm search returned {len(keys)}, expected exactly 1.')
 
-        for i, var in enumerate(variable_ids):
-            stores = (
-                col.search(
-                    activity_id=activity_ids,
-                    experiment_id=experiment_ids,
-                    member_id=member_ids,
-                    source_id=source_ids,
-                    table_id=table_ids,
-                    grid_label=grid_labels,
-                    variable_id=[var],
-                )
-                .df['zstore']
-                .to_list()
-            )
-
-            storage_options = config.get('data_catalog.era5.storage_options')
-            if len(stores) > 1:
-                raise ValueError('can only get 1 store at a time')
-            if return_type == 'zarr':
-                ds = zarr.open_consolidated(stores[0], mode='r', storage_options=storage_options)
-            elif return_type == 'xr':
-                ds = xr.open_zarr(stores[0], consolidated=True, storage_options=storage_options)
-
-            # flip the lats if necessary and drop the extra dims/vars like bnds
-            ds = gcm_munge(ds)
-            ds = lon_to_180(ds)
-
-            # convert to mm/day - helpful to prevent rounding errors from very tiny numbers
-            if var == 'pr':
-                ds['pr'] *= 86400
+        ds = col_subset[keys[0]].to_dask().pipe(postprocess)
 
-            if i == 0:
-                ds_out = ds
-            else:
-                ds_out[var] = ds[var]
+        # convert to mm/day - helpful to prevent rounding errors from very tiny numbers
+        if 'pr' in ds:
+            ds['pr'] *= 86400
 
-        return ds_out
-
-
-def convert_to_360(lon: Union[float, int]) -> Union[float, int]:
-    """Convert lons to 0-360 basis.
-    Parameters
-    ----------
-    lon : float or int
-        Longitude on -180 to 180 basis
-    Returns
-    -------
-    lon : float or int
-        Longitude on 0 to 360 basis
-    """
-    if lon > 0:
-        return lon
-    elif lon < 0:
-        return 360 + lon
-
-
-def gcm_munge(ds: xr.Dataset) -> xr.Dataset:
-    """Clean up GCM dataset by swapping lats if necessary to match ERA5 and
-    deleting unnecessary variables (e.g. height).
-    Parameters
-    ----------
-    ds : xr.Dataset
-        GCM dataset direct from catalog (though perhaps subsetted temporally)
-    Returns
-    -------
-    ds : xr.Dataset
-        Super clean GCM dataset
-    """
-    # TODO: check if we need to flip this to > now that we have a preprocessed version of ERA5
-    # TODO: for other gcm grids check the lons
-    if ds.lat[0] < ds.lat[-1]:
-        ds = ds.reindex({"lat": ds.lat[::-1]})
-    ds = maybe_drop_band_vars(ds)
-    if 'height' in ds:
-        ds = ds.drop('height')
-    ds = ds.squeeze(drop=True)
-    return ds
+        return ds
 
 
 def get_gcm(
-    gcm: str,
     scenario: str,
-    variables: Union[str, List[str]],
-    train_period: slice,
-    predict_period: slice,
-    bbox,
-    chunking_approach: Optional[str] = None,
-    cache_within_rechunk: Optional[bool] = True,
+    member_id: str,
+    table_id: str,
+    grid_label: str,
+    source_id: str,
+    variable: str,
 ) -> xr.Dataset:
     """
     Load and combine historical and future GCM into one dataset.
+
     Parameters
     ----------
-    gcm: str
-        Name of GCM
-    scenario: str
+    scenario : str
         Name of scenario
-    variables: str or list
-        Name of variable(s) to load
-    train_period_start: str
-        Start year of train/historical period
-    train_period_end: str
-        End year of train/historical period
-    predict_period_start: str
-        Start year of predict/future period
-    predict_period_end: str
-        End year of predict/future period
-    chunking_approach: Optional[str]
-        'full_space', 'full_time', or None
+    member_id : str
+        Name of member ID
+    table_id : str
+        Name of table ID
+    grid_label : str
+        Name of grid_label
+    source_id : str
+        Name of source_id
+    variable : str
+        Name of variable to load
+    bbox : BBox
+        Bounding box for subset
+
     Returns
     -------
-    ds_gcm: xr.Dataset
+    ds_gcm : xr.Dataset
         A dataset containing both historical and future period of GCM data
     """
+
     historical_gcm = load_cmip(
         activity_ids='CMIP',
         experiment_ids='historical',
-        source_ids=gcm,
-        variable_ids=variables,
-        return_type='xr',
+        member_ids=member_id,
+        table_ids=table_id,
+        grid_labels=grid_label,
+        source_ids=source_id,
+        variable_ids=variable,
     )
 
     future_gcm = load_cmip(
         activity_ids='ScenarioMIP',
         experiment_ids=scenario,
-        source_ids=gcm,
-        variable_ids=variables,
-        return_type='xr',
+        member_ids=member_id,
+        table_ids=table_id,
+        grid_labels=grid_label,
+        source_ids=source_id,
+        variable_ids=variable,
     )
 
     ds_gcm = xr.combine_by_coords([historical_gcm, future_gcm], combine_attrs='drop_conflicts')
 
-    ds_gcm_train = subset_dataset(
-        ds=ds_gcm,
-        variable=variables[0],
-        time_period=train_period,
-        bbox=bbox,
-    )
-    ds_gcm_predict = subset_dataset(
-        ds=ds_gcm,
-        variable=variables[0],
-        time_period=predict_period,
-        bbox=bbox,
-    )
-
-    ds_gcm = xr.combine_by_coords([ds_gcm_train, ds_gcm_predict], combine_attrs='drop_conflicts')
     ds_gcm = ds_gcm.reindex(time=sorted(ds_gcm.time.values))
 
-    if chunking_approach is None:
-        return ds_gcm
-
-    if cache_within_rechunk:
-        path_dict = {
-            'gcm': gcm,
-            'scenario': scenario,
-            'train_period': train_period,
-            'predict_period': predict_period,
-            'variables': variables,
-        }
-        rechunked_path = make_rechunked_gcm_path(chunking_approach=chunking_approach, **path_dict)
-    else:
-        rechunked_path = None
-    ds_gcm_rechunked = rechunk_zarr_array_with_caching(
-        zarr_array=ds_gcm,
-        chunking_approach=chunking_approach,
-        output_path=rechunked_path,
-    )
-    print('ds_gcm_rechunked:')
-    print(ds_gcm_rechunked.chunks)
-
-    return ds_gcm_rechunked
-
-
-def get_gcm_grid_spec(
-    gcm_name: Optional[str] = None, gcm_ds: Optional[Union[xr.Dataset, xr.DataArray]] = None
-) -> str:
-    with dask.config.set(**{'array.slicing.split_large_chunks': False}):
-
-        if gcm_ds is None:
-            """Silences the /srv/conda/envs/notebook/lib/python3.9/site-packages/xarray/core/indexing.py:1228: PerformanceWarning: Slicing is producing a large chunk. To accept the large
-            chunk and silence this warning, set the option >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):"""
-            if gcm_name is not None:
-                raise ValueError('one of gcm_ds or gcm_name has to be not empty')
-            gcm_grid = load_cmip(
-                source_ids=gcm_name,
-                return_type='xr',
-            ).isel(time=0)
-        else:
-            gcm_grid = gcm_ds.isel(time=0)
-
-    nlat = len(gcm_grid.lat)
-    nlon = len(gcm_grid.lon)
-    lat_spacing = int(np.round(abs(gcm_grid.lat[0] - gcm_grid.lat[1]), 1) * 10)
-    lon_spacing = int(np.round(abs(gcm_grid.lon[0] - gcm_grid.lon[1]), 1) * 10)
-    min_lat = int(np.round(gcm_grid.lat.min(), 1))
-    min_lon = int(np.round(gcm_grid.lon.min(), 1))
-
-    return f'{nlat:d}x{nlon:d}_gridsize_{lat_spacing:d}_{lon_spacing:d}_llcorner_{min_lat:d}_{min_lon:d}'
+    return ds_gcm
diff --git a/cmip6_downscaling/data/observations.py b/cmip6_downscaling/data/observations.py
index 2f7196ed..9a64ef6b 100644
--- a/cmip6_downscaling/data/observations.py
+++ b/cmip6_downscaling/data/observations.py
@@ -1,11 +1,9 @@
-from typing import List, Optional, Union
+from typing import List, Union
 
 import xarray as xr
 
-from cmip6_downscaling.workflows.paths import make_rechunked_obs_path
-from cmip6_downscaling.workflows.utils import lon_to_180, rechunk_zarr_array_with_caching
-
 from . import cat
+from .utils import lon_to_180
 
 variable_name_dict = {
     "tasmax": "air_temperature_at_2_metres_1hour_Maximum",
@@ -35,45 +33,15 @@ def open_era5(variables: Union[str, List[str]], time_period: slice) -> xr.Datase
     years = range(int(time_period.start), int(time_period.stop) + 1)
 
     ds = xr.concat([cat.era5(year=year).to_dask()[variables] for year in years], dim='time')
-    ds = lon_to_180(ds)
 
     if 'pr' in variables:
         # convert to mm/day - helpful to prevent rounding errors from very tiny numbers
         ds['pr'] *= 86400
-    return ds
-
 
-def get_obs(
-    obs: str,
-    train_period: slice,
-    variables: Union[str, List[str]],
-    chunking_approach: Optional[str] = None,
-    cache_within_rechunk: Optional[bool] = True,
-) -> xr.Dataset:
-    if obs == 'ERA5':
-        ds_obs = open_era5(variables=variables, time_period=train_period)
-    else:
-        raise NotImplementedError('only ERA5 is available as observation dataset right now')
-
-    if chunking_approach is None:
-        return ds_obs
+    ds = lon_to_180(ds)
 
-    if cache_within_rechunk:
-        path_dict = {
-            'obs': obs,
-            'train_period': train_period,
-            'variables': variables,
-        }
-        rechunked_path = make_rechunked_obs_path(
-            chunking_approach=chunking_approach,
-            **path_dict,
-        )
-    else:
-        rechunked_path = None
-    ds_obs_rechunked = rechunk_zarr_array_with_caching(
-        zarr_array=ds_obs,
-        chunking_approach=chunking_approach,
-        output_path=rechunked_path,
-    )
+    # Reorders latitudes to [-90, 90]
+    if ds.lat[0] > ds.lat[-1]:
+        ds = ds.reindex({"lat": ds.lat[::-1]})
 
-    return ds_obs_rechunked
+    return ds
diff --git a/cmip6_downscaling/data/utils.py b/cmip6_downscaling/data/utils.py
new file mode 100644
index 00000000..de3d7b2c
--- /dev/null
+++ b/cmip6_downscaling/data/utils.py
@@ -0,0 +1,76 @@
+import numpy as np
+import xclim
+from xarray.core.types import T_Xarray
+
+
+def to_standard_calendar(obj: T_Xarray) -> T_Xarray:
+    """Convert a Dataset's calendar to the "standard calendar"
+
+    When necessary, "missing" time points are filled in using linear interpolation.
+
+    Valid input dataset calendars include: `noleap`, `365_day`, `366_day`, and `all_leap`.
+
+    Parameters
+    ----------
+    obj : xr.Dataset or xr.DataArray
+        Xarray object with a `CFTimeIndex`.
+
+    Returns
+    -------
+    obj_new : xr.Dataset or xr.DataArray
+        Xarray object with standard calendar.
+
+    Raises
+    ------
+    ValueError
+        If an invalid calendar is supplied.
+    """
+
+    orig_calendar = getattr(obj.indexes["time"], "calendar", "standard")
+    if orig_calendar == "standard":
+        return obj
+    if orig_calendar == "360_day":
+        raise ValueError("360_day calendar is not supported")
+
+    # reindex / interpolate
+    obj_new = xclim.core.calendar.convert_calendar(obj, "standard", missing=np.nan).interpolate_na(
+        dim="time", method="linear"
+    )
+
+    # reset encoding
+    obj_new["time"].encoding["calendar"] = "standard"
+
+    return obj_new
+
+
+def lon_to_180(ds):
+    '''Converts longitude values to (-180, 180)
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        Input dataset with `lon` coordinate
+
+    Returns
+    -------
+    xr.Dataset
+        Copy of `ds` with updated coordinates
+
+    See also
+    --------
+    cmip6_preprocessing.preprocessing.correct_lon
+    '''
+
+    ds = ds.copy()
+
+    lon = ds["lon"].where(ds["lon"] < 180, ds["lon"] - 360)
+    ds = ds.assign_coords(lon=lon)
+
+    if not (ds["lon"].diff(dim="lon") > 0).all():
+        ds = ds.reindex(lon=np.sort(ds["lon"].data))
+
+    if "lon_bounds" in ds.variables:
+        lon_b = ds["lon_bounds"].where(ds["lon_bounds"] < 180, ds["lon_bounds"] - 360)
+        ds = ds.assign_coords(lon_bounds=lon_b)
+
+    return ds
diff --git a/cmip6_downscaling/disagg/terraclimate.py b/cmip6_downscaling/disagg/terraclimate.py
index 69bb2e1c..5d230353 100644
--- a/cmip6_downscaling/disagg/terraclimate.py
+++ b/cmip6_downscaling/disagg/terraclimate.py
@@ -147,7 +147,7 @@ def monthly_pet(
     delta = (4098 * es) / (tmean + 237.3) ** 2
 
     lhv = 2.501 - 2.361e-3 * tmean  # latent heat of vaporization
-    cp = 1.013 * 10 ** -3  # specific heat of air
+    cp = 1.013 * 10**-3  # specific heat of air
 
     gamma = cp * P / (0.622 * lhv)  # Psychrometer constant (kPa C-1)
     pa = P / (1.01 * (tmean + KELVIN) * 0.287)  # mean air density at constant pressure
diff --git a/cmip6_downscaling/methods/bcsd.py b/cmip6_downscaling/methods/bcsd.py
deleted file mode 100644
index a67fea2f..00000000
--- a/cmip6_downscaling/methods/bcsd.py
+++ /dev/null
@@ -1,620 +0,0 @@
-from typing import List, Union
-
-import xarray as xr
-from skdownscale.pointwise_models import PointWiseDownscaler
-from skdownscale.pointwise_models.bcsd import BcsdPrecipitation, BcsdTemperature
-
-from cmip6_downscaling import config
-from cmip6_downscaling.constants import ABSOLUTE_VARS, RELATIVE_VARS
-from cmip6_downscaling.data.cmip import get_gcm, load_cmip
-from cmip6_downscaling.data.observations import open_era5
-from cmip6_downscaling.workflows.paths import (
-    make_bias_corrected_path,
-    make_coarse_obs_path,
-    make_gcm_predict_subset_path,
-    make_gcm_train_subset_path,
-    make_interpolated_obs_path,
-    make_interpolated_prediction_path_full_space,
-    make_interpolated_prediction_path_full_time,
-    make_return_obs_path,
-    make_spatial_anomalies_path,
-)
-from cmip6_downscaling.workflows.utils import (
-    BBox,
-    lon_to_180,
-    rechunk_zarr_array_with_caching,
-    regrid_ds,
-    subset_dataset,
-)
-
-
-def return_obs(
-    obs: str,
-    variable: Union[str, List[str]],
-    train_period: slice,
-    bbox: BBox,
-    obs_identifier: str,
-    **kwargs
-) -> xr.Dataset:
-    """Loads ERA5 observation data for given time bounds and variable
-
-    Parameters
-    ----------
-    obs : str
-        Input obs
-    variable: str
-        The variable included in the dataset.
-    train_period: slice
-        Start and end year slice of training/historical period. Ex: slice('1990','1990')
-    predict_period: slice
-        Start and end year slice of predict period. Ex: slice('2020','2020')
-    bbox : BBox
-        Bounding box including latmin,latmax,lonmin,lonmax.
-    **kwargs : dict, optional
-            Other arguments to be used in generating the target path
-
-    Returns
-    -------
-    xr.Dataset
-        Loaded xarray dataset of ERA5 observation data. Chunked in time: 365
-    """
-    if isinstance(variable, str):
-        variable = [variable]
-
-    obs_load = open_era5(variable, train_period)
-    obs_load_180 = lon_to_180(obs_load)
-    obs_ds = subset_dataset(
-        obs_load_180,
-        variable[0],
-        train_period,
-        bbox,
-        chunking_schema={'time': 365, 'lat': 150, 'lon': 150},
-    )
-    del obs_ds[variable[0]].encoding['chunks']
-    return obs_ds
-
-
-def get_coarse_obs(
-    gcm: str,
-    scenario: str,
-    variable: Union[str, List[str]],
-    train_period: slice,
-    predict_period: slice,
-    bbox: BBox,
-    obs_identifier: str,
-    gcm_grid_spec: str,
-    chunking_approach: str,
-    **kwargs
-) -> xr.Dataset:
-    """Regrids the observation dataset to match the GCM resolution
-
-    Parameters
-    ----------
-
-    gcm : str
-        Input GCM
-    scenario: str
-        Input GCM scenario
-    variable: str
-        The variable included in the dataset.
-    train_period: slice
-        Start and end year slice of training/historical period. Ex: slice('1990','1990')
-    predict_period: slice
-        Start and end year slice of predict period. Ex: slice('2020','2020')
-    bbox : BBox
-        Bounding box including latmin,latmax,lonmin,lonmax.
-
-    **kwargs : dict, optional
-            Other arguments to be used in generating the target path
-
-    Returns
-    -------
-    xr.Dataset
-        observation dataset at coarse resolution
-    """
-    # Load single slice of target cmip6 dataset for target grid dimensions
-    # gcm_one_slice = load_cmip(return_type='xr', variable_ids=[variable]).isel(time=0)
-
-    obs_ds_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_return_obs_path(obs_identifier=obs_identifier)
-    )
-    if isinstance(variable, str):
-        variable = [variable]
-
-    gcm_ds = load_cmip(source_ids=gcm, return_type='xr', variable_ids=variable)
-
-    gcm_ds_180 = lon_to_180(gcm_ds)
-    gcm_subset = subset_dataset(gcm_ds_180, variable[0], train_period, bbox)
-
-    # rechunk and regrid observation dataset to target gcm resolution
-    coarse_obs_ds = regrid_ds(ds_path=obs_ds_path, target_grid_ds=gcm_subset)
-
-    return coarse_obs_ds
-
-
-def get_interpolated_obs(
-    target_grid_ds: xr.Dataset, gcm_grid_spec: str, chunking_approach: str, obs_identifier: str
-):
-
-    coarse_obs_ds_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_return_obs_path(obs_identifier=obs_identifier)
-    )
-    interpolated_obs = regrid_ds(ds_path=coarse_obs_ds_path, target_grid_ds=target_grid_ds)
-    return interpolated_obs
-
-
-def get_spatial_anomalies(
-    gcm: str,
-    scenario: str,
-    variable: str,
-    train_period: slice,
-    predict_period: slice,
-    bbox: BBox,
-    obs_identifier: str,
-    gcm_grid_spec: str,
-    chunking_approach: str,
-    **kwargs
-) -> xr.Dataset:
-
-    """Returns spatial anomalies
-    Calculate the seasonal cycle (12 timesteps) spatial anomaly associated
-    with aggregating the fine_obs to a given coarsened scale and then reinterpolating
-    it back to the original spatial resolution. The outputs of this function are
-    dependent on three parameters:
-    * a grid (as opposed to a specific GCM since some GCMs run on the same grid)
-    * the time period which fine_obs (and by construct coarse_obs) cover
-    * the variable
-    We will save these anomalies to use them in the post-processing. We will add them to the
-    spatially-interpolated coarse predictions to add the spatial heterogeneity back in.
-    Conceptually, this step figures out, for example, how much colder a finer-scale pixel
-    containing Mt. Rainier is compared to the coarse pixel where it exists. By saving those anomalies,
-    we can then preserve the fact that "Mt Rainier is x degrees colder than the pixels around it"
-    for the prediction. It is important to note that that spatial anomaly is the same for every month of the
-    year and the same for every day. So, if in January a finescale pixel was on average 4 degrees colder than
-    the neighboring pixel to the west, in every day in the prediction (historic or future) that pixel
-    will also be 4 degrees colder.
-
-    Parameters
-    ----------
-    gcm : str
-        Input GCM
-    scenario: str
-        Input GCM scenario
-    variable: str
-        The variable included in the dataset.
-    train_period: slice
-        Start and end year slice of training/historical period. Ex: slice('1990','1990')
-    predict_period: slice
-        Start and end year slice of predict period. Ex: slice('2020','2020')
-    bbox : BBox
-        Bounding box including latmin,latmax,lonmin,lonmax.
-    **kwargs : dict, optional
-            Other arguments to be used in generating the target path
-
-    Returns
-    -------
-    seasonal_cycle_spatial_anomalies : xr.Dataset
-        Spatial anomaly for each month (i.e. of shape (nlat, nlon, 12))
-    """
-
-    interpolated_obs_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_interpolated_obs_path(
-            gcm_grid_spec=gcm_grid_spec,
-            chunking_approach=chunking_approach,
-            obs_identifier=obs_identifier,
-        )
-    )
-    # take interpolated obs and rechunk into full_time -- returns path
-    coarse_obs_interpolated_rechunked_path = rechunk_zarr_array_with_caching(
-        interpolated_obs_path, chunking_approach='full_time', max_mem='2GB'
-    )
-
-    obs_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_return_obs_path(obs_identifier=obs_identifier)
-    )
-    # original obs rechunked into full-time
-    obs_rechunked_path = rechunk_zarr_array_with_caching(
-        obs_path, chunking_approach='full_time', max_mem='2GB'
-    )
-    coarse_obs_interpolated_rechunked_ds = xr.open_zarr(coarse_obs_interpolated_rechunked_path)
-    obs_rechunked_ds = xr.open_zarr(obs_rechunked_path)
-
-    # calculate the difference between the actual obs (with finer spatial heterogeneity)
-    # and the interpolated coarse obs this will be saved and added to the
-    # spatially-interpolated coarse predictions to add the spatial heterogeneity back in.
-
-    spatial_anomalies = obs_rechunked_ds - coarse_obs_interpolated_rechunked_ds
-    seasonal_cycle_spatial_anomalies = spatial_anomalies.groupby("time.month").mean()
-
-    return seasonal_cycle_spatial_anomalies
-
-
-def return_coarse_obs_full_time(
-    gcm: str,
-    scenario: str,
-    variable: str,
-    train_period: slice,
-    predict_period: slice,
-    bbox: BBox,
-    obs_identifier: str,
-    gcm_grid_spec: str,
-    **kwargs
-) -> xr.Dataset:
-
-    """
-
-    Return coarse observation dataset that has been chunked in time.
-
-    Parameters
-    ----------
-    gcm : str
-        Input GCM
-    scenario: str
-        Input GCM scenario
-    variable: str
-        The variable included in the dataset.
-    train_period: slice
-        Start and end year slice of training/historical period. Ex: slice('1990','1990')
-    predict_period: slice
-        Start and end year slice of predict period. Ex: slice('2020','2020')
-    bbox : BBox
-        Bounding box including latmin,latmax,lonmin,lonmax.
-    **kwargs : dict, optional
-
-    Returns
-    -------
-    xr.Dataset
-        coarse_obs_full_time_ds rechunked dataset
-    """
-    coarse_obs_path_full_space_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_coarse_obs_path(
-            obs_identifier=obs_identifier,
-            gcm_grid_spec=gcm_grid_spec,
-            chunking_approach='full_space',
-        )
-    )
-    # this is temp_path: az://flow-outputs/temporary/djnqgjlkbr.zarr. this is target_path: az://flow-outputs/test_intermediate_zarr/az://flow-outputs/test_intermediate_zarr/coarsened_obs/ERA5/tasmax/-90.0_90.0_-180.0_180.0/1990_1990/full_time_128x256_gridsize_14_14_llcorner_-88_-180.zarr
-    #    az://flow-outputs/test_intermediate_zarr/coarsened_obs/ERA5/tasmax/-90.0_90.0_-180.0_180.0/1990_1990/full_time_128x256_gridsize_14_14_llcorner_-88_-180.zarr
-    coarse_obs_path_full_time_path = make_coarse_obs_path(
-        obs_identifier=obs_identifier, gcm_grid_spec=gcm_grid_spec, chunking_approach='full_time'
-    )
-
-    coarse_obs_full_time_path = rechunk_zarr_array_with_caching(
-        coarse_obs_path_full_space_path,
-        output_path=coarse_obs_path_full_time_path,
-        chunking_approach='full_time',
-        max_mem='2GB',
-    )
-    return coarse_obs_full_time_path
-
-
-def return_gcm_train_full_time(
-    coarse_obs_full_time_path: str,
-    gcm: str,
-    scenario: str,
-    variable: str,
-    train_period: slice,
-    predict_period: slice,
-    bbox: BBox,
-    gcm_identifier: str,
-    **kwargs
-) -> xr.Dataset:
-    """Returns GCM training rechunked dataset in full time.
-
-    Parameters
-    ----------
-    coarse_obs_full_time_path : str
-        Output coarse observation dataset path rechunked in full_time
-    gcm : str
-        Input GCM
-    scenario: str
-        Input GCM scenario
-    variable: str
-        The variable included in the dataset.
-    train_period: slice
-        Start and end year slice of training/historical period. Ex: slice('1990','1990')
-    predict_period: slice
-        Start and end year slice of predict period. Ex: slice('2020','2020')
-    bbox : BBox
-        Bounding box including latmin,latmax,lonmin,lonmax.
-    **kwargs : dict, optional
-
-    Returns
-    -------
-    xr.Dataset
-        x_train rechunked dataset in full time.
-    """
-    gcm_train_ds = load_cmip(source_ids=gcm, variable_ids=[variable], return_type='xr')
-    gcm_train_ds_180 = lon_to_180(gcm_train_ds)
-
-    gcm_train_ds_subset = subset_dataset(
-        gcm_train_ds_180,
-        variable,
-        train_period,
-        bbox,
-        # is the chunking schema needed?
-        chunking_schema={'time': 365, 'lat': 150, 'lon': 150},
-    )
-    del gcm_train_ds_subset[variable].encoding['chunks']
-
-    # this call was to force the timestamps for the cmip data to use the friendlier era5 timestamps. (i forget which dataset used which time formats). i could picture this introducing a tricky bug though (for instance if gcm timestamp didn't align for some reason) so we could use another conversion system if that is better. Perhaps datetime equivilence test.
-    coarse_obs_full_time_ds = xr.open_zarr(coarse_obs_full_time_path)
-    gcm_train_ds_subset['time'] = coarse_obs_full_time_ds.time.values
-    gcm_train_subset_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_gcm_train_subset_path(gcm_identifier=gcm_identifier)
-    )
-    gcm_train_ds_subset.to_zarr(gcm_train_subset_path, mode='w')
-
-    gcm_train_subset_full_time_path = rechunk_zarr_array_with_caching(
-        gcm_train_subset_path,
-        chunking_approach='full_time',
-        max_mem='2GB',
-    )
-    gcm_train_ds_subset = xr.open_zarr(gcm_train_subset_full_time_path)
-
-    return gcm_train_ds_subset
-
-
-def return_gcm_predict_rechunked(
-    gcm_train_subset_full_time_ds: xr.Dataset,
-    gcm: str,
-    scenario: str,
-    variable: str,
-    train_period: slice,
-    predict_period: slice,
-    bbox: BBox,
-    gcm_identifier: str,
-    **kwargs
-) -> xr.Dataset:
-    """Returns GCM prediction rechunked dataset in full time.  Chunks are matched to chunks of gcm train. In the current use case, this means in full_time.
-
-    Parameters
-    ----------
-    gcm_train_subset_full_time_ds : xr.Dataset
-        Input gcm training rechunked dataset
-    gcm : str
-        Input GCM
-    scenario: str
-        Input GCM scenario
-    variable: str
-        The variable included in the dataset.
-    train_period: slice
-        Start and end year slice of training/historical period. Ex: slice('1990','1990')
-    predict_period: slice
-        Start and end year slice of predict period. Ex: slice('2020','2020')
-    bbox : BBox
-        Bounding box including latmin,latmax,lonmin,lonmax.
-    **kwargs : dict, optional
-
-    Returns
-    -------
-    xr.Dataset
-        gcm predict rechunked dataset
-    """
-    gcm_predict_ds = get_gcm(
-        gcm=gcm,
-        scenario=scenario,
-        variables=[variable],
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-    )
-
-    gcm_predict_ds_subset = subset_dataset(
-        gcm_predict_ds,
-        variable,
-        predict_period,
-        bbox,
-    )
-
-    del gcm_predict_ds_subset[variable].encoding['chunks']
-
-    gcm_predict_ds_subset_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_gcm_predict_subset_path(gcm_identifier=gcm_identifier)
-    )
-    # Note: Explain why this .chunk is needed!
-    gcm_predict_ds_subset.chunk({'time': 1500}).to_zarr(gcm_predict_ds_subset_path, mode='w')
-
-    gcm_predict_rechunked_path = rechunk_zarr_array_with_caching(
-        gcm_predict_ds_subset_path,
-        # note: predict dataset might need to match the spatial chunking of the training dataset -- time chunking can be differant.
-        template_chunk_array=gcm_train_subset_full_time_ds,
-        max_mem='2GB',
-    )
-    gcm_predict_rechunked_ds = xr.open_zarr(gcm_predict_rechunked_path)
-    return gcm_predict_rechunked_ds
-
-
-def fit_and_predict(
-    gcm_train_subset_full_time_ds: xr.Dataset,
-    gcm_predict_rechunked_ds: xr.Dataset,
-    gcm: str,
-    scenario: str,
-    variable: str,
-    train_period: slice,
-    predict_period: slice,
-    bbox: BBox,
-    chunking_approach: str,
-    obs_identifier: str,
-    gcm_grid_spec: str,
-    gcm_identifier: str,
-    dim: str = "time",
-    **kwargs
-) -> xr.Dataset:
-    """Fit bcsd model on prepared CMIP data with obs at corresponding spatial scale.
-    Then predict for a set of CMIP data (likely future).
-
-    Parameters
-    ----------
-    gcm_train_subset_full_time_ds : xr.Dataset
-        GCM training dataset chunked along space
-    gcm_predict_rechunked_ds : xr.Dataset
-        GCM prediction dataset chunked along space.
-    gcm : str
-        Input GCM
-    scenario: str
-        Input GCM scenario
-    variable: str
-        The variable included in the dataset.
-    train_period: slice
-        Start and end year slice of training/historical period. Ex: slice('1990','1990')
-    predict_period: slice
-        Start and end year slice of predict period. Ex: slice('2020','2020')
-    bbox : BBox
-        Bounding box including latmin,latmax,lonmin,lonmax.
-    dim : str, optional
-        dimension on which you want to do the modelling, by default "time"
-    **kwargs : dict, optional
-
-    Returns
-    -------
-    bias_corrected_ds : xr.Dataset
-        Bias-corrected dataset
-    """
-    if variable in ABSOLUTE_VARS:
-        bcsd_model = BcsdTemperature(return_anoms=False)
-    elif variable in RELATIVE_VARS:
-        bcsd_model = BcsdPrecipitation(return_anoms=False)
-
-    pointwise_model = PointWiseDownscaler(model=bcsd_model, dim=dim)
-    coarse_obs_full_time_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_coarse_obs_path(
-            chunking_approach=chunking_approach,
-            obs_identifier=obs_identifier,
-            gcm_grid_spec=gcm_grid_spec,
-        )
-    )
-    coarse_obs_rechunked_validated_path = rechunk_zarr_array_with_caching(
-        coarse_obs_full_time_path, template_chunk_array=gcm_train_subset_full_time_ds
-    )
-    coarse_obs_rechunked_validated_ds = xr.open_zarr(coarse_obs_rechunked_validated_path)
-
-    pointwise_model.fit(
-        gcm_train_subset_full_time_ds[variable], coarse_obs_rechunked_validated_ds[variable]
-    )
-    bias_corrected_da = pointwise_model.predict(gcm_predict_rechunked_ds[variable])
-
-    bias_corrected_ds = bias_corrected_da.astype('float32').to_dataset(name=variable)
-
-    return bias_corrected_ds
-
-
-def get_interpolated_prediction(
-    target_grid_obs_ds: xr.Dataset, gcm_grid_spec: str, gcm_identifier: str
-) -> xr.Dataset:
-
-    bias_corrected_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_bias_corrected_path(gcm_identifier=gcm_identifier)
-    )
-
-    interpolated_prediction_ds = regrid_ds(
-        ds_path=bias_corrected_path, target_grid_ds=target_grid_obs_ds
-    )
-    return interpolated_prediction_ds
-
-
-def rechunked_interpolated_prediciton_task_full_time(
-    gcm_identifier: str, gcm_grid_spec: str, **kwargs
-) -> str:
-
-    interpolated_prediction_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_interpolated_prediction_path_full_space(
-            gcm_identifier=gcm_identifier, gcm_grid_spec=gcm_grid_spec
-        )
-    )
-    rechunked_interpolated_prediction_path = make_interpolated_prediction_path_full_time(
-        gcm_identifier=gcm_identifier, gcm_grid_spec=gcm_grid_spec
-    )
-
-    rechunked_interpolated_prediction_path = rechunk_zarr_array_with_caching(
-        interpolated_prediction_path,
-        output_path=rechunked_interpolated_prediction_path,
-        chunking_approach='full_time',
-        max_mem='2GB',
-    )
-
-    return rechunked_interpolated_prediction_path
-
-
-def postprocess_bcsd(
-    gcm: str,
-    scenario: str,
-    variable: str,
-    train_period: slice,
-    predict_period: slice,
-    bbox: BBox,
-    obs_identifier: str,
-    gcm_identifier: str,
-    gcm_grid_spec: str,
-    **kwargs
-) -> xr.Dataset:
-    """Downscale the bias-corrected data by interpolating and then
-    adding the spatial anomalies back in.
-
-    Parameters
-    ----------
-    gcm : str
-        Input GCM
-    scenario: str
-        Input GCM scenario
-    variable: str
-        The variable included in the dataset.
-    train_period : slice
-        Start and end year slice of training/historical period. Ex: slice('1990', '1990')
-    predict_period : slice
-        Start and end year slice of predict period. Ex: slice('2020', '2020')
-    bbox : BBox
-        Bounding box including latmin,latmax,lonmin,lonmax.
-    **kwargs : dict, optional
-
-    Returns
-    -------
-    bcsd_results_ds : xr.Dataset
-        Final BCSD dataset
-    """
-
-    spatial_anomalies_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_spatial_anomalies_path(obs_identifier=obs_identifier, gcm_grid_spec=gcm_grid_spec)
-    )
-    rechunked_interpolated_prediction_path = (
-        config.get("storage.intermediate.uri")
-        + '/'
-        + make_interpolated_prediction_path_full_time(
-            gcm_identifier=gcm_identifier, gcm_grid_spec=gcm_grid_spec
-        )
-    )
-
-    spatial_anomalies_ds = xr.open_zarr(spatial_anomalies_path)
-    rechunked_interpolated_prediction_ds = xr.open_zarr(rechunked_interpolated_prediction_path)
-
-    bcsd_results_ds = (
-        rechunked_interpolated_prediction_ds.groupby("time.month") + spatial_anomalies_ds
-    )
-    del bcsd_results_ds['month'].encoding['chunks']
-
-    # failes without .chunk: #ValueError: Zarr requires uniform chunk sizes except for final chunk. Variable named 'tasmax' has incompatible dask chunks: ((31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31), (93, 93, 93, 93, 93, 93, 93, 70), (93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 93, 45)). Consider rechunking using `chunk()`.
-    rechunked_bcsd_results_ds = bcsd_results_ds.chunk({'time': 365, 'lat': -1, 'lon': -1})
-
-    return rechunked_bcsd_results_ds
diff --git a/cmip6_downscaling/methods/bcsd/flow.py b/cmip6_downscaling/methods/bcsd/flow.py
new file mode 100644
index 00000000..972b8984
--- /dev/null
+++ b/cmip6_downscaling/methods/bcsd/flow.py
@@ -0,0 +1,123 @@
+import warnings
+
+from prefect import Flow, Parameter
+
+from . import runtimes
+from .methods.bcsd.tasks import fit_and_predict, postprocess_bcsd, spatial_anomalies
+from .methods.common.tasks import (
+    annual_summary,
+    get_experiment,
+    get_obs,
+    make_run_parameters,
+    monthly_summary,
+    pyramid,
+    rechunk,
+    regrid,
+    run_analyses,
+)
+
+warnings.filterwarnings(
+    "ignore",
+    "(.*) filesystem path not explicitly implemented. falling back to default implementation. This filesystem may not be tested",
+    category=UserWarning,
+)
+
+runtime = runtimes.get_runtime()
+print(runtime)
+
+
+with Flow(
+    name="bcsd", storage=runtime.storage, run_config=runtime.run_config, executor=runtime.executor
+) as flow:
+
+    run_parameters = make_run_parameters(
+        method=Parameter("method"),
+        obs=Parameter("obs"),
+        model=Parameter("model"),
+        member=Parameter("member"),
+        grid_label=Parameter("grid_label"),
+        table_id=Parameter("table_id"),
+        scenario=Parameter("scenario"),
+        variable=Parameter("variable"),
+        latmin=Parameter("latmin"),
+        latmax=Parameter("latmax"),
+        lonmin=Parameter("lonmin"),
+        lonmax=Parameter("lonmax"),
+        train_dates=Parameter("train_period"),
+        predict_dates=Parameter("predict_period"),
+    )
+
+    # input datasets
+    obs_path = get_obs(run_parameters)
+    obs_full_space_path = rechunk(path=obs_path, pattern='full_space')
+    experiment_train_path = get_experiment(run_parameters, time_subset='train_period')
+    experiment_predict_path = get_experiment(run_parameters, time_subset='predict_period')
+
+    # after regridding coarse_obs will have smaller array size in space but still
+    # be chunked finely along time. but that's good to get it for regridding back to
+    # the interpolated obs in next task
+    coarse_obs_path = regrid(obs_full_space_path, experiment_train_path)
+
+    # interpolated obs should have same exact chunking schema as ds at `obs_full_space_path`
+    interpolated_obs_path = regrid(source_path=coarse_obs_path, target_grid_path=obs_path)
+
+    interpolated_obs_full_time_path = rechunk(path=interpolated_obs_path, pattern="full_time")
+    obs_full_time_path = rechunk(path=obs_path, pattern="full_time")
+    spatial_anomalies_path = spatial_anomalies(
+        obs_full_time_path, interpolated_obs_full_time_path, run_parameters
+    )
+    coarse_obs_full_time_path = rechunk(coarse_obs_path, pattern='full_time')
+    experiment_train_full_time_path = rechunk(experiment_train_path, pattern='full_time')
+    experiment_predict_full_time_path = rechunk(
+        experiment_predict_path,
+        pattern='full_time',
+        template=coarse_obs_full_time_path,
+    )
+    bias_corrected_path = fit_and_predict(
+        experiment_train_full_time_path=experiment_train_full_time_path,
+        experiment_predict_full_time_path=experiment_predict_full_time_path,
+        coarse_obs_full_time_path=coarse_obs_full_time_path,
+        run_parameters=run_parameters,
+    )
+    bias_corrected_full_space_path = rechunk(
+        bias_corrected_path,
+        pattern='full_space',
+        template=obs_full_space_path,
+    )
+    bias_corrected_fine_full_space_path = regrid(
+        source_path=bias_corrected_full_space_path, target_grid_path=obs_path
+    )
+
+    bias_corrected_fine_full_time_path = rechunk(
+        bias_corrected_fine_full_space_path,
+        pattern='full_time',
+        template=obs_full_time_path,
+    )
+    final_bcsd_full_time_path = postprocess_bcsd(
+        bias_corrected_fine_full_time_path, spatial_anomalies_path, run_parameters
+    )  # fine-scale maps (full_space) (time: 365)
+
+    # temporary aggregations - these come out in full time
+    monthly_summary_path = monthly_summary(final_bcsd_full_time_path, run_parameters)
+    annual_summary_path = annual_summary(final_bcsd_full_time_path, run_parameters)
+
+    # analysis notebook
+    analysis_location = run_analyses(final_bcsd_full_time_path, run_parameters)
+
+    # since pyramids require full space we now rechunk everything into full
+    # space before passing into pyramid step. we probably want to add a cleanup
+    # to this step in particular since otherwise we will have an exact
+    # duplicate of the daily, monthly, and annual datasets
+    final_bcsd_full_space_path = rechunk(final_bcsd_full_time_path, pattern='full_space')
+
+    # make temporal summaries
+    monthly_summary_full_space_path = rechunk(monthly_summary_path, pattern='full_space')
+    annual_summary_full_space_path = rechunk(annual_summary_path, pattern='full_space')
+
+    # pyramids
+    daily_pyramid_path = pyramid(final_bcsd_full_space_path, levels=4)
+    monthly_pyramid_path = pyramid(monthly_summary_full_space_path, levels=4)
+    annual_pyramid_path = pyramid(annual_summary_full_space_path, levels=4)
+
+    # # if config.get('run_options.cleanup_flag') is True:
+    # #     cleanup.run_rsfip(gcm_identifier, obs_identifier)
diff --git a/cmip6_downscaling/methods/bcsd/tasks.py b/cmip6_downscaling/methods/bcsd/tasks.py
new file mode 100644
index 00000000..4edcd550
--- /dev/null
+++ b/cmip6_downscaling/methods/bcsd/tasks.py
@@ -0,0 +1,201 @@
+import warnings
+from dataclasses import asdict
+
+import xarray as xr
+from prefect import task
+from skdownscale.pointwise_models import PointWiseDownscaler
+from skdownscale.pointwise_models.bcsd import BcsdPrecipitation, BcsdTemperature
+from upath import UPath
+
+import cmip6_downscaling
+
+from . import config
+from .constants import ABSOLUTE_VARS, RELATIVE_VARS
+from .methods.bcsd.utils import reconstruct_finescale
+from .methods.common.containers import RunParameters
+from .methods.common.utils import zmetadata_exists
+
+warnings.filterwarnings(
+    "ignore",
+    "(.*) filesystem path not explicitly implemented. falling back to default implementation. This filesystem may not be tested",
+    category=UserWarning,
+)
+
+
+code_version = cmip6_downscaling.__version__
+scratch_dir = UPath(config.get("storage.scratch.uri"))
+intermediate_dir = UPath(config.get("storage.intermediate.uri")) / cmip6_downscaling.__version__
+results_dir = UPath(config.get("storage.results.uri")) / cmip6_downscaling.__version__
+use_cache = config.get('run_options.use_cache')
+
+
+@task(log_stdout=True)
+def spatial_anomalies(
+    obs_full_time_path: UPath, interpolated_obs_full_time_path: UPath, run_parameters: RunParameters
+) -> UPath:
+    """Returns spatial anomalies
+    Calculate the seasonal cycle (12 timesteps) spatial anomaly associated
+    with aggregating the fine_obs to a given coarsened scale and then reinterpolating
+    it back to the original spatial resolution. The outputs of this function are
+    dependent on three parameters:
+    * a grid (as opposed to a specific GCM since some GCMs run on the same grid)
+    * the time period which fine_obs (and by construct coarse_obs) cover
+    * the variable
+    We will save these anomalies to use them in the post-processing. We will add them to the
+    spatially-interpolated coarse predictions to add the spatial heterogeneity back in.
+    Conceptually, this step figures out, for example, how much colder a finer-scale pixel
+    containing Mt. Rainier is compared to the coarse pixel where it exists. By saving those anomalies,
+    we can then preserve the fact that "Mt Rainier is x degrees colder than the pixels around it"
+    for the prediction. It is important to note that that spatial anomaly is the same for every month of the
+    year and the same for every day. So, if in January a finescale pixel was on average 4 degrees colder than
+    the neighboring pixel to the west, in every day in the prediction (historic or future) that pixel
+    will also be 4 degrees colder.
+
+    Parameters
+    ----------
+    obs_full_time_path : UPath
+        UPath to observation dataset chunked in full_time.
+    interpolated_obs_full_time_path : UPath
+        UPath to observation dataset interpolated to gcm grid and chunked in full time.
+    run_parameters : RunParameters
+        Prefect run parameters
+
+    Returns
+    -------
+    UPath
+        Path to spatial anomalies dataset.  (shape (nlat, nlon, 12))
+    """
+    target = (
+        intermediate_dir
+        / "spatial_anomalies"
+        / "{obs}_{model}_{variable}_{latmin}_{latmax}_{lonmin}_{lonmax}_{train_dates[0]}_{train_dates[1]}".format(
+            **asdict(run_parameters)
+        )
+    )
+    if use_cache and zmetadata_exists(target):
+        print(f"found existing target: {target}")
+        return target
+    interpolated_obs_full_time_ds = xr.open_zarr(interpolated_obs_full_time_path)
+    obs_full_time_ds = xr.open_zarr(obs_full_time_path)
+
+    # calculate the difference between the actual obs (with finer spatial heterogeneity)
+    # and the interpolated coarse obs this will be saved and added to the
+    # spatially-interpolated coarse predictions to add the spatial heterogeneity back in.
+    spatial_anomalies = obs_full_time_ds - interpolated_obs_full_time_ds
+    seasonal_cycle_spatial_anomalies = spatial_anomalies.groupby("time.month").mean()
+
+    seasonal_cycle_spatial_anomalies.to_zarr(target, mode="w")
+
+    return target
+
+
+@task(log_stdout=True)
+def fit_and_predict(
+    experiment_train_full_time_path: UPath,
+    experiment_predict_full_time_path: UPath,
+    coarse_obs_full_time_path: UPath,
+    run_parameters: RunParameters,
+) -> UPath:
+    """Fit bcsd model on prepared CMIP data with obs at corresponding spatial scale.
+    Then predict for a set of CMIP data (likely future).
+
+    Parameters
+    ----------
+    experiment_train_full_time_path : UPath
+        UPath to experiment training dataset chunked in full time
+    experiment_predict_full_time_path : UPath
+        UPath to experiment prediction dataset chunked in full time
+    coarse_obs_full_time_path : UPath
+        UPath to coarse observation dataset chunked in full time
+    run_parameters : RunParameters
+        Prefect run parameters
+
+    Returns
+    -------
+    UPath
+        UPath to prediction results dataset.
+
+    Raises
+    ------
+    ValueError
+        ValueError checking validity of input variables.
+    """
+
+    target = intermediate_dir / "fit_and_predict" / run_parameters.run_id
+    if use_cache and zmetadata_exists(target):
+        print(f"found existing target: {target}")
+        return target
+
+    if run_parameters.variable in ABSOLUTE_VARS:
+        bcsd_model = BcsdTemperature(return_anoms=False)
+    elif run_parameters.variable in RELATIVE_VARS:
+        bcsd_model = BcsdPrecipitation(return_anoms=False)
+    else:
+        raise ValueError('variable not found in ABSOLUTE_VARS OR RELATIVE_VARS.')
+
+    pointwise_model = PointWiseDownscaler(model=bcsd_model, dim="time")
+
+    coarse_obs_full_time_ds = xr.open_zarr(coarse_obs_full_time_path)
+    experiment_train_full_time_ds = xr.open_zarr(experiment_train_full_time_path)
+    experiment_predict_full_time_ds = xr.open_zarr(experiment_predict_full_time_path)
+
+    pointwise_model.fit(
+        experiment_train_full_time_ds[run_parameters.variable],
+        coarse_obs_full_time_ds[run_parameters.variable],
+    )
+    bias_corrected_da = pointwise_model.predict(
+        experiment_predict_full_time_ds[run_parameters.variable]
+    )
+
+    bias_corrected_ds = bias_corrected_da.astype('float32').to_dataset(name=run_parameters.variable)
+    bias_corrected_ds.to_zarr(target, mode='w')
+    return target
+
+
+@task
+def postprocess_bcsd(
+    bias_corrected_fine_full_time_path: UPath,
+    spatial_anomalies_path: UPath,
+    run_parameters: RunParameters,
+) -> UPath:
+    """Downscale the bias-corrected data (fit_and_predict results) by interpolating and then
+    adding the spatial anomalies back in.
+
+    Parameters
+    ----------
+    bias_corrected_fine_full_time_path : UPath
+        UPath to output dataset from the fit_and_predict task.
+    spatial_anomalies_path : UPath
+        UPath to the output of the spatial_anomalies task.
+    run_parameters : RunParameters
+        Prefect run parameters.
+
+    Returns
+    -------
+    UPath
+        UPath to post-processed dataset.
+    """
+
+    target = results_dir / "daily" / run_parameters.run_id
+    if use_cache and zmetadata_exists(target):
+        print(f"found existing target: {target}")
+        return target
+
+    bias_corrected_fine_full_time_ds = xr.open_zarr(bias_corrected_fine_full_time_path)
+    # hint for mapblocks about which month each day corresponds to
+    bias_corrected_fine_full_time_ds = bias_corrected_fine_full_time_ds.assign_coords(
+        {'month': bias_corrected_fine_full_time_ds['time.month']}
+    )
+    spatial_anomalies_ds = xr.open_zarr(spatial_anomalies_path)
+    # make all spatial anomalies into one chunk so that map_blocks has access to every month.
+    # otherwise it will only have access to one chunk and will only grab the last chunk (december)
+    # and result in nans in all months except for december
+    spatial_anomalies_ds = spatial_anomalies_ds.chunk({'month': -1}).persist()
+    bcsd_results_ds = xr.map_blocks(
+        reconstruct_finescale,
+        bias_corrected_fine_full_time_ds,
+        args=[spatial_anomalies_ds],
+        template=bias_corrected_fine_full_time_ds,
+    )
+    bcsd_results_ds.to_zarr(target, mode='w')
+    return target
diff --git a/cmip6_downscaling/methods/common/containers.py b/cmip6_downscaling/methods/common/containers.py
new file mode 100644
index 00000000..7539dbd4
--- /dev/null
+++ b/cmip6_downscaling/methods/common/containers.py
@@ -0,0 +1,94 @@
+from dataclasses import dataclass
+
+
+@dataclass
+class BBox:
+    """Lat/lon bounding box"""
+
+    latmin: float = -90.0
+    latmax: float = 90.0
+    lonmin: float = -180.0
+    lonmax: float = 180
+
+    @property
+    def lat_slice(self) -> slice:
+        return slice(float(self.latmax), float(self.latmin))
+
+    @property
+    def lon_slice(self) -> slice:
+        return slice(float(self.lonmin), float(self.lonmax))
+
+    def __str__(self) -> str:
+        return f"BBox([{self.latmin}, {self.latmax}, {self.lonmin}, {self.lonmax}])"
+
+
+@dataclass
+class TimePeriod:
+    """Time period object"""
+
+    start: str = None
+    stop: str = None
+
+    @property
+    def time_slice(self):
+        return slice(self.start, self.stop)
+
+    def __str__(self) -> str:
+        return f"TimePeriod([{self.start}, {self.stop}])"
+
+
+@dataclass
+class CMIP6Experiment:
+
+    model: str
+    scenario: str
+    member: str
+
+
+@dataclass
+class RunParameters:
+    """Run parameters object"""
+
+    method: str
+    obs: str
+    model: str
+    member: str
+    grid_label: str
+    table_id: str
+    scenario: str
+    variable: str
+    latmin: float
+    latmax: float
+    lonmin: float
+    lonmax: float
+    train_dates: list
+    predict_dates: list
+
+    @property
+    def bbox(self) -> BBox:
+        return BBox(
+            latmin=self.latmin,
+            latmax=self.latmax,
+            lonmin=self.lonmin,
+            lonmax=self.lonmax,
+        )
+
+    @property
+    def train_period(self) -> TimePeriod:
+        if len(self.train_dates) != 2:
+            raise ValueError("expected train_dates to be a list of length 2")
+        return TimePeriod(start=self.train_dates[0], stop=self.train_dates[1])
+
+    @property
+    def predict_period(self) -> TimePeriod:
+        if len(self.predict_dates) != 2:
+            raise ValueError("expected predict_dates to be a list of length 2")
+        return TimePeriod(start=self.predict_dates[0], stop=self.predict_dates[1])
+
+    @property
+    def experiment(self) -> CMIP6Experiment:
+        return CMIP6Experiment(model=self.model, scenario=self.scenario, member=self.member)
+
+    @property
+    def run_id(self):
+        return f"{self.method}_{self.obs}_{self.model}_{self.scenario}_{self.variable}_{self.latmin}_{self.latmax}_{self.lonmin}_{self.lonmax}_{self.train_dates[0]}_{self.train_dates[1]}_{self.predict_dates[0]}_{self.predict_dates[1]}"
diff --git a/cmip6_downscaling/methods/common/tasks.py b/cmip6_downscaling/methods/common/tasks.py
new file mode 100644
index 00000000..e96872ac
--- /dev/null
+++ b/cmip6_downscaling/methods/common/tasks.py
@@ -0,0 +1,557 @@
+import os
+import warnings
+from dataclasses import asdict
+from pathlib import PosixPath
+
+import dask
+import datatree as dt
+import fsspec
+import rechunker
+import xarray as xr
+import zarr
+from carbonplan_data.metadata import get_cf_global_attrs
+from carbonplan_data.utils import set_zarr_encoding
+from ndpyramid import pyramid_regrid
+from prefect import task
+from upath import UPath
+from xarray_schema import DataArraySchema, DatasetSchema
+from xarray_schema.base import SchemaError
+
+import cmip6_downscaling
+from cmip6_downscaling import config, version
+from cmip6_downscaling.data.cmip import get_gcm
+from cmip6_downscaling.data.observations import open_era5
+from cmip6_downscaling.methods.common.utils import (
+    calc_auspicious_chunks_dict,
+    subset_dataset,
+    zmetadata_exists,
+)
+
+from .containers import RunParameters
+
+warnings.filterwarnings(
+    "ignore",
+    "(.*) filesystem path not explicitly implemented. falling back to default implementation. This filesystem may not be tested",
+    category=UserWarning,
+)
+
+PIXELS_PER_TILE = 128
+code_version = cmip6_downscaling.__version__
+scratch_dir = UPath(config.get("storage.scratch.uri"))
+intermediate_dir = UPath(config.get("storage.intermediate.uri")) / cmip6_downscaling.__version__
+results_dir = UPath(config.get("storage.results.uri")) / cmip6_downscaling.__version__
+use_cache = config.get('run_options.use_cache')
+
+
+@task
+def make_run_parameters(**kwargs) -> RunParameters:
+    """Prefect task wrapper for RunParameters"""
+    return RunParameters(**kwargs)
+
+
+@task(log_stdout=True)
+def get_obs(run_parameters: RunParameters) -> UPath:
+    """Task to return observation data subset from input parameters.
+
+    Parameters
+    ----------
+    run_parameters : RunParameters
+        RunParameter dataclass defined in common/conatiners.py. Constructed from prefect parameters.
+
+    Returns
+    -------
+    UPath
+        Path to subset observation dataset.
+    """
+
+    ds_name = (
+        "get_obs"
+        + "/"
+        + "{obs}_{variable}_{latmin}_{latmax}_{lonmin}_{lonmax}_{train_dates[0]}_{train_dates[1]}".format(
+            **asdict(run_parameters)
+        )
+    )
+    target = str(intermediate_dir) + "/" + ds_name
+
+    if use_cache and zmetadata_exists(target):
+        print(f'found existing target: {target}')
+        return target
+
+    ds = open_era5(run_parameters.variable, run_parameters.train_period)
+
+    subset = subset_dataset(
+        ds,
+        run_parameters.variable,
+        run_parameters.train_period.time_slice,
+        run_parameters.bbox,
+        chunking_schema={'time': 365, 'lat': 150, 'lon': 150},
+    )
+    del subset[run_parameters.variable].encoding['chunks']
+
+    subset.attrs.update({'title': ds_name}, **get_cf_global_attrs(version=version))
+    subset.to_zarr(target, mode='w')
+    return target
+
+
+@task(log_stdout=True)
+def get_experiment(run_parameters: RunParameters, time_subset: str) -> UPath:
+    """Prefect task that returns cmip GCM data from input run parameters.
+
+    Parameters
+    ----------
+    run_parameters : RunParameters
+        RunParameter dataclass defined in common/conatiners.py. Constructed from prefect parameters.
+    time_subset : str
+        String describing time subset request. Either 'train_period' or 'predict_period'
+
+    Returns
+    -------
+    UPath
+        UPath to experiment dataset.
+    """
+    time_period = getattr(run_parameters, time_subset)
+    ds_name = (
+        "get_experiment"
+        + "/"
+        + "{model}_{scenario}_{variable}_{latmin}_{latmax}_{lonmin}_{lonmax}_{time_period.start}_{time_period.stop}".format(
+            time_period=time_period, **asdict(run_parameters)
+        )
+    )
+    target = str(intermediate_dir) + "/" + ds_name
+
+    if use_cache and zmetadata_exists(target):
+        print(f'found existing target: {target}')
+        return target
+
+    ds = get_gcm(
+        scenario=run_parameters.scenario,
+        member_id=run_parameters.member,
+        table_id=run_parameters.table_id,
+        grid_label=run_parameters.grid_label,
+        source_id=run_parameters.model,
+        variable=run_parameters.variable,
+    )
+
+    subset = subset_dataset(
+        ds, run_parameters.variable, time_period.time_slice, run_parameters.bbox
+    )
+
+    # Note: dataset is chunked into time:365 chunks to standardize leap-year chunking.
+    subset = subset.chunk({'time': 365})
+    del subset[run_parameters.variable].encoding['chunks']
+
+    subset.attrs.update({'title': ds_name}, **get_cf_global_attrs(version=version))
+    subset.to_zarr(target, mode='w')
+    return target
+
+
+@task(log_stdout=True)
+def rechunk(
+    path: UPath,
+    pattern: str = None,
+    template: UPath = None,
+    max_mem: str = "2GB",
+) -> UPath:
+    """Use `rechunker` package to adjust chunks of dataset to a form
+    conducive for your processing.
+
+    Parameters
+    ----------
+    path : UPath
+        path to zarr store
+    pattern : str
+        The pattern of chunking you want to use. If used together with `template` it will override the template
+        to ensure that the final dataset truly follows that `full_space` or `full_time` spec. This matters when you are passing
+        a template that is either a shorter time length or a template that is a coarser grid (and thus a shorter lat/lon chunksize)
+    template : UPath
+        The path to the file you want to use as a chunking template. The utility will grab the chunk sizes and use them as the chunk
+        target to feed to rechunker.
+    max_mem : str
+        The memory available for rechunking steps. Must look like "2GB". Optional, default is 2GB.
+
+    Returns
+    -------
+    target : UPath
+        Path to rechunked dataset
+    """
+    # print('rechunking dataset at {}'.format(path))
+    # if both defined then you'll take the spatial part of template and override one dimension with the specified pattern
+    if template is not None:
+        pattern_string = 'matched'
+        if pattern is not None:
+            pattern_string += '_' + pattern
+    # if only pattern specified then use that pattern
+    elif pattern is not None:
+        pattern_string = pattern
+    target = intermediate_dir / ("rechunk_" + pattern_string) / (str(path).split("/")[-1])
+    path_tmp = scratch_dir / ("rechunk_" + pattern_string) / (str(path).split("/")[-1])
+    print('writing rechunked dataset to {}'.format(target))
+    target_store = fsspec.get_mapper(str(target))
+    temp_store = fsspec.get_mapper(str(path_tmp))
+
+    if use_cache and zmetadata_exists(target):
+        print(f'found existing target: {target}')
+        # if we wanted to check that it was chunked correctly we could put this down below where
+        # the target_schema is validated. but that requires us going through the development
+        # of the schema would just hurt performance likely unnecessarily.
+        # nevertheless, as future note: if we encounter chunk issues i suggest putting a schema check here
+        return target
+    # if a cached target isn't found we'll go through the rechunking step
+    # open the zarr group
+    target_store.clear()
+    temp_store.clear()
+    group = zarr.open_consolidated(path)
+    # open the dataset to access the coordinates
+    ds = xr.open_zarr(path)
+    example_var = list(ds.data_vars)[0]
+    # if you have defined a template then use the chunks of that template
+    # to form the desired chunk definition
+    if template is not None:
+        template_ds = xr.open_zarr(template)
+        # define the chunk definition
+        chunk_def = {
+            'time': min(template_ds.chunks['time'][0], len(ds.time)),
+            'lat': min(template_ds.chunks['lat'][0], len(ds.lat)),
+            'lon': min(template_ds.chunks['lon'][0], len(ds.lon)),
+        }
+        # if you have also defined a pattern then override the dimension you've specified there
+        if pattern is not None:
+            # the chunking pattern will return the dimensions that you'll chunk along
+            # so `full_time` will return `('lat', 'lon')`
+            chunk_dims = config.get(f"chunk_dims.{pattern}")
+            for dim in chunk_def.keys():
+                if dim not in chunk_dims:
+                    print('correcting dim')
+                    # override the chunksize of those unchunked dimensions to be the complete length (like passing chunksize=-1
+                    chunk_def[dim] = len(ds[dim])
+    # if you don't have a target template then you'll just use the `full_time` or `full_space` approach
+    elif pattern is not None:
+        chunk_dims = config.get(f"chunk_dims.{pattern}")
+        chunk_def = calc_auspicious_chunks_dict(ds[example_var], chunk_dims=chunk_dims)
+    else:
+        raise AttributeError('must either define chunking pattern or template')
+    # Note:
+    # for rechunker v 0.3.3:
+    # initialize the chunks_dict that you'll pass in, filling the coordinates with
+    # `None`` because you don't want to rechunk the coordinate arrays. this works with
+    # for rechunker v 0.4.2:
+    # initialize chunks_dict using the `chunk_def`` above
+    chunks_dict = {
+        'time': (chunk_def['time'],),
+        'lon': (chunk_def['lon'],),
+        'lat': (chunk_def['lat'],),
+    }
+    for var in ds.data_vars:
+        chunks_dict[var] = chunk_def
+    # now that you have your chunks_dict you can check that the dataset at `path`
+    # you're passing in doesn't already match that schema. because if so, we don't
+    # need to bother with rechunking and we'll skip it!
+    schema_dict = {}
+    for var in ds.data_vars:
+        schema_dict[var] = DataArraySchema(chunks=chunk_def)
+    target_schema = DatasetSchema(schema_dict)
+    try:
+        # check to see if the initial dataset already matches the schema, in which case just
+        # return the initial path and work with that
+        target_schema.validate(ds)
+        return path
+    except SchemaError:
+        pass
+    rechunk_plan = rechunker.rechunk(
+        source=group,
+        target_chunks=chunks_dict,
+        max_mem=max_mem,
+        target_store=target_store,
+        temp_store=temp_store,
+    )
+
+    rechunk_plan.execute()
+
+    # consolidate_metadata here since when it comes out of rechunker it isn't consolidated.
+    zarr.consolidate_metadata(target_store)
+    temp_store.clear()
+
+    return target
+
+
+@task
+def monthly_summary(ds_path: UPath, run_parameters: RunParameters) -> UPath:
+    """Prefect task to create monthly resampled data. Takes mean of `tasmax` and `tasmin` and sum of `pr`.
+
+    Parameters
+    ----------
+    ds_path : UPath
+        UPath to input zarr store at daily timestep
+    run_parameters : RunParameters
+        prefect run parameters
+
+    Returns
+    -------
+    UPath
+        Path to resampled dataset.
+    """
+
+    ds_name = "monthly_summary" + "/" + str(run_parameters.run_id)
+    target = str(results_dir) + "/" + ds_name
+
+    if use_cache and zmetadata_exists(target):
+        print(f'found existing target: {target}')
+        return target
+
+    ds = xr.open_zarr(ds_path)
+
+    if run_parameters.variable in ['tasmax', 'tasmin']:
+        out_ds = ds.resample(time='1MS').mean(dim='time')
+    elif run_parameters.variable in ['pr']:
+        out_ds = ds.resample(time='1MS').sum(dim='time')
+    else:
+        print(f'{run_parameters.variable} not implemented')
+
+    out_ds.attrs.update({'title': ds_name}, **get_cf_global_attrs(version=version))
+
+    out_ds.to_zarr(target, mode='w')
+
+    return target
+
+
+@task
+def annual_summary(ds_path: UPath, run_parameters: RunParameters) -> UPath:
+    """Prefect task to create yearly resampled data. Takes mean of `tasmax` and `tasmin` and sum of `pr`.
+
+    Parameters
+    ----------
+    ds_path : UPath
+        UPath to input zarr store at daily timestep
+    run_parameters : RunParameters
+        prefect run parameters
+
+    Returns
+    -------
+    UPath
+        Path to resampled dataset.
+    """
+
+    ds_name = "annual_summary" + "/" + str(run_parameters.run_id)
+    target = str(results_dir) + "/" + ds_name
+
+    if use_cache and zmetadata_exists(target):
+        print(f'found existing target: {target}')
+        return target
+
+    ds = xr.open_zarr(ds_path)
+
+    if run_parameters.variable in ['tasmax', 'tasmin']:
+        out_ds = ds.resample(time='YS').mean(dim='time')
+    elif run_parameters.variable in ['pr']:
+        out_ds = ds.resample(time='YS').sum(dim='time')
+    else:
+        print(f'{run_parameters.variable} not implemented')
+
+    out_ds.attrs.update({'title': ds_name}, **get_cf_global_attrs(version=version))
+    out_ds.to_zarr(target, mode='w')
+
+    return target
+
+
+@task(tags=['dask-resource:TASKSLOTS=1'], log_stdout=True)
+def regrid(source_path: UPath, target_grid_path: UPath) -> UPath:
+    """Task to regrid a dataset to target grid.
+
+    Parameters
+    ----------
+    source_path : UPath
+        Path to dataset that will be regridded
+    target_grid_path : UPath
+        Path to template grid dataset
+
+    Returns
+    -------
+    UPath
+        Path to regridded output dataset.
+    """
+
+    import xesmf as xe
+
+    ds_name = (
+        "regrid"
+        + "/"
+        + "source_path"
+        + "/"
+        + (str(source_path).split("/")[-1])
+        + "/"
+        + "target_path"
+        + "/"
+        + (str(target_grid_path).split("/")[-1])
+    )
+    target = str(intermediate_dir) + "/" + ds_name
+
+    if use_cache and zmetadata_exists(target):
+        print(f'found existing target: {target}')
+        return target
+    source_ds = xr.open_zarr(source_path)
+    target_grid_ds = xr.open_zarr(target_grid_path)
+
+    regridder = xe.Regridder(source_ds, target_grid_ds, "bilinear", extrap_method="nearest_s2d")
+    regridded_ds = regridder(source_ds)
+    regridded_ds.attrs.update({'title': ds_name}, **get_cf_global_attrs(version=version))
+    regridded_ds.to_zarr(target, mode='w')
+
+    return target
+
+
+def _load_coords(ds: xr.Dataset) -> xr.Dataset:
+    '''Helper function to explicitly load all dataset coordinates'''
+    for var, da in ds.coords.items():
+        ds[var] = da.load()
+    return ds
+
+
+def _pyramid_postprocess(
+    dt: dt.DataTree, levels: int, other_chunks: dict = None, ds_name: str = None
+) -> dt.DataTree:
+    '''Postprocess data pyramid
+
+    Adds multiscales metadata and sets Zarr encoding
+
+    Parameters
+    ----------
+    dt : dt.DataTree
+        Input data pyramid
+    levels : int
+        Number of levels in pyramid
+    other_chunks : dict
+        Chunks for non-spatial dims
+
+    Returns
+    -------
+    dt.DataTree
+        Updated data pyramid with metadata / encoding set
+    '''
+    chunks = {"x": PIXELS_PER_TILE, "y": PIXELS_PER_TILE}
+    if other_chunks is not None:
+        chunks.update(other_chunks)
+
+    for level in range(levels):
+        slevel = str(level)
+        dt.ds.attrs['multiscales'][0]['datasets'][level]['pixels_per_tile'] = PIXELS_PER_TILE
+
+        # set dataset chunks
+        dt[slevel].ds = dt[slevel].ds.chunk(chunks)
+        if 'date_str' in dt[slevel].ds:
+            dt[slevel].ds['date_str'] = dt[slevel].ds['date_str'].chunk(-1)
+
+        # set dataset encoding
+        dt[slevel].ds = set_zarr_encoding(
+            dt[slevel].ds, codec_config={"id": "zlib", "level": 1}, float_dtype="float32"
+        )
+        for var in ['time', 'time_bnds']:
+            if var in dt[slevel].ds:
+                dt[slevel].ds[var].encoding['dtype'] = 'int32'
+
+    # set global metadata
+    dt.ds.attrs.update({'title': ds_name}, **get_cf_global_attrs(version=version))
+    return dt
+
+
+@task(log_stdout=True, tags=['dask-resource:TASKSLOTS=1'])
+def pyramid(
+    ds_path: UPath, run_parameters: RunParameters, levels: int = 2, other_chunks: dict = None
+) -> UPath:
+    '''Task to create a data pyramid from an xarray Dataset
+
+    Parameters
+    ----------
+    ds_path : UPath
+        Path to input dataset
+    levels : int, optional
+        Number of levels in pyramid, by default 2
+    uri : str, optional
+        Path to write output data pyamid to, by default None
+    other_chunks : dict
+        Chunks for non-spatial dims
+
+    Returns
+    -------
+    target : UPath
+    '''
+
+    ds_name = "pyarmid" + ds_path.path.replace('/', '_')
+    target = str(results_dir) + "/" + ds_name
+
+    if use_cache and zmetadata_exists(target):
+        print(f'found existing target: {target}')
+        return target
+
+    ds = xr.open_zarr(ds_path).pipe(_load_coords)
+
+    ds.coords['date_str'] = ds['time'].dt.strftime('%Y-%m-%d').astype('S10')
+    # note: this worked when 8 processors and 4 cores
+    with dask.config.set(scheduler='threads'):
+        # create pyramid
+        dta = pyramid_regrid(ds, target_pyramid=None, levels=levels)
+
+        dta = _pyramid_postprocess(dta, levels, other_chunks=other_chunks, ds_name=ds_name)
+
+    """Prefect task to run the analyses on results from a downscaling run.
+
+    Parameters
+    ----------
+    ds_path : UPath
+    run_parameters : RunParameters
+        Downscaling run parameter container
+
+    Returns
+    -------
+    PosixPath
+        The local location of an executed notebook path.
+    """
+
+    import papermill
+    from azure.storage.blob import BlobServiceClient, ContentSettings
+
+    from cmip6_downscaling.analysis import metrics
+
+    root = PosixPath(metrics.__file__)
+    template_path = root.parent / 'analyses_template.ipynb'
+    executed_notebook_path = root.parent / f'analyses_{run_parameters.run_id}.ipynb'
+    executed_html_path = root.parent / f'analyses_{run_parameters.run_id}.html'
+
+    parameters = asdict(run_parameters)
+    parameters['run_id'] = run_parameters.run_id
+    # TODO: figure out how to unpack these fields in the notebook
+    # asdict will return lists for train_dates and predict_dates
+    # parameters['train_period_start'] = train_period.start
+    # parameters['train_period_end'] = train_period.stop
+    # parameters['predict_period_start'] = predict_period.start
+    # parameters['predict_period_end'] = predict_period.stop
+
+    # execute notebook with papermill
+    papermill.execute_notebook(template_path, executed_notebook_path, parameters=parameters)
+
+    # convert from ipynb to html
+    # TODO: move this to stand alone function
+    # Q: can we control the output path name?
+    os.system(f"jupyter nbconvert {executed_notebook_path} --to html")
+
+    # TODO: move to stand alone function
+    connection_string = os.getenv('AZURE_STORAGE_CONNECTION_STRING', None)
+    if connection_string is not None:
+        # if you have a connection_string, copy the html to azure, if not just return
+        # because it is already in your local machine
+        blob_service_client = BlobServiceClient.from_connection_string(connection_string)
+        # TODO: fix b/c the run_id has slashes now!!!
+        blob_name = config.get('storage.web_results.blob') / parameters.run_id / 'analyses.html'
+        blob_client = blob_service_client.get_blob_client(container='$web', blob=blob_name)
+        # clean up before writing
+        try:
+            blob_client.delete_blob()
+        except:  # TODO: raise specific error
+            pass
+
+        #  need to specify html content type so that it will render and not download
+        with open(executed_html_path, "rb") as data:
+            blob_client.upload_blob(
+                data, content_settings=ContentSettings(content_type='text/html')
+            )
+
+    return executed_notebook_path
diff --git a/cmip6_downscaling/methods/common/utils.py b/cmip6_downscaling/methods/common/utils.py
new file mode 100644
index 00000000..27bb90fa
--- /dev/null
+++ b/cmip6_downscaling/methods/common/utils.py
@@ -0,0 +1,128 @@
+import pathlib
+import re
+from typing import Tuple
+
+import numpy as np
+import xarray as xr
+from upath import UPath
+from xarray_schema import DataArraySchema
+from xarray_schema.base import SchemaError
+
+from .containers import BBox
+
+
+def zmetadata_exists(path: UPath):
+    '''temporary workaround until path.exists() works'''
+
+    if isinstance(path, pathlib.PosixPath):
+        return (path / '.zmetadata').exists()
+    elif isinstance(path, UPath):
+        return path.fs.exists(str(path / '.zmetadata'))
+    else:
+        return (UPath(path) / '.zmetadata').exists()
+
+
+def subset_dataset(
+    ds: xr.Dataset,
+    variable: str,
+    time_period: slice,
+    bbox: BBox,
+    chunking_schema: dict = None,
+) -> xr.Dataset:
+    """Uses Xarray slicing to spatially subset a dataset based on input params.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+         Input Xarray dataset
+    time_period : slice
+        Start and end year slice. Ex: slice('2020','2020')
+    bbox : dataclass
+        dataclass containing the latmin,latmax,lonmin,lonmax. Class can be found in utils.
+    chunking_schema : str, optional
+        Desired chunking schema. ex: {'time': 365, 'lat': 150, 'lon': 150}
+
+    Returns
+    -------
+    Xarray Dataset
+        Spatially subsetted Xarray dataset.
+    """
+
+    subset_ds = ds.sel(
+        time=time_period,
+        lon=bbox.lon_slice,
+        lat=bbox.lat_slice,
+    )
+    if chunking_schema is not None:
+        target_schema = DataArraySchema(chunks=chunking_schema)
+        try:
+            target_schema.validate(subset_ds[variable])
+        except SchemaError:
+            subset_ds = subset_ds.chunk(chunking_schema)
+
+    return subset_ds
+
+
+def calc_auspicious_chunks_dict(
+    da: xr.DataArray,
+    chunk_dims: Tuple = ("lat", "lon"),
+) -> dict:
+    """Figure out a chunk size that, given the size of the dataset, the dimension(s) you want to chunk on
+    and the data type, will fit under the target_size. Currently only works for 100mb which
+    is the recommended chunk size for dask processing.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+        Dataset or data array you're wanting to chunk
+    chunk_dims : tuple, optional
+        Dimension(s) you want to chunk along, by default ('lat', 'lon')
+
+    Returns
+    -------
+    chunks_dict : dict
+        Dictionary of chunk sizes with the dimensions not listed in `chunk_dims` being the
+        length of that dimension (avoiding the shorthand -1 in order to play nice
+        with rechunker)
+    """
+    if not isinstance(chunk_dims, tuple):
+        raise TypeError(
+            "Your chunk_dims likely includes one string but needs a comma after it! to be a tuple!"
+        )
+    # setting target_size_bytes to the 100mb chunksize recommended by dask. could modify in future.
+    target_size_bytes = 100e6
+    dim_sizes = dict(zip(da.dims, da.shape))
+
+    # initialize chunks_dict
+    chunks_dict = {}
+
+    # dims not in chunk_dims should be one chunk (length -1), since these ones are going to
+    # be contiguous while the dims in chunk_dims will be chunked
+    for dim in dim_sizes.keys():
+        if dim not in chunk_dims:
+            # we'll only add the unchunked dimensions to chunks_dict right now
+            # rechunker doesn't like the the shorthand of -1 meaning the full length
+            # so we'll always just give it the full length of the dimension
+            chunks_dict[dim] = dim_sizes[dim]
+    # calculate the bytesize given the dtype bitsize and divide by 8
+    data_bytesize = int(re.findall(r"\d+", str(da.dtype))[0]) / 8
+    # calculate the size of the smallest minimum chunk based upon dtype and the
+    # length of the unchunked dim(s). chunks_dict currently only has unchunked dims right now
+    smallest_size_one_chunk = data_bytesize * np.prod(
+        [dim_sizes[dim] for dim in chunks_dict.keys()]
+    )
+
+    # the dims in chunk_dims should be of an array size (as in number of elements in the array)
+    # that creates ~100 mb. `perfect_chunk` is the how many of the smallest_size_chunks you can
+    # handle at once while still staying below the `target_size_bytes`
+    perfect_chunk = target_size_bytes / smallest_size_one_chunk
+
+    # then make reasonable chunk size by rounding up (avoids corner case of it rounding down to 0...)
+    # but if the array is oblong it might get big (? is logic right there- might it get small??)
+    perfect_chunk_length = int(np.ceil(perfect_chunk ** (1 / len(chunk_dims))))
+    for dim in chunk_dims:
+        # check that the rounding up as part of the `perfect_chunk_length` calculation
+        # didn't make the chunk sizes bigger than the array itself, and if so
+        # clip it to that size
+        chunks_dict[dim] = min(perfect_chunk_length, dim_sizes[dim])
+    return chunks_dict
diff --git a/cmip6_downscaling/runtimes.py b/cmip6_downscaling/runtimes.py
index 8e10a144..44844ffd 100644
--- a/cmip6_downscaling/runtimes.py
+++ b/cmip6_downscaling/runtimes.py
@@ -1,6 +1,7 @@
 import os
 from abc import abstractmethod
 
+import dask
 from dask_kubernetes import KubeCluster, make_pod_spec
 from prefect.executors import DaskExecutor, Executor, LocalDaskExecutor, LocalExecutor
 from prefect.run_configs import KubernetesRun, LocalRun, RunConfig
@@ -17,7 +18,14 @@
 
 
 class BaseRuntime:
-    """Base configuration class that defines abstract methods (storage, run_config and executor) for subclasses."""
+    """
+    Base runtime configuration class
+
+    Subclasses must define methods (storage, run_config and executor)
+    """
+
+    def __init__(self):
+        pass
 
     @property
     @abstractmethod
@@ -34,71 +42,22 @@ def run_config(self) -> RunConfig:  # pragma: no cover
     def executor(self) -> Executor:  # pragma: no cover
         pass
 
-
-class CloudRuntime(BaseRuntime):
-    def __init__(
-        self,
-        runtime_name="cloud",
-        storage_options=None,
-        agent=None,
-        extra_pip_packages=None,
-        kubernetes_cpu=None,
-        kubernetes_memory=None,
-        image=None,
-        pod_memory_limit=None,
-        pod_memory_request=None,
-        pod_threads_per_worker=None,
-        pod_cpu_limit=None,
-        pod_cpu_request=None,
-        deploy_mode=None,
-        adapt_min=None,
-        adapt_max=None,
-        dask_distributed_worker_resources_taskslots=None,
-    ):
-        self._runtime_name = runtime_name
-        self._storage_options = storage_options is not None or config.get(
-            "runtime.cloud.storage_options.container"
-        )
-
-        self._agent = agent is not None or config.get("runtime.cloud.agent")
-        self._extra_pip_packages = extra_pip_packages is not None or config.get(
-            "runtime.cloud.extra_pip_packages"
-        )
-        self._kubernetes_cpu = kubernetes_cpu is not None or config.get(
-            "runtime.cloud.kubernetes_cpu"
-        )
-        self._kubernetes_memory = kubernetes_memory is not None or config.get(
-            "runtime.cloud.kubernetes_memory"
-        )
-        self._image = image is not None or config.get("runtime.cloud.image")
-        self._pod_memory_limit = pod_memory_limit is not None or config.get(
-            "runtime.cloud.pod_memory_limit"
-        )
-        self._pod_memory_request = pod_memory_request is not None or config.get(
-            "runtime.cloud.pod_memory_request"
-        )
-        self._pod_threads_per_worker = pod_threads_per_worker is not None or config.get(
-            "runtime.cloud.pod_threads_per_worker"
-        )
-        self._pod_cpu_limit = pod_cpu_limit is not None or config.get("runtime.cloud.pod_cpu_limit")
-        self._pod_cpu_request = pod_cpu_request is not None or config.get(
-            "runtime.cloud.pod_cpu_request"
-        )
-        self._deploy_mode = deploy_mode is not None or config.get("runtime.cloud.deploy_mode")
-        self._adapt_min = adapt_min is not None or config.get("runtime.cloud.adapt_min")
-        self._adapt_max = adapt_max is not None or config.get("runtime.cloud.adapt_max")
-        self._dask_distributed_worker_resources_taskslots = (
-            dask_distributed_worker_resources_taskslots is not None
-            or config.get("runtime.cloud.dask_distributed_worker_resources_taskslots")
+    def __repr__(self):
+        return (
+            f"{type(self)}\n"
+            f"  Storage    : {type(self.storage)}\n"
+            f"  Run Config : {type(self.run_config)}\n"
+            f"  Executor   : {type(self.executor)}\n"
         )
 
-    def __repr__(self):
-        return """CloudRuntime configuration for running on prefect-cloud.  storage is `Azure("prefect")`, run_config is `KubernetesRun() and executor is ` DaskExecutor(KubeCluster())`"""
 
+class CloudRuntime(BaseRuntime):
     def _generate_env(self):
         env = {
-            "EXTRA_PIP_PACKAGES": self._extra_pip_packages,
-            "DASK_DISTRIBUTED__WORKER__RESOURCES__TASKSLOTS": self._dask_distributed_worker_resources_taskslots,
+            "EXTRA_PIP_PACKAGES": config.get("runtime.cloud.extra_pip_packages"),
+            "DASK_DISTRIBUTED__WORKER__RESOURCES__TASKSLOTS": config.get(
+                "runtime.cloud.dask_distributed_worker_resources_taskslots"
+            ),
             "PREFECT__FLOWS__CHECKPOINTING": "true",
             **_threadsafe_env_vars,
         }
@@ -110,53 +69,51 @@ def _generate_env(self):
 
     @property
     def storage(self) -> Storage:
-        return Azure(self._storage_options)
+        return Azure(config.get("runtime.cloud.storage_options.container"))
 
     @property
     def run_config(self) -> RunConfig:
-        kube_run_config = KubernetesRun(
-            cpu_request=self._kubernetes_cpu,
-            memory_request=self._kubernetes_memory,
-            image=self._image,
-            labels=[self._agent],
+        return KubernetesRun(
+            cpu_request=config.get("runtime.cloud.kubernetes_cpu"),
+            memory_request=config.get("runtime.cloud.kubernetes_memory"),
+            image=config.get("runtime.cloud.image"),
+            labels=[config.get("runtime.cloud.agent")],
             env=self._generate_env(),
         )
-        return kube_run_config
 
     @property
     def executor(self) -> Executor:
 
         pod_spec = make_pod_spec(
-            image=self._image,
-            memory_limit=self._pod_memory_limit,
-            memory_request=self._pod_memory_request,
-            threads_per_worker=self._pod_threads_per_worker,
-            cpu_limit=self._pod_cpu_limit,
-            cpu_request=self._pod_cpu_request,
+            image=config.get("runtime.cloud.image"),
+            memory_limit=config.get("runtime.cloud.pod_memory_limit"),
+            memory_request=config.get("runtime.cloud.pod_memory_request"),
+            threads_per_worker=config.get("runtime.cloud.pod_threads_per_worker"),
+            cpu_limit=config.get("runtime.cloud.pod_cpu_limit"),
+            cpu_request=config.get("runtime.cloud.pod_cpu_request"),
             env=self._generate_env(),
         )
+
+        # TODO: this can done in the constructor now
         pod_spec.spec.containers[0].args.extend(["--resources", "TASKSLOTS=1"])
 
         executor = DaskExecutor(
-            cluster_class=lambda: KubeCluster(pod_spec, deploy_mode=self._deploy_mode),
-            adapt_kwargs={"minimum": self._adapt_min, "maximum": self._adapt_max},
+            cluster_class=lambda: KubeCluster(
+                pod_spec, deploy_mode=config.get("runtime.cloud.deploy_mode")
+            ),
+            adapt_kwargs={
+                "minimum": config.get("runtime.cloud.adapt_min"),
+                "maximum": config.get("runtime.cloud.adapt_max"),
+            },
         )
 
         return executor
 
 
 class LocalRuntime(BaseRuntime):
-    def __init__(self, storage_options: dict = None):
-        self._storage_options = storage_options is not None or config.get(
-            "runtime.local.storage_options"
-        )
-
-    def __repr__(self):
-        return "LocalRuntime configuration is for running on local machines. storage is `Local()`, run_config is `LocalRun() and executor is `LocalExecutor()` "
-
     @property
     def storage(self) -> Storage:
-        return Local(**self._storage_options)
+        return Local(**config.get("runtime.local.storage_options"))
 
     @property
     def run_config(self) -> RunConfig:
@@ -171,39 +128,22 @@ def _generate_env(self):
 
 
 class CIRuntime(LocalRuntime):
-    def __init__(self, storage_options: dict = None):
-        self._storage_options = storage_options is not None or config.get(
-            "runtime.test.storage_options"
-        )
+    @property
+    def storage(self) -> Storage:
+        return Local(**config.get("runtime.test.storage_options"))
 
     @property
     def executor(self) -> Executor:
         return LocalExecutor()
 
-    def __repr__(self):
-        return "CIRuntime configuration is for running on CI machines. storage is `Local()`, run_config is `LocalRun() and executor is `LocalExecutor()` "
-
 
 class PangeoRuntime(LocalRuntime):
-    def __init__(
-        self, storage_options: dict = None, n_workers: int = None, threads_per_worker: int = None
-    ):
-        self._storage_options = storage_options is not None or config.get(
-            "runtime.pangeo.storage_options"
-        )
-        self._n_workers = (
-            n_workers if n_workers is not None else config.get("runtime.pangeo.n_workers")
-        )
-        self._threads_per_worker = (
-            n_workers if n_workers is not None else config.get("runtime.pangeo.threads_per_worker")
-        )
-
-    def __repr__(self):
-        return "PangeoRuntime configuration is for running on jupyter-hubs. storage is `Local()`, run_config is `LocalRun() and executor is `LocalExecutor()` "
+    def __init__(self):
+        dask.config.set({'temporary_directory': '/tmp/dask'})
 
     @property
     def storage(self) -> Storage:
-        return Local(**self._storage_options)
+        return Local(**config.get("runtime.pangeo.storage_options"))
 
     @property
     def run_config(self) -> RunConfig:
@@ -214,37 +154,31 @@ def executor(self) -> Executor:
         return DaskExecutor(
             cluster_kwargs={
                 'resources': {'TASKSLOTS': 1},
-                'n_workers': self._n_workers,
-                'threads_per_worker': self._threads_per_worker,
+                'n_workers': config.get("runtime.pangeo.n_workers"),
+                'threads_per_worker': config.get("runtime.pangeo.threads_per_worker"),
             }
         )
 
-    # def executor(self) -> Executor:
-    #     return LocalExecutor()
     def _generate_env(self):
         return _threadsafe_env_vars
 
 
-def get_runtime(**kwargs):
+def get_runtime():
     if config.get("run_options.runtime") == "ci":
-        runtime = CIRuntime(**kwargs)
+        runtime = CIRuntime()
     elif config.get("run_options.runtime") == "local":
-        runtime = LocalRuntime(**kwargs)
+        runtime = LocalRuntime()
     elif config.get("run_options.runtime") == "cloud":
-        runtime = CloudRuntime(**kwargs)
+        runtime = CloudRuntime()
     elif config.get("run_options.runtime") == "pangeo":
-        runtime = PangeoRuntime(**kwargs)
+        runtime = PangeoRuntime()
     elif os.environ.get("CI") == "true":
-        runtime = CIRuntime(**kwargs)
+        runtime = CIRuntime()
         print("CIRuntime selected from os.environ")
     elif os.environ.get("PREFECT__BACKEND") == "cloud":
-        runtime = CloudRuntime(**kwargs)
-        print("PrefectCloudRuntime selected from os.environ")
+        runtime = CloudRuntime()
     elif "JUPYTER_IMAGE" in os.environ:
-        runtime = PangeoRuntime(**kwargs)
-        print("PangeoRuntime selected from os.environ")
+        runtime = PangeoRuntime()
     else:
-        raise ValueError(
-            "Name not in ['ci', 'local', 'cloud', 'pangeo'] and environment variable not found for: [CI, PREFECT__BACKEND, PANGEO__BACKEND or TEST]"
-        )
+        raise ValueError('unable to determine runtime')
     return runtime
diff --git a/cmip6_downscaling/tasks/common_tasks.py b/cmip6_downscaling/tasks/common_tasks.py
deleted file mode 100644
index a2a7cf65..00000000
--- a/cmip6_downscaling/tasks/common_tasks.py
+++ /dev/null
@@ -1,473 +0,0 @@
-from dataclasses import dataclass
-from typing import Any, Callable, Dict, List, Optional, Tuple, Union
-
-import numpy as np
-import xarray as xr
-import xclim
-from prefect import task
-from skdownscale.pointwise_models.utils import default_none_kwargs
-from xarray.core.types import T_Xarray
-from xpersist import CacheStore
-from xpersist.prefect.result import XpersistResult
-
-from cmip6_downscaling import config
-from cmip6_downscaling.data.cmip import get_gcm, get_gcm_grid_spec
-from cmip6_downscaling.data.observations import get_obs
-from cmip6_downscaling.methods.bcsd import get_coarse_obs, return_obs
-from cmip6_downscaling.methods.bias_correction import (
-    bias_correct_gcm_by_method,
-    bias_correct_obs_by_method,
-)
-from cmip6_downscaling.workflows.paths import (
-    build_gcm_identifier,
-    build_obs_identifier,
-    make_annual_pyramid_path,
-    make_bias_corrected_gcm_path,
-    make_bias_corrected_obs_path,
-    make_daily_pyramid_path,
-    make_interpolated_gcm_path,
-    make_interpolated_obs_path,
-    make_monthly_pyramid_path,
-)
-from cmip6_downscaling.workflows.utils import BBox, rechunk_zarr_array_with_caching, regrid_ds
-
-get_obs_task = task(get_obs)
-get_gcm_task = task(get_gcm)
-
-
-@task
-def build_bbox(latmin: str, latmax: str, lonmin: str, lonmax: str) -> dataclass:
-    """Build bounding box out of lat/lon inputs using BBox data class defined in /utils.py
-
-    Args:
-    Paramters
-    ---------
-    latmin : float
-         Minimum latitude
-    latmax : float
-         Maximum latitude
-    lonmin : float
-         Minimum longitude
-    lonmax : float
-        Maximum longitude
-
-    Returns
-    -------
-    BBox
-    """
-
-    return BBox(
-        latmin=float(latmin), latmax=float(latmax), lonmin=float(lonmin), lonmax=float(lonmax)
-    )
-
-
-@task
-def build_time_period_slices(time_period: list) -> slice:
-    """Return slice from list containing two time strings
-
-    Parameters
-    ----------
-    time_period : list
-        Input time period list. Ex. ['1990','1991']
-
-    Returns
-    -------
-    slice
-    """
-    return slice(*time_period)
-
-
-@task
-def to_standard_calendar(obj: T_Xarray) -> T_Xarray:
-    """Convert a Dataset's calendar to the "standard calendar"
-
-    When necessary, "missing" time points are filled in using linear interpolation.
-
-    Valid input dataset calendars include: `noleap`, `365_day`, `366_day`, and `all_leap`.
-
-    Parameters
-    ----------
-    obj : xr.Dataset or xr.DataArray
-        Xarray object with a `CFTimeIndex`.
-
-    Returns
-    -------
-    obj_new : xr.Dataset or xr.DataArray
-        Xarray object with standard calendar.
-
-    Raises
-    ------
-    ValueError
-        If an invalid calendar is supplied.
-    """
-
-    orig_calendar = getattr(obj.indexes["time"], "calendar", "standard")
-    if orig_calendar == "standard":
-        return obj
-    if orig_calendar == "360_day":
-        raise ValueError("360_day calendar is not supported")
-
-    # reindex / interpolate
-    obj_new = xclim.core.calendar.convert_calendar(obj, "standard", missing=np.nan).interpolate_na(
-        dim="time", method="linear"
-    )
-
-    # reset encoding
-    obj_new["time"].encoding["calendar"] = "standard"
-
-    return obj_new
-
-
-@task
-def rechunker_task(
-    zarr_array: xr.Dataset,
-    chunking_approach: Optional[str] = None,
-    template_chunk_array: Optional[xr.Dataset] = None,
-    naming_func: Optional[Callable] = None,
-    **kwargs,
-):
-    """
-    Task to rechunk a dataset
-
-    Parameters
-    ----------
-    zarr_array : zarr or xarray dataset
-        Dataset you want to rechunk.
-    chunking_approach : str, optional
-        Has to be one of `full_space` or `full_time`. If `full_space`, the data will be rechunked such that the space dimensions are contiguous (i.e. each chunk
-        will contain full maps). If `full_time`, the data will be rechunked such that the time dimension is contiguous (i.e. each chunk will contain full time
-        series). Either the chunking approach or the template chunk array must be provided.
-    template_chunk_array: zarr or xarray dataset, optional
-        A template dataset with the desired chunksizes. Either the chunking approach or the template chunk array must be provided.
-    naming_func: callable, optional
-        A function that returns a string that represents the output caching location that the rechunk task should save to.
-        The input arguments of this naming func must be provided as kwargs to this method
-
-    Returns
-    -------
-    rechunked_ds : xr.Dataset
-        Rechunked dataset
-    """
-    if naming_func is not None:
-        output_path = naming_func(chunking_approach=chunking_approach, **kwargs)
-    else:
-        output_path = None
-
-    rechunked = rechunk_zarr_array_with_caching(
-        zarr_array=zarr_array,
-        chunking_approach=chunking_approach,
-        template_chunk_array=template_chunk_array,
-        output_path=output_path,
-    )
-
-    return rechunked
-
-
-@task(log_stdout=True)
-def path_builder_task(
-    obs: str,
-    gcm: str,
-    scenario: str,
-    variable: str,
-    train_period: slice,
-    predict_period: slice,
-    bbox: BBox,
-) -> Tuple[str, str, str, str, str, str]:
-    """
-    Take in input parameters and make string patterns that identifies the obs dataset, gcm dataset, and the gcm grid. These
-    strings will then be used to identify cached files.
-    Parameters
-    ----------
-    obs: str
-        Name of obs dataset
-    gcm: str
-        Name of gcm model
-    scenario: str
-        Name of future emission scenario
-    variable: str
-        Name of the variable used in obs and gcm dataset (including features and label)
-    train_period: slice
-        Start and end year slice of training/historical period. Ex: slice('1990','1990')
-    predict_period: slice
-        Start and end year slice of predict period. Ex: slice('2020','2020')
-    bbox: dataclass
-        dataclass containing the latmin,latmax,lonmin,lonmax. Class can be found in utils.
-
-    Returns
-    -------
-    gcm_grid_spec: str
-        A string of parameters defining the grid of GCM, including number of lat/lon points, interval between points, lower left corner, etc.
-    obs_identifier: str
-        A string of parameters defining the obs dataset used, including variables, start/end year, etc
-    gcm_identifier : str
-        A string of parameters defining the GCM dataset used, including variables, start/end year for historical and future periods, etc
-    pyramid_path_daily : str
-        A string of parameters used to build the daily pyramid path.
-    pyramid_path_monthly : str
-        A string of parameters used to build the monthly pyramid path.
-    pyramid_path_annual : str
-        A string of parameters used to build the annual pyramid path.
-    """
-    gcm_grid_spec = get_gcm_grid_spec(gcm_name=gcm)
-    obs_identifier = build_obs_identifier(
-        obs=obs,
-        variable=variable,
-        train_period=train_period,
-        bbox=bbox,
-    )
-    gcm_identifier = build_gcm_identifier(
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-    )
-    pyramid_path_daily = make_daily_pyramid_path(gcm_identifier)
-    pyramid_path_monthly = make_monthly_pyramid_path(gcm_identifier)
-    pyramid_path_annual = make_annual_pyramid_path(gcm_identifier)
-
-    return (
-        gcm_grid_spec,
-        obs_identifier,
-        gcm_identifier,
-        pyramid_path_daily,
-        pyramid_path_monthly,
-        pyramid_path_annual,
-    )
-
-
-@task(
-    checkpoint=True,
-    result=XpersistResult(
-        CacheStore(config.get('storage.intermediate.uri')),
-        serializer='xarray.zarr',
-    ),
-    target=make_interpolated_obs_path,
-)
-def coarsen_and_interpolate_obs_task(
-    obs, train_period, predict_period, variables, gcm, scenario, chunking_approach, bbox, **kwargs
-):
-    """
-    Coarsen the observation dataset to the grid of the GCM model specified in inputs then
-    interpolate back into the observation grid. Rechunk the final output according to chunking approach.
-    Parameters
-    ----------
-    obs: str
-        Name of obs dataset
-    gcm: str
-        Name of GCM model
-    training_period_start: str
-        Start year of training/historical period
-    training_period_end: str
-        End year of training/historical period
-    variables: List[str]
-        List of variables to get in obs dataset
-    chunking_approach: str
-        'full_space', 'full_time', or None
-    **kwargs: Dict
-        Other arguments to be used in generating the target path
-    Returns
-    -------
-    ds_obs_interpolated_rechunked: xr.Dataset
-        An observation dataset that has been coarsened, interpolated back to original grid, and then rechunked.
-    """
-    # get obs
-    ds_obs = return_obs(obs=obs, train_period=train_period, variable=variables, bbox=bbox)
-
-    # regrid to coarse scale
-    ds_obs_coarse = get_coarse_obs(
-        obs_ds=ds_obs,
-        gcm=gcm,
-        scenario=scenario,
-        variable=variables,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        **kwargs,
-    )
-
-    # interpolate to fine scale again
-    ds_obs_interpolated = regrid_ds(
-        ds=ds_obs_coarse,
-        target_grid_ds=ds_obs.isel(time=0).chunk({'lat': -1, 'lon': -1}),
-        chunking_approach="full_space",
-    )
-
-    # rechunked to final output chunking approach if needed
-    ds_obs_interpolated_rechunked = rechunk_zarr_array_with_caching(
-        zarr_array=ds_obs_interpolated,
-        output_path=None,
-        chunking_approach=chunking_approach,
-    )
-
-    return ds_obs_interpolated_rechunked
-
-
-@task(
-    checkpoint=True,
-    result=XpersistResult(
-        CacheStore(config.get('storage.intermediate.uri')),
-        serializer='xarray.zarr',
-    ),
-    target=make_interpolated_gcm_path,
-)
-def interpolate_gcm_task(
-    obs: str,
-    gcm: str,
-    scenario: str,
-    train_period: slice,
-    predict_period: slice,
-    variables: Union[str, List[str]],
-    chunking_approach: str,
-    bbox,
-):
-    """
-    Interpolate the GCM dataset to the grid of the observation dataset.
-    Rechunk the final output according to chunking approach.
-    Parameters
-    ----------
-    obs: str
-        Name of obs dataset
-    gcm: str
-        Name of the GCM model
-    scenario: str
-        Name of the emission scenario
-    training_period: slice
-        Training/historical period bounds
-    predict_period: slice
-        Prediction (historical and/or future) period bounds
-    variables: List[str]
-        List of variables to get in obs dataset
-    chunking_approach: str
-        'full_space', 'full_time', or None
-
-    Returns
-    -------
-    ds_gcm_interpolated_rechunked: xr.Dataset
-        The GCM dataset that has been interpolated to the obs grid then rechunked.
-    """
-    # get gcm in full space chunks
-    ds_gcm_full_space = get_gcm(
-        gcm=gcm,
-        scenario=scenario,
-        variables=variables,
-        train_period=train_period,
-        predict_period=predict_period,
-        chunking_approach="full_space",
-        cache_within_rechunk=False,
-        bbox=bbox,
-    )
-
-    # get obs as a template
-    ds_obs = return_obs(obs=obs, train_period=train_period, variable=variables, bbox=bbox)
-
-    # interpolate gcm to obs resolution
-    ds_gcm_interpolated = regrid_ds(
-        ds=ds_gcm_full_space,
-        target_grid_ds=ds_obs.isel(time=0).load(),
-        chunking_approach="full_space",
-    )
-
-    # rechunked to final output chunking approach if needed
-    ds_gcm_interpolated_rechunked = rechunk_zarr_array_with_caching(
-        zarr_array=ds_gcm_interpolated,
-        output_path=None,
-        chunking_approach=chunking_approach,
-    )
-
-    return ds_gcm_interpolated_rechunked
-
-
-@task(
-    log_stdout=True,
-    result=XpersistResult(
-        CacheStore(config.get('storage.intermediate.uri')),
-        serializer='xarray.zarr',
-    ),
-    target=make_bias_corrected_obs_path,
-)
-def bias_correct_obs_task(
-    ds_obs: xr.Dataset, method: str, bc_kwargs: Optional[Dict[str, Any]] = None, **kwargs
-) -> xr.DataArray:
-    """
-    Bias correct observation data according to methods and kwargs.
-    Parameters
-    ----------
-    ds_obs : xr.Dataset
-        Observation dataset
-    method : str
-        Bias correction method to be used.
-    bc_kwargs: dict or None
-        Keyword arguments to be used with the bias correction method
-    kwargs: dict
-        Other arguments to be used in generating the target path
-    Returns
-    -------
-    ds_obs_bias_corrected : xr.Dataset
-        Bias corrected observation dataset
-    """
-    kws = default_none_kwargs(bc_kwargs, copy=True)
-    bias_corrected = bias_correct_obs_by_method(
-        da_obs=ds_obs, method=method, bc_kwargs=kws
-    ).to_dataset(dim="variable")
-
-    return bias_corrected
-
-
-@task(
-    result=XpersistResult(
-        CacheStore(config.get('storage.intermediate.uri')),
-        serializer='xarray.zarr',
-    ),
-    target=make_bias_corrected_gcm_path,
-)
-def bias_correct_gcm_task(
-    ds_gcm: xr.Dataset,
-    ds_obs: xr.Dataset,
-    historical_period: slice,
-    method: str,
-    bc_kwargs: Optional[Dict[str, Any]] = None,
-    **kwargs,
-) -> xr.DataArray:
-    """
-    Bias correct gcm data to the provided observation data according to methods and kwargs.
-
-    Parameters
-    ----------
-    ds_gcm : xr.Dataset
-        GCM dataset to be bias corrected
-    ds_obs : xr.Dataset
-        Observation dataset to bias correct to
-    historical_period_start : str
-        Start year of the historical/training period
-    historical_period_end : str
-        End year of the historical/training period
-    method : str
-        Bias correction method to be used.
-    bc_kwargs: dict or None
-        Keyword arguments to be used with the bias correction method
-
-    Returns
-    -------
-    ds_gcm_bias_corrected : xr.Dataset
-        Bias corrected GCM dataset
-    """
-    kws = default_none_kwargs(bc_kwargs, copy=True)
-
-    for v in ds_gcm.data_vars:
-        assert v in ds_obs.data_vars
-    bias_corrected = bias_correct_gcm_by_method(
-        da_gcm=ds_gcm,
-        da_obs=ds_obs,
-        historical_period=historical_period,
-        method=method,
-        bc_kwargs=kws,
-    ).to_dataset(dim="variable")
-
-    return bias_corrected
-
-
-def target_grid_obs_ds_task(obs_ds: xr.Dataset) -> xr.Dataset:
-    target_grid_obs_ds = obs_ds.isel(time=0)
-    return target_grid_obs_ds
diff --git a/cmip6_downscaling/tasks/pyramid.py b/cmip6_downscaling/tasks/pyramid.py
deleted file mode 100644
index fefd0c27..00000000
--- a/cmip6_downscaling/tasks/pyramid.py
+++ /dev/null
@@ -1,107 +0,0 @@
-from datetime import timedelta
-
-import dask
-import datatree as dt
-import fsspec
-import xarray as xr
-from carbonplan_data.metadata import get_cf_global_attrs
-from carbonplan_data.utils import set_zarr_encoding
-from ndpyramid import pyramid_regrid
-from prefect import task
-
-PIXELS_PER_TILE = 128
-
-
-def _load_coords(ds: xr.Dataset) -> xr.Dataset:
-    '''Helper function to explicitly load all dataset coordinates'''
-    for var, da in ds.coords.items():
-        ds[var] = da.load()
-    return ds
-
-
-def _postprocess(dt: dt.DataTree, levels: int, other_chunks: dict = None) -> dt.DataTree:
-    '''Postprocess data pyramid
-
-    Adds multiscales metadata and sets Zarr encoding
-
-    Parameters
-    ----------
-    dt : dt.DataTree
-        Input data pyramid
-    levels : int
-        Number of levels in pyramid
-    other_chunks : dict
-        Chunks for non-spatial dims
-
-    Returns
-    -------
-    dt.DataTree
-        Updated data pyramid with metadata / encoding set
-    '''
-    chunks = {"x": PIXELS_PER_TILE, "y": PIXELS_PER_TILE}
-    if other_chunks is not None:
-        chunks.update(other_chunks)
-
-    for level in range(levels):
-        slevel = str(level)
-        dt.ds.attrs['multiscales'][0]['datasets'][level]['pixels_per_tile'] = PIXELS_PER_TILE
-
-        # set dataset chunks
-        dt[slevel].ds = dt[slevel].ds.chunk(chunks)
-        if 'date_str' in dt[slevel].ds:
-            dt[slevel].ds['date_str'] = dt[slevel].ds['date_str'].chunk(-1)
-
-        # set dataset encoding
-        dt[slevel].ds = set_zarr_encoding(
-            dt[slevel].ds, codec_config={"id": "zlib", "level": 1}, float_dtype="float32"
-        )
-        for var in ['time', 'time_bnds']:
-            if var in dt[slevel].ds:
-                dt[slevel].ds[var].encoding['dtype'] = 'int32'
-
-    # set global metadata
-    dt.ds.attrs.update(**get_cf_global_attrs())
-    return dt
-
-
-@task(
-    log_stdout=True,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=10,
-    retry_delay=timedelta(seconds=10),
-)
-def regrid(ds: xr.Dataset, levels: int = 2, uri: str = None, other_chunks: dict = None) -> str:
-    '''Task to create a data pyramid from an xarray Dataset
-
-    Parameters
-    ----------
-    ds : xr.Dataset
-        Input Dataset
-    levels : int, optional
-        Number of levels in pyramid, by default 2
-    uri : str, optional
-        Path to write output data pyamid to, by default None
-    other_chunks : dict
-        Chunks for non-spatial dims
-    '''
-    with dask.config.set(scheduler='single-threaded'):
-
-        ds.coords['date_str'] = ds['time'].dt.strftime('%Y-%m-%d').astype('S10')
-
-        ds = _load_coords(ds)
-
-        mapper = fsspec.get_mapper(uri)
-        # Question: Is this needed or will prefect handle this for us?
-        if '.zmetadata' in mapper:
-            return
-
-        # create
-        dta = pyramid_regrid(ds, target_pyramid=None, levels=levels)
-
-        # postprocess
-        dta = _postprocess(dta, levels, other_chunks=other_chunks)
-
-        # write to uri
-        dta.to_zarr(mapper, mode='w')
-
-    return uri
diff --git a/cmip6_downscaling/workflows/bcsd_flow.py b/cmip6_downscaling/workflows/bcsd_flow.py
deleted file mode 100644
index 5ecd5b70..00000000
--- a/cmip6_downscaling/workflows/bcsd_flow.py
+++ /dev/null
@@ -1,431 +0,0 @@
-from datetime import timedelta
-
-from prefect import Flow, Parameter, task
-from xpersist import CacheStore
-from xpersist.prefect.result import XpersistResult
-
-from cmip6_downscaling import config, runtimes
-from cmip6_downscaling.analysis.analysis import annual_summary, monthly_summary, run_analyses
-from cmip6_downscaling.methods.bcsd import (
-    fit_and_predict,
-    get_coarse_obs,
-    get_interpolated_obs,
-    get_interpolated_prediction,
-    get_spatial_anomalies,
-    postprocess_bcsd,
-    rechunked_interpolated_prediciton_task_full_time,
-    return_coarse_obs_full_time,
-    return_gcm_predict_rechunked,
-    return_gcm_train_full_time,
-    return_obs,
-)
-from cmip6_downscaling.tasks import cleanup, pyramid
-from cmip6_downscaling.tasks.common_tasks import (
-    build_bbox,
-    build_time_period_slices,
-    path_builder_task,
-    target_grid_obs_ds_task,
-)
-from cmip6_downscaling.workflows.paths import (
-    make_annual_summary_path,
-    make_bcsd_output_path,
-    make_bias_corrected_path,
-    make_coarse_obs_path,
-    make_gcm_predict_path,
-    make_interpolated_obs_path,
-    make_interpolated_prediction_path_full_space,
-    make_monthly_summary_path,
-    make_rechunked_gcm_path,
-    make_return_obs_path,
-    make_spatial_anomalies_path,
-)
-
-runtime = runtimes.get_runtime()
-
-intermediate_cache_store = CacheStore(
-    config.get("storage.intermediate.uri"),
-    storage_options=config.get("storage.intermediate.storage_options"),
-)
-
-results_cache_store = CacheStore(
-    config.get("storage.results.uri"),
-    storage_options=config.get("storage.results.storage_options"),
-)
-
-serializer_dump_kwargs = config.get('storage.xpersist_overwrite')
-# Transform Functions into Tasks -----------------------------------------------------------
-
-
-return_obs_task = task(
-    return_obs,
-    result=XpersistResult(
-        intermediate_cache_store,
-        serializer="xarray.zarr",
-        serializer_dump_kwargs=serializer_dump_kwargs,
-    ),
-    target=make_return_obs_path,
-)
-
-
-get_coarse_obs_task = task(
-    get_coarse_obs,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-    result=XpersistResult(
-        intermediate_cache_store,
-        serializer="xarray.zarr",
-        serializer_dump_kwargs=serializer_dump_kwargs,
-    ),
-    target=make_coarse_obs_path,
-)
-
-target_grid_obs_ds_task = task(
-    target_grid_obs_ds_task,
-)
-
-
-interpolated_obs_task = task(
-    get_interpolated_obs,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-    result=XpersistResult(
-        intermediate_cache_store,
-        serializer="xarray.zarr",
-        serializer_dump_kwargs=serializer_dump_kwargs,
-    ),
-    target=make_interpolated_obs_path,
-)
-
-
-interpolated_prediction_task = task(
-    get_interpolated_prediction,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-    result=XpersistResult(
-        intermediate_cache_store,
-        serializer="xarray.zarr",
-        serializer_dump_kwargs=serializer_dump_kwargs,
-    ),
-    target=make_interpolated_prediction_path_full_space,
-)
-
-rechunked_interpolated_prediciton_task_full_time_task = task(
-    rechunked_interpolated_prediciton_task_full_time,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-)
-
-get_spatial_anomalies_task = task(
-    get_spatial_anomalies,
-    log_stdout=True,
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-    result=XpersistResult(
-        intermediate_cache_store,
-        serializer="xarray.zarr",
-        serializer_dump_kwargs=serializer_dump_kwargs,
-    ),
-    target=make_spatial_anomalies_path,
-)
-
-
-return_coarse_obs_full_time_task = task(
-    return_coarse_obs_full_time,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-)
-
-
-return_gcm_train_full_time_task = task(
-    return_gcm_train_full_time,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-    result=XpersistResult(
-        intermediate_cache_store,
-        serializer="xarray.zarr",
-        serializer_dump_kwargs=serializer_dump_kwargs,
-    ),
-    target=make_rechunked_gcm_path,
-)
-
-return_gcm_predict_rechunked_task = task(
-    return_gcm_predict_rechunked,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-    result=XpersistResult(
-        intermediate_cache_store,
-        serializer="xarray.zarr",
-        serializer_dump_kwargs=serializer_dump_kwargs,
-    ),
-    target=make_gcm_predict_path,
-)
-
-fit_and_predict_task = task(
-    fit_and_predict,
-    log_stdout=True,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-    result=XpersistResult(
-        intermediate_cache_store,
-        serializer="xarray.zarr",
-        serializer_dump_kwargs=serializer_dump_kwargs,
-    ),
-    target=make_bias_corrected_path,
-)
-postprocess_bcsd_task = task(
-    postprocess_bcsd,
-    tags=['dask-resource:TASKSLOTS=1'],
-    max_retries=5,
-    retry_delay=timedelta(seconds=5),
-    log_stdout=True,
-    result=XpersistResult(
-        results_cache_store, serializer="xarray.zarr", serializer_dump_kwargs=serializer_dump_kwargs
-    ),
-    target=make_bcsd_output_path,
-)
-
-monthly_summary_task = task(
-    monthly_summary,
-    log_stdout=True,
-    result=XpersistResult(
-        results_cache_store, serializer="xarray.zarr", serializer_dump_kwargs=serializer_dump_kwargs
-    ),
-    target=make_monthly_summary_path,
-)
-
-annual_summary_task = task(
-    annual_summary,
-    result=XpersistResult(
-        results_cache_store, serializer="xarray.zarr", serializer_dump_kwargs=serializer_dump_kwargs
-    ),
-    target=make_annual_summary_path,
-)
-
-
-# # Main Flow -----------------------------------------------------------
-
-with Flow(
-    name="bcsd",
-    storage=runtime.storage,
-    run_config=runtime.run_config,
-    executor=runtime.executor,
-) as bcsd_flow:
-    obs = Parameter("obs")
-    gcm = Parameter("gcm")
-    scenario = Parameter("scenario")
-    variable = Parameter("variable")
-
-    bbox = build_bbox(
-        latmin=Parameter("latmin"),
-        latmax=Parameter("latmax"),
-        lonmin=Parameter("lonmin"),
-        lonmax=Parameter("lonmax"),
-    )
-    train_period = build_time_period_slices(Parameter('train_period'))
-    predict_period = build_time_period_slices(Parameter('predict_period'))
-    (
-        gcm_grid_spec,
-        obs_identifier,
-        gcm_identifier,
-        pyramid_path_daily,
-        pyramid_path_monthly,
-        pyramid_path_annual,
-    ) = path_builder_task(
-        obs=obs,
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-    )
-    if config.get('run_options.cleanup_flag') is True:
-        cleanup.run_rsfip(gcm_identifier, obs_identifier)
-
-    # preprocess_bcsd_tasks(s):
-
-    obs_ds = return_obs_task(
-        obs=obs,
-        variable=variable,
-        train_period=train_period,
-        bbox=bbox,
-        obs_identifier=obs_identifier,
-    )
-
-    coarse_obs_ds = get_coarse_obs_task(
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        obs_identifier=obs_identifier,
-        gcm_grid_spec=gcm_grid_spec,
-        chunking_approach='full_space',
-        upstream_tasks=[obs_ds],
-    )
-
-    target_grid_obs_ds = target_grid_obs_ds_task(obs_ds)
-
-    interpolated_obs_ds = interpolated_obs_task(
-        target_grid_ds=target_grid_obs_ds,
-        gcm_grid_spec=gcm_grid_spec,
-        chunking_approach='full_space',
-        obs_identifier=obs_identifier,
-        upstream_tasks=[coarse_obs_ds, target_grid_obs_ds],
-    )
-
-    spatial_anomalies_ds = get_spatial_anomalies_task(
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        obs_identifier=obs_identifier,
-        gcm_grid_spec=gcm_grid_spec,
-        chunking_approach='full_space',
-        upstream_tasks=[interpolated_obs_ds, obs_ds],
-    )
-
-    # the next three tasks prepare the inputs required by bcsd
-    coarse_obs_full_time_path = return_coarse_obs_full_time_task(
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        obs_identifier=obs_identifier,
-        gcm_grid_spec=gcm_grid_spec,
-        upstream_tasks=[coarse_obs_ds],
-    )
-
-    gcm_train_subset_full_time_ds = return_gcm_train_full_time_task(
-        coarse_obs_full_time_path=coarse_obs_full_time_path,
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        gcm_identifier=gcm_identifier,
-        upstream_tasks=[coarse_obs_full_time_path],
-    )
-
-    gcm_predict_rechunked_ds = return_gcm_predict_rechunked_task(
-        gcm_train_subset_full_time_ds=gcm_train_subset_full_time_ds,
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        gcm_identifier=gcm_identifier,
-        upstream_tasks=[gcm_train_subset_full_time_ds],
-    )
-
-    # fit and predict tasks(s):
-    bias_corrected_ds = fit_and_predict_task(
-        gcm_train_subset_full_time_ds=gcm_train_subset_full_time_ds,
-        gcm_predict_rechunked_ds=gcm_predict_rechunked_ds,
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        chunking_approach='full_time',
-        obs_identifier=obs_identifier,
-        gcm_grid_spec=gcm_grid_spec,
-        gcm_identifier=gcm_identifier,
-        upstream_tasks=[gcm_train_subset_full_time_ds, gcm_predict_rechunked_ds],
-    )
-
-    interpolated_prediction_ds = interpolated_prediction_task(
-        target_grid_obs_ds=target_grid_obs_ds,
-        gcm_identifier=gcm_identifier,
-        gcm_grid_spec=gcm_grid_spec,
-        upstream_tasks=[bias_corrected_ds, target_grid_obs_ds],
-    )
-
-    rechunked_interpolated_prediction_path = rechunked_interpolated_prediciton_task_full_time_task(
-        interpolated_prediction_ds=interpolated_prediction_ds,
-        gcm_identifier=gcm_identifier,
-        gcm_grid_spec=gcm_grid_spec,
-        upstream_tasks=[interpolated_prediction_ds],
-    )
-
-    # # postprocess_bcsd_task(s):
-    postprocess_bcsd_ds = postprocess_bcsd_task(
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        obs_identifier=obs_identifier,
-        gcm_identifier=gcm_identifier,
-        gcm_grid_spec=gcm_grid_spec,
-        upstream_tasks=[spatial_anomalies_ds, rechunked_interpolated_prediction_path],
-    )
-
-    monthly_summary_ds = monthly_summary_task(
-        postprocess_bcsd_ds,
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        gcm_identifier=gcm_identifier,
-    )
-
-    annual_summary_ds = annual_summary_task(
-        postprocess_bcsd_ds,
-        gcm=gcm,
-        scenario=scenario,
-        variable=variable,
-        train_period=train_period,
-        predict_period=predict_period,
-        bbox=bbox,
-        gcm_identifier=gcm_identifier,
-    )
-
-    analysis_location = run_analyses(
-        {
-            'gcm_identifier': gcm_identifier,
-            'obs_identifier': obs_identifier,
-            'gcm_grid_spec': gcm_grid_spec,
-            'result_dir': config.get('storage.results.uri'),
-            'intermediate_dir': config.get('storage.intermediate.uri'),
-            'var': variable,
-            'gcm': gcm,
-            'scenario': scenario,
-        },
-        web_blob=config.get('storage.web_results.blob'),
-        bbox=bbox,
-        train_period=train_period,
-        predict_period=predict_period,
-        upstream_tasks=[postprocess_bcsd_ds],
-    )
-
-    pyramid_location_daily = pyramid.regrid(
-        postprocess_bcsd_ds, uri=config.get('storage.results.uri') + pyramid_path_daily, levels=4
-    )
-
-    pyramid_location_monthly = pyramid.regrid(
-        monthly_summary_ds, uri=config.get('storage.results.uri') + pyramid_path_monthly, levels=4
-    )
-
-    pyramid_location_annual = pyramid.regrid(
-        annual_summary_ds, uri=config.get('storage.results.uri') + pyramid_path_annual, levels=4
-    )
diff --git a/cmip6_downscaling/workflows/bcsd_params.json b/cmip6_downscaling/workflows/bcsd_params.json
index 45fe799e..fffb7482 100644
--- a/cmip6_downscaling/workflows/bcsd_params.json
+++ b/cmip6_downscaling/workflows/bcsd_params.json
@@ -1,18 +1,22 @@
 {
+  "method": "bcsd",
   "obs": "ERA5",
-  "gcm": "MIROC6",
+  "model": "MIROC6",
   "scenario": "ssp370",
+  "grid_label": "gn",
+  "table_id": "day",
+  "member": "r1i1p1f1",
   "variable": "tasmax",
   "train_period": [
-    "1990",
-    "1990"
+    "1981",
+    "2010"
   ],
   "predict_period": [
-    "2020",
-    "2020"
+    "1981",
+    "2099"
   ],
-  "latmin": "40",
-  "latmax": "50",
-  "lonmin": "80",
-  "lonmax": "90"
+  "latmin": "90",
+  "latmax": "-90",
+  "lonmin": "-180",
+  "lonmax": "180"
 }
diff --git a/cmip6_downscaling/workflows/paths.py b/cmip6_downscaling/workflows/paths.py
deleted file mode 100644
index 867e0c49..00000000
--- a/cmip6_downscaling/workflows/paths.py
+++ /dev/null
@@ -1,788 +0,0 @@
-import pathlib
-from typing import Optional, Tuple
-
-from cmip6_downscaling import config
-from cmip6_downscaling.workflows.utils import BBox
-
-
-def build_obs_identifier(
-    obs: str,
-    variable: str,
-    train_period: slice,
-    bbox: BBox,
-    **kwargs,
-) -> str:
-    """
-    Build the common identifier for observation related data: the same pattern is used for: 1) chunked raw obs, 2) coarsened obs, and 3) coarsened then interpolated obs
-
-    Parameters
-    ----------
-    obs: str
-        Name of obs dataset
-    variable: str
-        Name of the variable used in obs and gcm dataset (including features and label)
-    train_period: slice
-        Start and end year slice of training/historical period. Ex: slice('1990','1990')
-    bbox: BBox
-        dataclass containing the latmin,latmax,lonmin,lonmax. Class can be found in utils.
-
-
-    Returns
-    -------
-    identifier : str
-        string to be used in obs related paths as specified by the params
-    """
-    if isinstance(variable, str):
-        variable = [variable]
-    var_string = '_'.join(sorted(variable))
-
-    obs_identifier = config.get('storage.obs_identifier_template').format(
-        obs=obs,
-        train_period=f'{train_period.start}_{train_period.stop}',
-        bbox=bbox,
-        variable=var_string,
-    )
-
-    return obs_identifier
-
-
-def build_gcm_identifier(
-    gcm: str,
-    scenario: str,
-    variable: str,
-    train_period: slice,
-    predict_period: slice,
-    bbox: BBox,
-    **kwargs,
-) -> str:
-    """
-    Build the common identifier for GCM related data
-
-    Parameters
-    ----------
-    gcm : str
-        Name of the GCM model
-    scenario : str
-        Name of the future emission scenario to load
-    variable : str
-        Name of the variable used in obs and gcm dataset (including features and label)
-    train_period : slice
-        Start and end year slice of training/historical period. Ex: slice('1990', '1990')
-    predict_period: slice
-        Start and end year slice of prediction period. Ex: slice('2020', '2040')
-    bbox: BBox
-        Bounding box including latmin,latmax,lonmin,lonmax.
-
-    Returns
-    -------
-    identifier : str
-        string to be used in gcm related paths as specified by the params
-    """
-    if isinstance(variable, str):
-        variable = [variable]
-    var_string = '_'.join(sorted(variable))
-
-    gcm_identifier = config.get('storage.gcm_identifier_template').format(
-        gcm=gcm,
-        scenario=scenario,
-        variable=var_string,
-        train_period=f'{train_period.start}_{train_period.stop}',
-        predict_period=f'{predict_period.start}_{predict_period.stop}',
-        bbox=bbox,
-    )
-    return gcm_identifier
-
-
-def make_rechunked_obs_path(
-    chunking_approach: Optional[str] = '', obs_identifier: Optional[str] = None, **kwargs
-) -> str:
-    """Build the path for rechunked observation
-
-    Parameters
-    ----------
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-    chunking_approach : str
-        How the data is chunked. Eg. `full_time`, `full_space`, `matched` or None
-
-    Returns
-    -------
-    rechunked_obs_path: str
-        Path to which rechunked observation defined by the parameters should be stored
-    """
-    if obs_identifier is None:
-        obs_identifier = build_obs_identifier(**kwargs)
-
-    return f"rechunked_obs/{obs_identifier}{chunking_approach}.zarr"
-
-
-def make_coarse_obs_path(
-    chunking_approach: Optional[str] = '',
-    obs_identifier: Optional[str] = '',
-    gcm_grid_spec: Optional[str] = '',
-    **kwargs,
-) -> str:
-    """Build the path for coarsened observation
-
-    Parameters
-    ----------
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-
-    Returns
-    -------
-    coarse_obs_path: str
-        Path to which coarsened observation defined by the parameters should be stored
-    """
-
-    return f"coarsened_obs/{obs_identifier}{chunking_approach}_{gcm_grid_spec}.spatialzarr"
-
-
-def make_coarse_obs_path_full_space(
-    obs_identifier: Optional[str] = '', gcm_grid_spec: Optional[str] = '', **kwargs
-) -> str:
-    """Build the path for coarsened observation for full_time
-
-    Parameters
-    ----------
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-
-    Returns
-    -------
-    coarse_obs_path: str
-        Path to which coarsened observation defined by the parameters should be stored
-    """
-
-    return f"coarsened_obs/{obs_identifier}full_space_{gcm_grid_spec}.zarr"
-
-
-def make_coarse_obs_path_full_time(
-    obs_identifier: Optional[str] = '', gcm_grid_spec: Optional[str] = '', **kwargs
-) -> str:
-    """Build the path for coarsened observation for full_time
-
-    Parameters
-    ----------
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-
-    Returns
-    -------
-    coarse_obs_path: str
-        Path to which coarsened observation defined by the parameters should be stored
-    """
-
-    return f"coarsened_obs/{obs_identifier}full_time_{gcm_grid_spec}.zarr"
-
-
-def make_interpolated_obs_path(
-    gcm_grid_spec: str, chunking_approach: str, obs_identifier: Optional[str] = None, **kwargs
-) -> str:
-    """Build the path for coarsened observation that has then been interpolated back to the observation grid
-
-    Parameters
-    ----------
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-    chunking_approach : str
-        How the data is chunked. Eg. `full_time`, `full_space`, `matched` or None
-    gcm_grid_spec: str
-        Output of get_gcm_grid_spec. String to identify the GCM grid
-
-    Returns
-    -------
-    interpolated_obs_path: str
-        Path to which interpolated observation defined by the parameters should be stored
-    """
-    if obs_identifier is None:
-        obs_identifier = build_obs_identifier(**kwargs)
-    return f"interpolated_obs/{obs_identifier}{chunking_approach}_{gcm_grid_spec}.zarr"
-
-
-def make_interpolated_gcm_path(
-    obs: str, chunking_approach: Optional[str] = '', gcm_identifier: Optional[str] = None, **kwargs
-) -> str:
-    """Build the path for GCM that has then been interpolated back to the observation grid
-
-    Parameters
-    ----------
-    obs: str
-        name of obs dataset
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-
-    Returns
-    -------
-    interpolated_obs_path: str
-        Path to which interpolated observation defined by the parameters should be stored
-    """
-    if gcm_identifier is None:
-        gcm_identifier = build_gcm_identifier(**kwargs)
-
-    return f"interpolated_gcm/{gcm_identifier}{chunking_approach}_{obs}.zarr"
-
-
-def make_bias_corrected_obs_path(
-    method: str,
-    gcm_grid_spec: str,
-    obs_identifier: Optional[str] = None,
-    chunking_approach: Optional[str] = None,
-    **kwargs,
-):
-    """
-    Path to save "bias_corrected" or pre-processed (e.g. normalized or quantile transformed) observation dataset.
-    The dataset might be coarsened hence the need for gcm_grid_spec in the path.
-
-    Parameters
-    ----------
-    method: str
-        String to identify bias correction method used
-    gcm_grid_spec: str
-        Output of get_gcm_grid_spec. String to identify the GCM grid
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-    chunking_approach : str
-        How the data is chunked. Eg. `full_time`, `full_space`, `matched` or None
-
-    Returns
-    -------
-    bias_corrected_obs_path: str
-        Path to which bias corrected observation defined by the parameters should be stored
-    """
-    if obs_identifier is None:
-        obs_identifier = build_obs_identifier(**kwargs)
-    if chunking_approach is None:
-        chunking_approach = ''
-
-    return f"bias_corrected_obs/{obs_identifier}{chunking_approach}_{gcm_grid_spec}_{method}.zarr"
-
-
-def make_rechunked_gcm_path(
-    gcm_identifier: Optional[str] = None, chunking_approach: Optional[str] = None, **kwargs
-) -> str:
-    """Build the path for rechunked GCM
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-
-    Returns
-    -------
-    rechunked_gcm_path: str
-        Path to which rechunked GCM defined by the parameters should be stored
-    """
-    if gcm_identifier is None:
-        gcm_identifier = build_gcm_identifier(**kwargs)
-
-    if chunking_approach is None:
-        chunking_approach = ''
-
-    return f"rechunked_gcm/{gcm_identifier}{chunking_approach}.zarr"
-
-
-def make_bias_corrected_gcm_path(
-    method: str,
-    gcm_identifier: Optional[str] = None,
-    chunking_approach: Optional[str] = None,
-    **kwargs,
-):
-    """
-    Path to save bias_corrected or pre-processed (e.g. normalized or quantile transformed) GCM dataset
-
-    Parameters
-    ----------
-    method: str
-        String to identify bias correction method used
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-    chunking_approach : str
-        How the data is chunked. Eg. `full_time`, `full_space`, `matched` or None
-
-    Returns
-    -------
-    bias_corrected_obs_path: str
-        Path to which bias corrected observation defined by the parameters should be stored
-    """
-    if gcm_identifier is None:
-        gcm_identifier = build_gcm_identifier(**kwargs)
-    if chunking_approach is None:
-        chunking_approach = ''
-
-    return f"bias_corrected_gcm/{gcm_identifier}{chunking_approach}_{method}.zarr"
-
-
-def make_gard_predict_output_path(
-    gcm_identifier: str, bias_correction_method: str, model_type: str, label: str, **kwargs
-):
-    """
-    Path to save GARD prediction output
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-    bias_correction_method: str
-        String to identify bias correction method used
-    model_type : str
-        The type of GARD model (`AnalogRegression`, `PureAnalog` or `PureRegression`)
-    label: str
-        The variable being predicted
-
-    Returns
-    -------
-    gard_prediction_output_path: str
-        Path to which GARD prediction output defined by the parameters should be stored
-    """
-    return f"gard_pred_output/{gcm_identifier}_{bias_correction_method}_{model_type}_{label}.zarr"
-
-
-def make_scrf_path(
-    obs: str,
-    label: str,
-    start_year: str,
-    end_year: str,
-    **kwargs,
-):
-    """
-    Path where spatially-temporally correlated random fields (SCRF) are saved.
-
-    Parameters
-    ----------
-    obs : str
-        Name of observation dataset
-    label: str
-        The variable being predicted
-    start_year: str
-        Start year of the dataset
-    end_year: str
-        End year of the dataset
-
-    Returns
-    -------
-    scrf_path : str
-        Path of SCRF
-    """
-    return f"scrf/{obs}_{label}_{start_year}_{end_year}.zarr"
-
-
-def make_gard_post_processed_output_path(
-    gcm_identifier: str, bias_correction_method: str, model_type: str, label: str, **kwargs
-):
-    """
-    Path to save GARD post processed output
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-    bias_correction_method: str
-        String to identify bias correction method used
-    model_type : str
-        The type of GARD model (`AnalogRegression`, `PureAnalog` or `PureRegression`)
-    label: str
-        The variable being predicted
-
-    Returns
-    -------
-    gard_post_processed_output_path: str
-        Path to which GARD post processed output defined by the parameters should be stored
-    """
-    return f"gard/{gcm_identifier}{bias_correction_method}_{model_type}_{label}.zarr"
-
-
-def make_epoch_trend_path(
-    gcm_identifier: str, day_rolling_window: int, year_rolling_window: int, **kwargs
-):
-    """
-    Path to save epoch trend from MACA
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-    day_rolling_window: int
-        Number of days to include when calculating the rolling average
-    year_rolling_window: int
-        Number of years to include when calculating the rolling average
-
-    Returns
-    -------
-    epoch_trend_path : str
-        Path to save the epoch trend
-    """
-    return f'gcm_epoch_trend/{gcm_identifier}_{day_rolling_window}_day_{year_rolling_window}_year_rolling.zarr'
-
-
-def make_epoch_adjusted_gcm_path(
-    gcm_identifier: str, day_rolling_window: int, year_rolling_window: int, **kwargs
-):
-    """
-    Path to save epoch adjusted GCM data from MACA
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-    day_rolling_window: int
-        Number of days to include when calculating the rolling average
-    year_rolling_window: int
-        Number of years to include when calculating the rolling average
-
-    Returns
-    -------
-    epoch_adjusted_gcm_path : str
-        Path to save epoch adjusted GCM data from MACA
-    """
-    return f'epoch_adjusted_gcm/{gcm_identifier}_{day_rolling_window}_day_{year_rolling_window}_year_rolling.zarr'
-
-
-def make_epoch_adjusted_downscaled_gcm_path(
-    gcm_identifier: str,
-    label: str,
-    subdomain_bound: Optional[Tuple[float, float, float, float]] = (0, -90, 360, 90),
-    **kwargs,
-):
-    """
-    Path to save the downscaled GCM data where the epoch trend has not been repalced. Includes subdomain
-    bound because the downscaling step in MACA is performed on subdomains.
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-    label: str
-        The variable being predicted
-    subdomain_boudn: Tuple of floats
-        min_lon, min_lat, max_lon, max_lat
-
-    Returns
-    -------
-    maca_epoch_adjusted_downscaled_output_path : str
-        Path to save the downscaled GCM data where the epoch trend has not been repalced.
-    """
-    min_lon, min_lat, max_lon, max_lat = subdomain_bound
-    return f'maca_epoch_adjusted_downscaled_output/{gcm_identifier}_{label}_{min_lon:0.1f}_{min_lat:0.1f}_{max_lon:0.1f}_{max_lat:0.1f}.zarr'
-
-
-def make_epoch_replaced_downscaled_gcm_path(
-    gcm_identifier: str, day_rolling_window: int, year_rolling_window: int, **kwargs
-):
-    """
-    Path to save the downscaled GCM data where the epoch trend has already been repalced.
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-    day_rolling_window: int
-        Number of days to include when calculating the rolling average
-    year_rolling_window: int
-        Number of years to include when calculating the rolling average
-
-    Returns
-    -------
-    epoch_replaced_gcm : str
-        Path to save the downscaled GCM data where the epoch trend has already been repalced.
-    """
-    return f'epoch_replaced_gcm/{gcm_identifier}_{day_rolling_window}_day_{year_rolling_window}_year_rolling.zarr'
-
-
-def make_maca_output_path(gcm_identifier: str, label: str, **kwargs):
-    """
-    Path to save final MACA output after fine scale bias correction.
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the GCM dataset used
-    label: str
-        The variable being predicted
-
-    Returns
-    -------
-    maca_output_path : str
-        MACA output path
-    """
-    return f"maca_output/{gcm_identifier}_{label}.zarr"
-
-
-def make_daily_pyramid_path(gcm_identifier: str, **kwargs) -> str:
-    """Build the path for the daily pyramid
-
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-
-    Returns
-    -------
-    pyramid_path : str
-        Path to pyramid location
-    """
-    return f"/pyramid_daily/{gcm_identifier}.pyr"
-
-
-def make_monthly_pyramid_path(gcm_identifier: str, **kwargs) -> str:
-    """Build the path for the monthly pyramid
-
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-
-    Returns
-    -------
-    pyramid_path : str
-        Path to the monthly pyramid
-    """
-    return f"/pyramid_monthly/{gcm_identifier}.pyr"
-
-
-def make_annual_pyramid_path(gcm_identifier: str, **kwargs) -> str:
-    """Build the path for the annual pyramid
-
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-
-    Returns
-    -------
-    pyramid_path : str
-        Path to the annual pyramid
-    """
-
-    return f"/pyramid_annual/{gcm_identifier}.pyr"
-
-
-# ---addl bcsd paths
-
-
-def make_return_obs_path(obs_identifier: str, **kwargs) -> str:
-    """Build the path for return obs task
-
-    Parameters
-    ----------
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-
-    Returns
-    -------
-    obs_ds path : str
-        Path to bcsd obs_ds file location
-    """
-    return f"obs_ds/{obs_identifier}.zarr"
-
-
-def make_spatial_anomalies_path(obs_identifier: str, gcm_grid_spec: str, **kwargs) -> str:
-    """Build the path for spatial anomalies
-
-    Parameters
-    ----------
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-
-    Returns
-    -------
-    spatial_anomalies_path : str
-        Path to bcsd spatial anomalies file location
-    """
-    return f"spatial_anomalies/{obs_identifier}{gcm_grid_spec}.zarr"
-
-
-def make_spatial_anomalies_rechunked_full_time_path(
-    obs_identifier: str, gcm_grid_spec: str, **kwargs
-):
-    """
-        Parameters
-    ----------
-    obs_identifier : str
-        Output from build_obs_identifier. String to identify the observation dataset used
-
-    Returns
-    -------
-    spatial_anomalies_rechunked_full_time_path : str
-        Path to rechunked spatial anomalies file location
-    """
-    return f"spatial_anomalies_rechunked_full_time/{obs_identifier}{gcm_grid_spec}.zarr"
-
-
-def make_gcm_train_subset_path(gcm_identifier: str = None, **kwargs) -> str:
-    return f"gcm_train_subset/{gcm_identifier}.zarr"
-
-
-def make_gcm_predict_subset_path(gcm_identifier: str, **kwargs) -> str:
-    return f"gcm_predict_subset/{gcm_identifier}.zarr"
-
-
-def make_gcm_predict_path(gcm_identifier: str = None, **kwargs) -> str:
-    """Build the path for the gcm predict dataset
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-
-    Returns
-    -------
-    gcm_predict_path : str
-        Path to gcm_predict ds file location
-    """
-    return f"gcm_predict/{gcm_identifier}.zarr"
-
-
-def make_bias_corrected_path(gcm_identifier: str = None, **kwargs) -> str:
-    """Build the path for the bias corrected bcsd dataset
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-
-    Returns
-    -------
-    bias_corrected_path : str
-        Path to bcsd bias_corrected ds file location
-    """
-    return f"bias_corrected/{gcm_identifier}.zarr"
-
-
-def make_bcsd_output_path(gcm_identifier: str = None, **kwargs) -> str:
-    """Build the path for the bcsd output dataset
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-
-    Returns
-    -------
-    bcsd_output_path : str
-        Path to bcsd output ds file location
-    """
-    return f"bcsd_output/{gcm_identifier}.zarr"
-
-
-def make_monthly_summary_path(gcm_identifier: str = None, **kwargs) -> str:
-    """Build the path for the monthly summary dataset
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-
-    Returns
-    -------
-    bcsd_output_monthly_path : str
-        Path to monthly bcsd output ds file location
-    """
-    return f"bcsd_output_monthly/{gcm_identifier}.zarr"
-
-
-def make_annual_summary_path(gcm_identifier: str = None, **kwargs) -> str:
-    """Build the path for the annual summary dataset
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-
-    Returns
-    -------
-    bcsd_output_annual_path : str
-        Path to annual bcsd output ds file location
-    """
-    return f"bcsd_output_annual/{gcm_identifier}.zarr"
-
-
-def make_interpolated_prediction_path_full_space(
-    gcm_identifier: str = None, gcm_grid_spec: str = None, **kwargs
-) -> str:
-    """Build the path for the interpolated prediction in full space
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-    gcm_grid_spec: str
-        Output of get_gcm_grid_spec. String to identify the GCM grid
-
-    Returns
-    -------
-    interpolated_prediction_path_full_space_location : str
-        Path to interpolated obs full space ds file location
-    """
-    return f"interpolated_prediction_path_full_space/{gcm_identifier}_{gcm_grid_spec}.zarr"
-
-
-def make_interpolated_prediction_path_full_time(
-    gcm_identifier: str = None, gcm_grid_spec: str = None, **kwargs
-) -> str:
-    """Build the path for the interpolated prediction in full time
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-    gcm_grid_spec: str
-        Output of get_gcm_grid_spec. String to identify the GCM grid
-
-    Returns
-    -------
-    interpolated_prediction_path_full_time_location : str
-        Path to interpolated obs full time ds file location
-    """
-    return f"interpolated_prediction_full_time/{gcm_identifier}_{gcm_grid_spec}.zarr"
-
-
-def make_rechunked_spatial_anomalies_path_full_time(
-    gcm_identifier: str = None, gcm_grid_spec: str = None, **kwargs
-) -> str:
-    """Build the path for the rechunked spatial anomlaies in full time
-
-    Parameters
-    ----------
-    gcm_identifier : str
-        Output from build_gcm_identifier. String to identify the gcm dataset used
-    gcm_grid_spec: str
-        Output of get_gcm_grid_spec. String to identify the GCM grid
-
-    Returns
-    -------
-    make_rechunked_spatial_anomalies_path_full_time_location : str
-        Path to rechunked_spatial_anomalies full time ds file location
-    """
-    return f"rechunked_spatial_anomalies_full_time/{gcm_identifier}_{gcm_grid_spec}.zarr"
-
-
-def get_notebook_paths(
-    identifier: str,
-) -> Tuple[pathlib.PosixPath, pathlib.PosixPath, pathlib.PosixPath]:
-    """Create the paths for where the notebook template is, where
-    the output notebook will go, and where the converted html will go
-
-    Parameters
-    ----------
-    identifier : str
-        unique identifier including the gcm and domain
-
-    Returns
-    -------
-    Tuple[pathlib.PosixPath, pathlib.PosixPath, pathlib.PosixPath]
-        local paths to the three analysis files
-    """
-    # just using this metrics as a shortcut to getting the location
-    from cmip6_downscaling.analysis import metrics
-
-    path = pathlib.PosixPath(metrics.__file__)
-    template_path = path.parent / 'analyses_template.ipynb'
-    executed_path = path.parent / f'analyses_{identifier}.ipynb'
-    executed_html_path = path.parent / f'analyses_{identifier}.html'
-    return template_path, executed_path, executed_html_path
diff --git a/cmip6_downscaling/workflows/utils.py b/cmip6_downscaling/workflows/utils.py
index cf64300d..ffde3145 100644
--- a/cmip6_downscaling/workflows/utils.py
+++ b/cmip6_downscaling/workflows/utils.py
@@ -1,4 +1,3 @@
-import os
 import random
 import re
 import string
@@ -20,42 +19,6 @@
 from dataclasses import dataclass
 
 
-@dataclass
-class BBox:
-    """Lat/lon bounding box"""
-
-    latmin: float = -90.0
-    latmax: float = 90.0
-    lonmin: float = -180.0
-    lonmax: float = 180
-
-    @property
-    def lat_slice(self) -> slice:
-        return slice(float(self.latmax), float(self.latmin))
-
-    @property
-    def lon_slice(self) -> slice:
-        return slice(float(self.lonmin), float(self.lonmax))
-
-    def __str__(self) -> str:
-        return f'{self.latmin}_{self.latmax}_{self.lonmin}_{self.lonmax}'
-
-
-def get_store(prefix, account_key=None):
-    """helper function to create a zarr store"""
-
-    if account_key is None:
-        account_key = os.environ.get("BLOB_ACCOUNT_KEY", None)
-
-    store = zarr.storage.ABSStore(
-        "carbonplan-downscaling",
-        prefix=prefix,
-        account_name="carbonplan",
-        account_key=account_key,
-    )
-    return store
-
-
 def subset_dataset(
     ds: xr.Dataset,
     variable: str,
@@ -175,37 +138,37 @@ def delete_chunks_encoding(ds: Union[xr.Dataset, xr.DataArray]):
             del ds[coord].encoding["chunks"]
 
 
-def lon_to_180(ds):
-    '''Converts longitude values to (-180, 180)
+# def lon_to_180(ds):
+#     '''Converts longitude values to (-180, 180)
 
-    Parameters
-    ----------
-    ds : xr.Dataset
-        Input dataset with `lon` coordinate
+#     Parameters
+#     ----------
+#     ds : xr.Dataset
+#         Input dataset with `lon` coordinate
 
-    Returns
-    -------
-    xr.Dataset
-        Copy of `ds` with updated coordinates
+#     Returns
+#     -------
+#     xr.Dataset
+#         Copy of `ds` with updated coordinates
 
-    See also
-    --------
-    cmip6_preprocessing.preprocessing.correct_lon
-    '''
+#     See also
+#     --------
+#     cmip6_preprocessing.preprocessing.correct_lon
+#     '''
 
-    ds = ds.copy()
+#     ds = ds.copy()
 
-    lon = ds["lon"].where(ds["lon"] < 180, ds["lon"] - 360)
-    ds = ds.assign_coords(lon=lon)
+#     lon = ds["lon"].where(ds["lon"] < 180, ds["lon"] - 360)
+#     ds = ds.assign_coords(lon=lon)
 
-    if not (ds["lon"].diff(dim="lon") > 0).all():
-        ds = ds.reindex(lon=np.sort(ds["lon"].data))
+#     if not (ds["lon"].diff(dim="lon") > 0).all():
+#         ds = ds.reindex(lon=np.sort(ds["lon"].data))
 
-    if "lon_bounds" in ds.variables:
-        lon_b = ds["lon_bounds"].where(ds["lon_bounds"] < 180, ds["lon_bounds"] - 360)
-        ds = ds.assign_coords(lon_bounds=lon_b)
+#     if "lon_bounds" in ds.variables:
+#         lon_b = ds["lon_bounds"].where(ds["lon_bounds"] < 180, ds["lon_bounds"] - 360)
+#         ds = ds.assign_coords(lon_bounds=lon_b)
 
-    return ds
+#     return ds
 
 
 def make_rechunker_stores(
diff --git a/notebooks/deepsd_data_prep.ipynb b/notebooks/deepsd_data_prep.ipynb
index 14bac391..e68de8bd 100644
--- a/notebooks/deepsd_data_prep.ipynb
+++ b/notebooks/deepsd_data_prep.ipynb
@@ -1227,7 +1227,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.9"
+   "version": "3.9.10"
   }
  },
  "nbformat": 4,
diff --git a/setup.cfg b/setup.cfg
index 8a813b3e..d608c4ae 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -7,7 +7,7 @@ select = B,C,E,F,W,T4,B9
 
 [isort]
 
-known_third_party=ESMF,azure,carbonplan,carbonplan_data,cartopy,climate_indices,dask,dask_kubernetes,datatree,donfig,fsspec,intake,intake_esm,matplotlib,ndpyramid,numba,numpy,pandas,papermill,prefect,pytest,rechunker,regionmask,scipy,seaborn,setuptools,skdownscale,sklearn,xarray,xarray_schema,xclim,xesmf,xpersist,zarr
+known_third_party=ESMF,carbonplan,carbonplan_data,cartopy,climate_indices,dask,dask_kubernetes,datatree,donfig,fsspec,intake,intake_esm,matplotlib,ndpyramid,numba,numpy,pandas,pkg_resources,prefect,pytest,rechunker,regionmask,scipy,seaborn,setuptools,skdownscale,sklearn,upath,xarray,xarray_schema,xclim,xesmf,xpersist,zarr
 
 multi_line_output=3
 include_trailing_comma=True


Test patch:
diff --git a/cmip6_downscaling/tests/methods/common/test_tasks.py b/cmip6_downscaling/tests/methods/common/test_tasks.py
new file mode 100644
index 00000000..7199cf14
--- /dev/null
+++ b/cmip6_downscaling/tests/methods/common/test_tasks.py
@@ -0,0 +1,174 @@
+from dataclasses import asdict
+
+import numpy as np
+import pytest
+import xarray as xr
+from upath import UPath
+from xarray_schema import DataArraySchema, DatasetSchema
+
+from cmip6_downscaling import config
+
+# TODO: figure out how to make this a fixture or pytest config option
+config.set(
+    {
+        'storage.scratch.uri': '/tmp/cmip6_downscaling_tests/scratch',
+        'storage.intermediate.uri': '/tmp/cmip6_downscaling_tests/intermediate',
+        'storage.results.uri': '/tmp/cmip6_downscaling_tests/results',
+    }
+)
+
+from cmip6_downscaling.methods.common.containers import RunParameters
+from cmip6_downscaling.methods.common.tasks import (
+    get_experiment,
+    get_obs,
+    make_run_parameters,
+    rechunk,
+    regrid,
+)
+
+params = [
+    {
+        'method': 'bcsd',
+        'obs': 'ERA5',
+        'model': 'MIROC6',
+        'member': 'r1i1p1f1',
+        'grid_label': 'gn',
+        'table_id': 'day',
+        'scenario': 'ssp370',
+        'variable': 'tasmax',
+        'latmin': '-2',
+        'latmax': '2',
+        'lonmin': '14.5',
+        'lonmax': '18.5',
+        'train_dates': ['1980', '1981'],
+        'predict_dates': ['2050', '2051'],
+    },
+    {
+        'method': 'gard',
+        'obs': 'ERA5',
+        'model': 'MIROC6',
+        'member': 'r1i1p1f1',
+        'grid_label': 'gn',
+        'table_id': 'day',
+        'scenario': 'ssp370',
+        'variable': 'tasmax',
+        'latmin': '-2',
+        'latmax': '2',
+        'lonmin': '14.5',
+        'lonmax': '18.5',
+        'train_dates': ['1980', '1981'],
+        'predict_dates': ['2050', '2051'],
+    },
+]
+
+
+rechunk_params = [
+    {
+        'chunking_method': 'full_space',
+        'chunking_schema': {'time': (2359, 561), 'lat': -1, 'lon': -1},
+    },
+    {'chunking_method': 'full_time', 'chunking_schema': {'time': -1, 'lat': 25, 'lon': (33, 20)}},
+]
+
+
+@pytest.fixture(scope="module", params=params)
+def run_parameters(request):
+    return RunParameters(**request.param)
+
+
+@pytest.mark.parametrize("params", params)
+def test_make_run_parameters(params):
+    rps = make_run_parameters.run(**params)
+    assert isinstance(rps, RunParameters)
+    assert asdict(rps) == params
+
+
+def check_global_attrs(ds):
+    for key in ['history', 'hostname', 'institution', 'source', 'title', 'username', 'version']:
+        assert key in ds.attrs
+
+
+def test_get_obs(run_parameters):
+    obs_path = get_obs.run(run_parameters)
+
+    ds = xr.open_zarr(obs_path)
+
+    check_global_attrs(ds)
+    schema = DatasetSchema(
+        {
+            run_parameters.variable: DataArraySchema(
+                dtype=np.floating, name=run_parameters.variable, dims=['time', 'lat', 'lon']
+            )
+        }
+    )
+    schema.validate(ds)
+
+
+@pytest.mark.xfail
+def test_get_experiment(run_parameters):
+    get_experiment_path = get_experiment.run(run_parameters, time_subset='train_period')
+
+    ds = xr.open_zarr(get_experiment_path)
+    check_global_attrs(ds)
+    schema = DatasetSchema(
+        {
+            run_parameters.variable: DataArraySchema(
+                dtype=np.floating, name=run_parameters.variable, dims=['time', 'lat', 'lon']
+            )
+        }
+    )
+    schema.validate(ds)
+
+
+def test_regrid(tmp_path):
+
+    pytest.importorskip('xesmf')
+
+    ds = xr.tutorial.open_dataset('air_temperature').chunk({'time': 10, 'lat': -1, 'lon': -1})
+    source_path = UPath(tmp_path) / UPath('regrid_source.zarr')
+    ds.to_zarr(source_path)
+
+    target_ds = ds.isel(time=0).coarsen(lat=4, lon=4, boundary='trim').mean()
+    target_grid_path = UPath(tmp_path) / UPath('regrid_target.zarr')
+    target_ds.to_zarr(target_grid_path)
+
+    actual_path = regrid.run(source_path, target_grid_path)
+    actual_ds = xr.open_zarr(actual_path)
+
+    check_global_attrs(actual_ds)
+    expected_shape = (ds.dims['time'], target_ds.dims['lat'], target_ds.dims['lon'])
+    schema = DatasetSchema(
+        {
+            'air': DataArraySchema(
+                dtype=np.floating, name='air', dims=['time', 'lat', 'lon'], shape=expected_shape
+            )
+        }
+    )
+    schema.validate(actual_ds)
+
+
+@pytest.mark.parametrize('rechunk_params', rechunk_params)
+@pytest.mark.xfail
+def test_rechunk(rechunk_params, tmp_path):
+    # TODO Add testing parameterization to check full_space (done), full_time (done) and template match
+    ds = xr.tutorial.open_dataset('air_temperature').chunk({'time': 100, 'lat': 10, 'lon': 10})
+    source_path = str(tmp_path) + "/" + 'rechunk.zarr'
+    ds.to_zarr(source_path, mode='w')
+
+    actual_path = rechunk.run(
+        source_path,
+        pattern=rechunk_params['chunking_method'],
+    )
+    actual_ds = xr.open_zarr(actual_path)
+    print(actual_ds)
+
+    expected_chunks = rechunk_params['chunking_schema']
+    schema = DatasetSchema(
+        {
+            'air': DataArraySchema(
+                dtype=np.floating, name='air', dims=['time', 'lat', 'lon'], chunks=expected_chunks
+            )
+        }
+    )
+
+    schema.validate(actual_ds)
diff --git a/cmip6_downscaling/tests/test_prefect_tasks.py b/cmip6_downscaling/tests/test_prefect_tasks.py
deleted file mode 100644
index d245ffc9..00000000
--- a/cmip6_downscaling/tests/test_prefect_tasks.py
+++ /dev/null
@@ -1,23 +0,0 @@
-import numpy as np
-import pytest
-import xarray as xr
-
-from cmip6_downscaling.tasks.common_tasks import to_standard_calendar
-
-
-@pytest.fixture
-def da_noleap(val=1.0):
-    time = xr.cftime_range(start='2020-01-01', end='2020-12-31', freq='1D', calendar='noleap')
-    return xr.DataArray(
-        val,
-        dims=['lat', 'lon', 'time'],
-        coords={'lat': np.arange(19, 56, 1), 'lon': np.arange(-133, -61, 2), 'time': time},
-    )
-
-
-def test_to_standard_calendar(da_noleap):
-
-    da_std = to_standard_calendar.run(da_noleap)
-    assert da_noleap.sizes['time'] == 365
-    assert da_std.sizes['time'] == 366
-    assert not da_std.isnull().any().item()
diff --git a/cmip6_downscaling/tests/test_utils.py b/cmip6_downscaling/tests/test_utils.py
index 3aea9dea..f1ad9aea 100644
--- a/cmip6_downscaling/tests/test_utils.py
+++ b/cmip6_downscaling/tests/test_utils.py
@@ -2,7 +2,7 @@
 import pytest
 import xarray as xr
 
-from cmip6_downscaling.workflows.utils import lon_to_180
+from cmip6_downscaling.data.utils import lon_to_180, to_standard_calendar
 
 
 def create_test_ds(xname, yname, zname, xlen, ylen, zlen):
@@ -23,6 +23,16 @@ def create_test_ds(xname, yname, zname, xlen, ylen, zlen):
     return ds
 
 
+@pytest.fixture
+def da_noleap(val=1.0):
+    time = xr.cftime_range(start='2020-01-01', end='2020-12-31', freq='1D', calendar='noleap')
+    return xr.DataArray(
+        val,
+        dims=['lat', 'lon', 'time'],
+        coords={'lat': np.arange(19, 56, 1), 'lon': np.arange(-133, -61, 2), 'time': time},
+    )
+
+
 @pytest.mark.parametrize(
     "shift",
     [
@@ -42,3 +52,11 @@ def test_lon_to_180(shift):
     assert ds_lon_corrected.lon.min() < -1
     assert ds_lon_corrected.lon.max() <= 180
     assert (ds_lon_corrected.lon.diff(dim='lon') > 0).all()
+
+
+def test_to_standard_calendar(da_noleap):
+
+    da_std = to_standard_calendar(da_noleap)
+    assert da_noleap.sizes['time'] == 365
+    assert da_std.sizes['time'] == 366
+    assert not da_std.isnull().any().item()