workflow/2_TC_application.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "410f15c3-784d-4547-9dc4-ab424c898332",
   "metadata": {},
   "source": [
    "# Triple Collocation Uncertainty Analysis\n",
    "\n",
    "Now that we have all of our monthly ET datasets spatially collocated from the [regridding notebook](1_regrid.ipynb), we are ready to perform a Triple Collocation (TC) analysis on the common date ranges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cab12df-6303-4efe-9729-c39b2b0c96e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.xarray\n",
    "import holoviews as hv\n",
    "import panel as pn\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import itertools\n",
    "import warnings\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a370a119-c6d8-4063-b501-53672e816e32",
   "metadata": {},
   "source": [
    "First, we will run in the Extended Collocation notebook to create our TC function. (We use the EC function as it can do TC and runs each spatial point simultaneously. See the [EC notebook](../TC/EC_function.ipynb) and the [TC notebook](../TC/TC_function.ipynb) for details on the TC method.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ff608c6-411d-4af9-b135-75ccd9bc6a04",
   "metadata": {},
   "outputs": [],
   "source": [
    "%run ../TC/EC_function.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10712275-5f82-4a6b-878f-999873578ac0",
   "metadata": {},
   "source": [
    "## Combine Data Sets in Xarray\n",
    "\n",
    "Next, we need to load in our data sets and limit them to a common date range. Since we need at least three data sets to utilized TC, we will restrict the data ranges of all data sets to have the beginning date of the third oldest starting date and ending data of the third most recent ending date. This choice allows us to save memory usage, while also utilizing the largest amount of data. For triplets with a more restricted date range, due to one data set having a smaller date range, we will limit the date range further at the time of the TC computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2c070ad-347a-4f15-8540-3f180f0d4e19",
   "metadata": {},
   "outputs": [],
   "source": [
    "files = [\n",
    "    \"../Data/ssebop/ssebop_aet_regridded.nc\",\n",
    "    \"../Data/gleam/gleam_aet.nc\",\n",
    "    \"../Data/era5/era5_aet_regridded.nc\",\n",
    "    \"../Data/nldas/nldas_aet_regridded.nc\",\n",
    "    \"../Data/terraclimate/terraclimate_aet_regridded.nc\",\n",
    "    \"../Data/wbet/wbet_aet_regridded.nc\",\n",
    "]\n",
    "dataset_name = [\"SSEBop\", \"GLEAM\", \"ERA5\", \"NLDAS\", \"TerraClimate\", \"WBET\"]\n",
    "\n",
    "date_ranges = {}\n",
    "for file, name in zip(files, dataset_name):\n",
    "    ds_temp = xr.open_dataset(\n",
    "        file, engine=\"netcdf4\", chunks={\"lon\": -1, \"lat\": -1, \"time\": -1}\n",
    "    )\n",
    "    date_ranges[name] = [ds_temp.time.min().values, ds_temp.time.max().values]\n",
    "\n",
    "# Take the third oldest start and third most recent end dates\n",
    "date_range = [\n",
    "    np.sort(np.array(list(date_ranges.values()))[:, 0])[2],\n",
    "    np.sort(np.array(list(date_ranges.values()))[:, 1])[-3],\n",
    "]\n",
    "date_range"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78d9e767-3d97-4f32-a0a3-36c2c36b8920",
   "metadata": {},
   "source": [
    "Using the date range, we can now combine all of the data sets into a single `xarray.DataSet` for easy computations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8c50291-2fc7-4996-852b-8e28822214c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def preprocess(ds):\n",
    "    \"\"\"\n",
    "    Keep only the specified time range for each file.\n",
    "    \"\"\"\n",
    "    return ds.sel(time=slice(date_range[0], date_range[1]))\n",
    "\n",
    "\n",
    "ds = xr.open_mfdataset(\n",
    "    files,\n",
    "    engine=\"netcdf4\",\n",
    "    preprocess=preprocess,\n",
    "    combine=\"nested\",\n",
    "    concat_dim=\"dataset_name\",\n",
    ")\n",
    "ds = ds.assign_coords({\"dataset_name\": dataset_name})\n",
    "ds.dataset_name.attrs[\"description\"] = \"Dataset name\"\n",
    "\n",
    "# Need time as first index for TC computation\n",
    "ds = ds.transpose(\"time\", ...)\n",
    "# The data set is less than 1GiB, so let's read it into memory vs keeping as a dask array\n",
    "ds = ds.compute()\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3501d7d-950f-4e60-9332-c3cd4dc657dd",
   "metadata": {},
   "source": [
    "## Time Series Exploration\n",
    "In case we want to explore the time series of each pixel later, let's make an interactive figure where we can select the latitude and longitude and plot the time series. As we may also want to explore how the time series vary with season, let's add the functionality to select the full year or just a certain season. (Seasons will be denoted by the first letter of each month within the season. For example winter contains December, January, and February. So, it will be denoted by `DJF`. The full year with all seasons will be denoted with `All`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2f7f805-15dc-44f4-9041-46357b42005b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_timseries(lat=40.125, lon=-100.125, season=\"All\"):\n",
    "    location_map = ds.aet.isel(time=0, dataset_name=2)\n",
    "    location_map = location_map * 0\n",
    "    location_map.loc[{\"lat\": lat, \"lon\": lon}] = 1\n",
    "    if season == \"All\":\n",
    "        ds_season = ds\n",
    "    else:\n",
    "        ds_season = ds.isel(time=(ds.time.dt.season == season))\n",
    "\n",
    "    plt = location_map.hvplot(\n",
    "        geo=True,\n",
    "        coastline=True,\n",
    "        title=\"Time Series Location  (Red dot indicates current pixel)\",\n",
    "        colorbar=False,\n",
    "        cmap=\"kr\",\n",
    "    ).opts(frame_height=250) + ds_season.sel(lat=lat, lon=lon, method=\"nearest\").hvplot(\n",
    "        x=\"time\", groupby=\"dataset_name\", title=\"Datasets' ET Time Series\"\n",
    "    ).overlay().opts(legend_position=\"right\", frame_height=250)\n",
    "\n",
    "    return plt\n",
    "\n",
    "\n",
    "lat_widget = pn.widgets.FloatSlider(\n",
    "    name=\"lat\",\n",
    "    start=ds.lat.min().item(),\n",
    "    end=ds.lat.max().item(),\n",
    "    step=0.25,\n",
    "    value=40.125,\n",
    ")\n",
    "lon_widget = pn.widgets.FloatSlider(\n",
    "    name=\"lon\",\n",
    "    start=ds.lon.min().item(),\n",
    "    end=ds.lon.max().item(),\n",
    "    step=0.25,\n",
    "    value=-100.125,\n",
    ")\n",
    "season_widget = pn.widgets.Select(\n",
    "    name=\"season\", value=\"All\", options=[\"All\", \"DJF\", \"MAM\", \"JJA\", \"SON\"]\n",
    ")\n",
    "\n",
    "bound_plot = pn.bind(\n",
    "    create_timseries, lat=lat_widget, lon=lon_widget, season=season_widget\n",
    ")\n",
    "\n",
    "pn.Column(pn.Row(pn.Column(lat_widget, lon_widget), season_widget), bound_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f235cab-dd47-4927-8ccc-d060a2381f28",
   "metadata": {},
   "source": [
    "## TC Estimation\n",
    "\n",
    "Time to compute the TC uncertainty estimates. To do that, we first need to decide on data sets that have \"independent\" errors in order to group them together into TC sets.\n",
    "\n",
    "Here is a table of the data and method used for calculating each ET data set:\n",
    "\n",
    "| Data Sets      | Reference Data Type | Calculation Method | Date Range       | Resolution       | Input Data |\n",
    "| -------------- | ------------------- | ------------------ | ---------------- | ---------------- | ---------- |\n",
    "| SSEBop         | Ex situ             | Energy balance     | 2001/01-2022/12 | 1 km             | **STRM** elevn; **PRISM** Ta; **MODIS** Ts, emissivity, albedo, and NDVI; **GDAS** ETo |\n",
    "| GLEAM v3.7b    | Ex situ             | Energy balance     | 2003/01-2022/12 | 0.25$^{\\circ}$   | **CERES** radiation; **TMPA** precip; **AIRS** Ta; **GLOBSNOW** snow-water equiv; **CCI** vegetation optical depth; **GLDAS** and **CCI** Soil moisture; **MODIS** GVCF (global vegetation continuous fields); **IGBP-DIS** soil properties; **CGLFRD** lightning flash rate for rainfall inference |\n",
    "| ERA5-Land      | Reanalysis          | Energy balance     | 1950/01-2022/12 | 0.1$^{\\circ}$    | **CHTESSEL** Land surface model using model cycle Cy45r1 (2018) |\n",
    "| NLDAS-2 (Noah) | Reanalysis          | Energy balance     | 1979/02-2022/12 | 0.125$^{\\circ}$  | **NARR** (North American Regional Reanalysis) atmospheric forcing data; **PRISM** precip |\n",
    "| TerraClimate   | Interpolated        | Water balance      | 1958/01-2022/12 | 0.0416$^{\\circ}$ | **WorldClim** Ta, vapor, precip, solar radiation, wind (Uses **MODIS** Ts, cloud cover; **STRM** elevn); **CRU** Ts4.0, Tmax, Tmin, vapor, precip, Ta; **JRA-55** Ta, vapor, precip, radiation, wind |\n",
    "| WBET           | Interpolated        | Water balance      | 1895/10-2018/09 | 800 m            | **PRISM** precip, mean Ta, max Ta, min Ta; **USGS** water use irrigaion, national elevation dataset, NWIS gage II discharge; **EROS** land cover (1938-1999); **Landsat** NLCD land cover (2000-2018); **gridMT** wind; **Koppen-Geiger** climate classification; **Fenneman & Johnson** physiographic province classification; **EPA** level III ecoregions; **STATSGO2** soil saturated hydraulic conductivity, porosity, field capacity, thickness, available water capacity |\n",
    "\n",
    "From this table, we can group the data sets into measurement systems that \"should\" be independent:\n",
    "\n",
    "1) Ex situ\n",
    "2) Reanalysis\n",
    "3) Interpolated\n",
    "\n",
    "This give us 8 different possible combinations of datasets. However, since the computation is fast and resulting TC error estimates will be small in memory (~65MiB), we will just compute all 20 combinations and can filter them out later if we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41ddf1e9-ddfc-4a57-a08a-06767d9f7f49",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate a list of the combinations\n",
    "combos = list(itertools.combinations(dataset_name, 3))\n",
    "combos = [list(combo) for combo in combos]\n",
    "combos"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4adec689-7c6b-42a6-9c5e-fd361cd520b3",
   "metadata": {},
   "source": [
    "Since we have data sets with different date ranges, we will need to trim the date ranges here before computing the TC error variances. This will be slightly complicated. So, let's make it the date range selection its own function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f82c2479-0bca-451f-bfa1-d81aca05c85a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def common_date_range(ds, combo):\n",
    "    \"\"\"Return the common date slice of the datasets.\"\"\"\n",
    "    old_common_date = []\n",
    "    recent_common_date = []\n",
    "    for name in combo:\n",
    "        old_common_date.append(date_ranges[name][0])\n",
    "        recent_common_date.append(date_ranges[name][1])\n",
    "\n",
    "    return slice(np.max(old_common_date), np.min(recent_common_date))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7f1bdec-224f-4bbb-963c-5597675b786d",
   "metadata": {},
   "source": [
    "Now that we have the ability to select the common date range, let's compute the TC error standard deviations. We will do this for each season independently along with the full year. (The season and full year will be denoted with the monthly abbreviations or `All` as above.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68004cc6-3c96-4b16-9503-4b5c13ed312d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We want to ignore all of the sqrt and log warnings with negative values\n",
    "warnings.filterwarnings(\"ignore\", category=RuntimeWarning)\n",
    "\n",
    "# Create list of seasons\n",
    "seasons = [\"All\"] + list(np.unique(ds.time.dt.season))\n",
    "\n",
    "tc_est = []\n",
    "tc_est_season = []\n",
    "for combo in combos:\n",
    "    for season in seasons:\n",
    "        if season == \"All\":\n",
    "            ds_season = ds\n",
    "        else:\n",
    "            ds_season = ds.isel(time=(ds.time.dt.season == season))\n",
    "\n",
    "        ds_combo = ds_season.sel(time=common_date_range(ds, combo), dataset_name=combo)\n",
    "\n",
    "        tc_covar, snr_temp = ec_covar_multi(\n",
    "            ds_combo.aet.data, corr_sets=[1, 2, 3], return_snr=True\n",
    "        )\n",
    "\n",
    "        tc_est_season.append(\n",
    "            xr.Dataset(\n",
    "                data_vars={\n",
    "                    \"error\": (\n",
    "                        [\"dataset_combo\", \"season\", \"combo_idx\", \"lat\", \"lon\"],\n",
    "                        np.sqrt(np.diagonal(tc_covar)).transpose((2, 0, 1))[\n",
    "                            None, None, ...\n",
    "                        ],\n",
    "                    ),\n",
    "                    \"snr\": (\n",
    "                        [\"dataset_combo\", \"season\", \"combo_idx\", \"lat\", \"lon\"],\n",
    "                        # Perform a logarithm as a simple method for setting\n",
    "                        # negatives to NaN\n",
    "                        (10 ** np.log10(snr_temp[None, None, ...])),\n",
    "                    ),\n",
    "                },\n",
    "                coords={\n",
    "                    \"dataset_combo\": [\" \".join(combo)],\n",
    "                    \"season\": [season],\n",
    "                    \"combo_idx\": [0, 1, 2],\n",
    "                    \"lat\": ds.lat,\n",
    "                    \"lon\": ds.lon,\n",
    "                },\n",
    "            )\n",
    "        )\n",
    "\n",
    "    tc_est.append(xr.concat(tc_est_season, dim=\"season\"))\n",
    "    tc_est_season = []\n",
    "\n",
    "tc_est = xr.concat(tc_est, dim=\"dataset_combo\")\n",
    "\n",
    "tc_est.error.attrs[\"description\"] = (\n",
    "    \"TC error standard deviation estimate for the dataset_combo triplet.\"\n",
    ")\n",
    "tc_est.snr.attrs[\"description\"] = (\n",
    "    \"TC unbiased SNR estimate for the dataset_combo triplet.\"\n",
    ")\n",
    "tc_est.dataset_combo.attrs[\"description\"] = \"Dataset combination used in TC evaluation.\"\n",
    "tc_est.combo_idx.attrs[\"description\"] = (\n",
    "    'Name index of \"dataset_combo\" coordinate associated with the data set.'\n",
    ")\n",
    "tc_est.season.attrs[\"description\"] = (\n",
    "    \"Season of the year given by the first letter of each month within the \"\n",
    "    'season. The full year is given by \"All\".'\n",
    ")\n",
    "tc_est.error.attrs[\"units\"] = \"mm.month-1\"\n",
    "tc_est = tc_est.compute()\n",
    "\n",
    "tc_est"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd78e15a-38cc-4365-8e2b-4620af280714",
   "metadata": {},
   "source": [
    "For convenience, let's rearrange the ET `DataSet` to be by data set rather than combination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b84031d5-3e41-4827-9ff4-b90868d7c74e",
   "metadata": {},
   "outputs": [],
   "source": [
    "tc_est_by_dataset = []\n",
    "est_pair = []\n",
    "for name in dataset_name:\n",
    "    idx_loc = np.char.find(tc_est.dataset_combo.data, name)\n",
    "    dataset_loc = np.where(idx_loc != -1)[0]\n",
    "    combos_single_dataset = tc_est.isel(dataset_combo=dataset_loc)\n",
    "\n",
    "    tc_est_dataset = []\n",
    "    for i in range(len(combos_single_dataset.dataset_combo.values)):\n",
    "        tc_est_single_dataset = combos_single_dataset.isel(dataset_combo=i)\n",
    "        idx = str(tc_est_single_dataset.dataset_combo.data).split(\" \").index(name)\n",
    "        tc_est_single_dataset = tc_est_single_dataset.sel(combo_idx=idx)\n",
    "        tc_est_single_dataset = tc_est_single_dataset.drop_vars(\n",
    "            [\"dataset_combo\", \"combo_idx\"]\n",
    "        )\n",
    "\n",
    "        tc_est_dataset.append(\n",
    "            xr.Dataset(\n",
    "                data_vars=tc_est_single_dataset,\n",
    "                coords={\n",
    "                    \"dataset_name\": name,\n",
    "                    \"est_idx\": i,\n",
    "                    \"lat\": ds.lat,\n",
    "                    \"lon\": ds.lon,\n",
    "                },\n",
    "            )\n",
    "        )\n",
    "\n",
    "    tc_est_dataset = xr.concat(tc_est_dataset, dim=\"est_idx\")\n",
    "\n",
    "    tc_est_by_dataset.append(tc_est_dataset)\n",
    "    est_pair.append(\n",
    "        [\n",
    "            combinations.replace(name, \"\").strip()\n",
    "            for combinations in tc_est.dataset_combo.data[dataset_loc]\n",
    "        ]\n",
    "    )\n",
    "\n",
    "tc_est_by_dataset = xr.concat(tc_est_by_dataset, dim=\"dataset_name\")\n",
    "\n",
    "tc_est_by_dataset = tc_est_by_dataset.assign_coords(\n",
    "    est_pair=([\"dataset_name\", \"est_idx\"], np.array(est_pair))\n",
    ")\n",
    "\n",
    "tc_est_by_dataset.dataset_name.attrs[\"description\"] = \"Dataset names.\"\n",
    "tc_est_by_dataset.est_idx.attrs[\"description\"] = (\n",
    "    \"Index of the other two data sets used in the TC triplet as contained in est_pair.\"\n",
    ")\n",
    "tc_est_by_dataset.est_pair.attrs[\"description\"] = (\n",
    "    \"Names of the other two data sets used in the TC triplet.\"\n",
    ")\n",
    "\n",
    "tc_est_by_dataset = tc_est_by_dataset.compute()\n",
    "\n",
    "# Save the results for later use in the 4_dataset_agreement and 5_regional notebooks\n",
    "if not os.path.isfile(\"../Data/TC_errs.nc\"):\n",
    "    _ = tc_est_by_dataset.to_netcdf(\n",
    "        path=\"../Data/TC_errs.nc\", format=\"NETCDF4\", engine=\"netcdf4\"\n",
    "    )\n",
    "\n",
    "tc_est_by_dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "242842d8-0674-49ec-be7f-de8971ba240c",
   "metadata": {},
   "source": [
    "Let's see how the resulting error estimates look along with the unbiased SNR of the data sets. This should give us an indication of how much non-noise information the data sets actually contain.\n",
    "\n",
    "> Since we suppressed the `sqrt` and `log` run time warnings, we can expect to see `NaN`s throughout the maps, where the TC calculation resulted in negative values. This works as intended as any negative error variances should be flagged as incorrect (i.e., this is done with `NaN`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c84d633-fc79-4323-86d6-4e88c8037510",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def tc_plts(dataset_name=\"SSEBop\", est_idx=0, season=\"All\"):\n",
    "    tc_data = tc_est_by_dataset.sel(\n",
    "        dataset_name=dataset_name, est_idx=est_idx, season=season\n",
    "    )\n",
    "\n",
    "    est_pairs = str(\n",
    "        tc_est_by_dataset.est_pair.sel(dataset_name=dataset_name, est_idx=est_idx).data\n",
    "    )\n",
    "    plt = tc_data.error.hvplot(\n",
    "        geo=True,\n",
    "        coastline=True,\n",
    "        clim=(0, 50),\n",
    "        title=\"Error Standard Deviation (other triplet datasets: \" + est_pairs + \")\",\n",
    "    ).opts(frame_width=500) + tc_data.snr.hvplot(\n",
    "        geo=True,\n",
    "        coastline=True,\n",
    "        clim=(0.1, 50),\n",
    "        cnorm=\"log\",\n",
    "        title=\"Unbiased SNR (other triplet datasets: \" + est_pairs + \")\",\n",
    "    ).opts(frame_width=500)\n",
    "\n",
    "    return plt\n",
    "\n",
    "\n",
    "dataset_name_widget = pn.widgets.Select(\n",
    "    name=\"dataset_name\",\n",
    "    value=\"SSEBop\",\n",
    "    options=list(tc_est_by_dataset.dataset_name.values),\n",
    ")\n",
    "est_idx_widget = pn.widgets.IntSlider(name=\"est_idx\", start=0, end=9, step=1, value=0)\n",
    "season_widget = pn.widgets.Select(\n",
    "    name=\"season\", value=\"All\", options=[\"All\", \"DJF\", \"MAM\", \"JJA\", \"SON\"]\n",
    ")\n",
    "\n",
    "bound_plot = pn.bind(\n",
    "    tc_plts,\n",
    "    dataset_name=dataset_name_widget,\n",
    "    season=season_widget,\n",
    "    est_idx=est_idx_widget,\n",
    ")\n",
    "\n",
    "pn.Column(dataset_name_widget, est_idx_widget, season_widget, bound_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13c9715b-4976-4e3a-9a01-e4e8d075138d",
   "metadata": {},
   "source": [
    "## TC Discussion\n",
    "\n",
    "Some of the datasets (Mainly GLEAM and NLDAS, and some ERA5 and WBET) have large swaths of `NaN` values caused by negative variances. This is typically caused by one of two things. (1) Covariances in the errors, which are assumed to not be present, or (2) two datasets have approximately order of magnitude larger error variances compared to the third.\n",
    "\n",
    "After carefully looking at the maps for each data set and pair, we can see that SSEBop and TerraClimate have the largest estimated errors across all combinations, followed by ERA5 and WBET, then GLEAM and NLDAS. Therefore, it appears that (2) is likely the culprit in causing the `NaN` swaths, since the `NaN` density appears to be inversely proportional to the estimate error standard deviation.\n",
    "\n",
    "To overcome these swaths of `NaN` values, let's combine each error variance to find the mean, median, and standard deviation for each computation, excluding the `NaN`s.\n",
    "\n",
    "> Note we use the error variances in the mean calculation and then take the sqrt for the mean error standard deviation. Using the error standard deviation in the mean is not statistically correct, as standard deviations should be added in quadrature (i.e., sum of variances). This does not effect the median as it is just a center value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf37246c-6dfa-4aba-a9a6-22e2bb1c1ec3",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_tc_est = np.sqrt(\n",
    "    (tc_est_by_dataset.error**2).mean(dim=\"est_idx\", skipna=True, keep_attrs=True)\n",
    ")\n",
    "mean_tc_est.name = \"mean_error\"\n",
    "mean_tc_est.attrs[\"description\"] = (\n",
    "    \"Mean TC error estimate across all combinations with other datasets.\"\n",
    ")\n",
    "mean_tc_est.attrs[\"units\"] = \"mm.month-1\"\n",
    "\n",
    "median_tc_est = tc_est_by_dataset.error.median(\n",
    "    dim=\"est_idx\", skipna=True, keep_attrs=True\n",
    ")\n",
    "median_tc_est.name = \"median_error\"\n",
    "median_tc_est.attrs[\"description\"] = (\n",
    "    \"Median TC error estimate across all combinations with other datasets.\"\n",
    ")\n",
    "median_tc_est.attrs[\"units\"] = \"mm.month-1\"\n",
    "\n",
    "std_tc_est = tc_est_by_dataset.error.std(\n",
    "    dim=\"est_idx\", ddof=1, skipna=True, keep_attrs=True\n",
    ")\n",
    "std_tc_est.name = \"std_error\"\n",
    "std_tc_est.attrs[\"description\"] = (\n",
    "    \"Standard Deviation of TC error estimate across all combinations with other datasets.\"\n",
    ")\n",
    "std_tc_est.attrs[\"units\"] = \"mm.month-1\"\n",
    "\n",
    "median_tc_snr = tc_est_by_dataset.snr.median(\n",
    "    dim=\"est_idx\", skipna=True, keep_attrs=True\n",
    ")\n",
    "median_tc_snr.name = \"median_snr\"\n",
    "median_tc_snr.attrs[\"description\"] = (\n",
    "    \"Median TC unbiased SNR estimate across all combinations with other datasets.\"\n",
    ")\n",
    "\n",
    "std_tc_snr = tc_est_by_dataset.snr.std(dim=\"est_idx\", skipna=True, keep_attrs=True)\n",
    "std_tc_snr.name = \"std_snr\"\n",
    "std_tc_snr.attrs[\"description\"] = (\n",
    "    \"Standard Deviation of TC unbiased SNR estimate across all \"\n",
    "    \"combinations with other datasets.\"\n",
    ")\n",
    "\n",
    "count_tc_est = np.isfinite(tc_est_by_dataset.error).sum(dim=\"est_idx\")\n",
    "count_tc_est.name = \"counts\"\n",
    "count_tc_est.attrs[\"description\"] = (\n",
    "    \"Number of datasets used in the average TC error \"\n",
    "    \"estimates (i.e., number of finite values in a given pixel).\"\n",
    ")\n",
    "count_tc_est.attrs[\"units\"] = \"counts\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39e024c6-2fe6-49d0-82e8-b207cb90217b",
   "metadata": {},
   "source": [
    "However, before we look at these average error estimates, let's check the fractional difference between each data set combination and the average to make sure the majority of the combinations are close to the average. This way we are not biasing the average by some extreme estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d720f1b5-4107-40e5-a516-fc2adc244e99",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "frac_diff = (tc_est_by_dataset.error - mean_tc_est) / mean_tc_est\n",
    "quant = frac_diff.quantile(np.linspace(0, 1, 101), dim=[\"lat\", \"lon\", \"est_idx\"])\n",
    "quant.name = \"Fractional Difference\"\n",
    "plt = (\n",
    "    quant.hvplot(groupby=[\"dataset_name\", \"season\"])\n",
    "    * hv.HLine(-0.5).opts(color=\"red\", line_width=1, line_dash=\"dashed\")\n",
    "    * hv.HLine(0.5).opts(color=\"red\", line_width=1, line_dash=\"dashed\")\n",
    "    * hv.HLine(-0.25).opts(color=\"orange\", line_width=1, line_dash=\"dashed\")\n",
    "    * hv.HLine(0.25).opts(color=\"orange\", line_width=1, line_dash=\"dashed\")\n",
    "    * hv.HLine(0.0).opts(color=\"black\", line_width=1, line_dash=\"dashed\")\n",
    "    * hv.VLine(0.5).opts(color=\"black\", line_width=1, line_dash=\"dashed\")\n",
    ")\n",
    "pn.panel(plt, widget_location=\"top\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1f249e0-1ea0-498a-8b4a-761471ab1154",
   "metadata": {},
   "source": [
    "From the fraction differences, we can see that at least 60% (86%) of the data is within a fractional difference of 0.25 (0.5) of the average, with the GLEAM having the lowest percentage due it having the smallest estimated error variances. Additionally, the quantile distribution is relatively symmetric for each data set with the median close to a zero fractional difference. Therefore, averaging the different combinations should result in a robust mean error estimate.\n",
    "\n",
    "Now that we have confirmed the average is not biased, let's take a look at them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99861940-d810-44cc-aca8-dfbb0f20c463",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt = (\n",
    "    median_tc_est.hvplot(\n",
    "        groupby=[\"dataset_name\", \"season\"],\n",
    "        geo=True,\n",
    "        coastline=True,\n",
    "        clim=(0, 50),\n",
    "        title=\"Median Error Standard Deviation\",\n",
    "    ).opts(frame_width=500)\n",
    "    + std_tc_est.hvplot(\n",
    "        groupby=[\"dataset_name\", \"season\"],\n",
    "        geo=True,\n",
    "        coastline=True,\n",
    "        clim=(0, 15),\n",
    "        title=\"Std of Error Standard Deviation\",\n",
    "    ).opts(frame_width=500)\n",
    "    + median_tc_snr.hvplot(\n",
    "        groupby=[\"dataset_name\", \"season\"],\n",
    "        geo=True,\n",
    "        coastline=True,\n",
    "        clim=(0, 50),\n",
    "        title=\"Median Unbiased SNR\",\n",
    "    ).opts(frame_width=500)\n",
    "    + std_tc_snr.hvplot(\n",
    "        groupby=[\"dataset_name\", \"season\"],\n",
    "        geo=True,\n",
    "        coastline=True,\n",
    "        clim=(0, 25),\n",
    "        title=\"Std of Unbiased SNR\",\n",
    "    ).opts(frame_width=500)\n",
    "    + count_tc_est.hvplot(\n",
    "        groupby=[\"dataset_name\", \"season\"],\n",
    "        geo=True,\n",
    "        coastline=True,\n",
    "        title=\"Number of data points used in calculation\",\n",
    "    ).opts(frame_width=500)\n",
    ")\n",
    "\n",
    "pn.panel(plt.cols(2), widget_location=\"top\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "000c3705-03a5-4a4a-ac3b-f17f2f947fc3",
   "metadata": {},
   "source": [
    "Now that we have some quality error standard deviation estimates, let's discuss how the TC assumptions may affect the results. To remind ourselves of these assumptions, they are:\n",
    "\n",
    "1. The signal and random errors are stationary (i.e., the mean of each is constant with time).\n",
    "2. No cross-correlation of errors (i.e., measurement system errors are independent of each other).\n",
    "3. Error orthogonality (i.e., the measurement system errors are independent of the true value).\n",
    "4. No error autocorrelation (i.e., the error estimates are not correlated with time).\n",
    "5. All data sets are represent exactly the same ET state (i.e., the three data sets have the same spatial resolution and sampling intervals).\n",
    "\n",
    "Of these assumptions, it is possible that each is influencing the result. First, it is possible that our signal and random errors are not stationary. We have performed some stationarity tests (Augmented Dickey Fuller and Kwiatkowski-Phillips-Schmidt-Shin tests) off-hand and found that the signal is stationary for all pixels in one test or the other, indicating the data sets are difference stationary. (This is expected since ET has has a clear seasonal pattern.) Additionally, the error likely has a non-stationary seasonal component. As discussed in [Gruber et al. (2016)](http://dx.doi.org/10.1016/j.jag.2015.09.002), this is not an issue if the data sets all have the same non-stationarity effect, which they likely would if the primary driver is seasonality. However, determining this for the error is difficult and other studies typically avoid the issue by using anomaly space or not discussing it at all. Therefore, gauging its effect is not possible, but likely not a major issue in the error estimates.\n",
    "\n",
    "Of these assumptions, the largest issue in our error variance estimates is likely the inclusion of data sets with cross-correlated errors. Meeting this assumption has been shown to be the most influential in getting correct error variance estimates, more than error orthogonality and error autocorrelation [(Yilmax & Crow 2014)](https://doi.org/10.1175/JHM-D-13-0158.1). From a basic assumption standpoint, as discussed above, we would have expected that SSEBop and GLEAM would be correlated along with ERA5 and NLDAS, and TerraClimate and WBET, as they are generated from similar measurement systems. However, this does not seem to affect the error standard deviation estimates drastically as the estimates are relatively consistent between combinations. However, there is still some significant variation. Therefore, this TC application demonstrates that using basic TC is likely not a reliable measure of **exact** error variances for any arbitrarily chosen data set triplet, as knowing what data sets are truly independent is almost an insurmountable task. However, at the very least, it gives us a lower bound on the error variances as triplets with cross-correlated errors result in underestimated error variances for the correlated pair.\n",
    "\n",
    "Next, error orthogonality and autocorrelation are expected to have a minimal influence on the resulting error estimates. For the orthogonality assumption, the change in magnitude of the estimated errors with season is causing the assumption to not be met. However, its effect is negligible if it is approximately equal for at least two of the three ``independent'' data sets in EC ([Yilmax & Crow 2014](https://doi.org/10.1175/JHM-D-13-0158.1); [Gruber et al. 2016)](http://dx.doi.org/10.1016/j.jag.2015.09.002). As the same seasonal variation in the errors is present in all data sets, we concluded that the error orthogonality assumption, while not met for our data sets, can be expected to have insignificant impact on the error estimates. As for autocorrelation in the errors, the change in errors with season is likely also causing this assumption to not be met. However, it is expected that a violation of this assumption will not affect the accuracy of the error estimates, but it will increase the variation (i.e., precision) of the error estimates [(Zwieback et al. 2012)](https://doi.org/10.5194/npg-19-69-2012).\n",
    "\n",
    "Finally, representativeness can bias one or two of the error variances in the triplet if the spatial or physical representativeness is highly different between the data sets [(Gruber et al. 2016)](http://dx.doi.org/10.1016/j.jag.2015.09.002). In our case in spatial terms, the data sets all originally had their own native resolution, which we degraded to match the GLEAM resolution. Therefore, it is possible that data sets with high native resolution may be penalized for not resolving coarse-scale features in the data sets with lower resolution. Of the six data sets, SSEBop and WBET were similar in resolution, with TerraClimate close to their resolution as well. ERA5 and NLDAS are both also similar in resolution, with them being almost double GLEAM and 10x lower than SSEBop and WBET. However, by conservative regridding, we are aggregating the higher resolution data sets to a common lower resolution such that coarse-scale features outside of a collocated point in the higher resolution data are now present by aggregating multiple points during regridding. Therefore, the effects should be minimal, but if they are not completely suppressed in the data after regridding, it will only result in overestimated error variances. Additionally, the physical representativeness is likely not a large concern either. While each uses different methods to estimate ET, they are attempting to estimate the same physical process. So, they should have similar physical representativeness.\n",
    "\n",
    "| Data Set | SSEBop | GLEAM v3b | ERA5 | NLDAS | TerraClimate | WBET |\n",
    "| ------  |  ----  | -----     | ---- | ----- | ----         | ---- |\n",
    "| Resolution | 0.01 deg (1 km) | 0.25 deg (22.5 km) | 0.1 deg (9 km) | 0.125 deg (11.25 km) | 0.04166 deg (3.75 km) | 0.009 deg (800 m) |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9ad53c9-6a40-4f62-8390-156172c5c30f",
   "metadata": {},
   "source": [
    "Finally, let's look that how the error variance estimations vary between each combination in an aggregate form. We will do this by plotting the whole map of ET estimates as a histogram (i.e., each pixel becomes a count in the histogram). If we find that the error variances are similar between each combination, it may be that cross-correlation is not prominent or it may help us see certain combinations that result in larger error variances, which indicate strong cross-correlations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "047989db-eadd-4e90-8d55-bc2e98dc6569",
   "metadata": {},
   "outputs": [],
   "source": [
    "def histogram_plts(dataset_name=\"SSEBop\", season=\"All\"):\n",
    "    da_pair = []\n",
    "    for i in range(10):\n",
    "        da_pair.append(\n",
    "            tc_est_by_dataset.error.sel(\n",
    "                dataset_name=dataset_name, est_idx=i, season=season\n",
    "            )\n",
    "        )\n",
    "        da_pair[i].name = tc_est_by_dataset.est_pair.sel(\n",
    "            dataset_name=dataset_name, est_idx=i\n",
    "        ).data.item()\n",
    "    da_mean = mean_tc_est.sel(dataset_name=dataset_name, season=season)\n",
    "    da_median = median_tc_est.sel(dataset_name=dataset_name, season=season)\n",
    "    da_mean.name = \"Mean\"\n",
    "    da_median.name = \"Median\"\n",
    "\n",
    "    plt = (\n",
    "        da_pair[0]\n",
    "        .hvplot.hist(\n",
    "            bins=50,\n",
    "            bin_range=(0, 50),\n",
    "            xlabel=\"Error (mm.month-1)\",\n",
    "            ylabel=\"Counts\",\n",
    "            title=\"TC Error Standard Deviation Distribution of \"\n",
    "            + dataset_name\n",
    "            + \" Data Set\",\n",
    "            alpha=0.7,\n",
    "            normed=True,\n",
    "        )\n",
    "        .opts(height=400)\n",
    "    )\n",
    "    for i in range(1, 10):\n",
    "        plt *= da_pair[i].hvplot.hist(\n",
    "            bins=50,\n",
    "            bin_range=(0, 50),\n",
    "            alpha=0.7,\n",
    "            normed=True,\n",
    "            title=\"TC Error Standard Deviation Distribution of \"\n",
    "            + dataset_name\n",
    "            + \" Data Set\",\n",
    "        )\n",
    "    plt *= da_mean.hvplot.hist(\n",
    "        bins=50,\n",
    "        bin_range=(0, 50),\n",
    "        alpha=0.7,\n",
    "        normed=True,\n",
    "        title=\"TC Error Standard Deviation Distribution of \"\n",
    "        + dataset_name\n",
    "        + \" Data Set\",\n",
    "    )\n",
    "    plt *= da_median.hvplot.hist(\n",
    "        bins=50,\n",
    "        bin_range=(0, 50),\n",
    "        alpha=0.7,\n",
    "        normed=True,\n",
    "        title=\"TC Error Standard Deviation Distribution of \"\n",
    "        + dataset_name\n",
    "        + \" Data Set\",\n",
    "    )\n",
    "\n",
    "    return plt\n",
    "\n",
    "\n",
    "def median_table(dataset_name=\"SSEBop\", season=\"All\"):\n",
    "    da_pair = []\n",
    "    for i in range(10):\n",
    "        da_pair.append(\n",
    "            tc_est_by_dataset.error.sel(\n",
    "                dataset_name=dataset_name, est_idx=i, season=season\n",
    "            )\n",
    "        )\n",
    "    da_mean = mean_tc_est.sel(dataset_name=dataset_name, season=season)\n",
    "    da_median = median_tc_est.sel(dataset_name=dataset_name, season=season)\n",
    "\n",
    "    medians = [da.median().data.item() for da in da_pair + [da_mean, da_median]]\n",
    "    table = hv.Table(\n",
    "        {\n",
    "            \"Independent Pair\": list(\n",
    "                tc_est_by_dataset.est_pair.sel(dataset_name=dataset_name).data\n",
    "            )\n",
    "            + [\"Mean\", \"Median\"],\n",
    "            \"Median\": medians,\n",
    "        },\n",
    "        [\"Independent Pair\", \"Median\"],\n",
    "    ).opts(width=250, height=450)\n",
    "\n",
    "    return table\n",
    "\n",
    "\n",
    "dataset_name_widget = pn.widgets.Select(\n",
    "    name=\"dataset_name\",\n",
    "    value=\"SSEBop\",\n",
    "    options=list(tc_est_by_dataset.dataset_name.values),\n",
    ")\n",
    "season_widget = pn.widgets.Select(\n",
    "    name=\"season\", value=\"All\", options=[\"All\", \"DJF\", \"MAM\", \"JJA\", \"SON\"]\n",
    ")\n",
    "\n",
    "bound_plot = pn.bind(\n",
    "    histogram_plts, dataset_name=dataset_name_widget, season=season_widget\n",
    ")\n",
    "bound_table = pn.bind(\n",
    "    median_table, dataset_name=dataset_name_widget, season=season_widget\n",
    ")\n",
    "\n",
    "pn.Column(dataset_name_widget, season_widget, pn.Row(bound_plot, bound_table))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f18ae410-abda-4f76-a7e7-000292fb42ec",
   "metadata": {},
   "source": [
    "From these results we can see that the ET error estimates are actually quite similar between combinations. We can see some interesting results. Like ERA5 and TerraClimate give larger error variances for other data sets when in the combination together, but lower for themselves. This strongly indicates that these two data sets may actually have correlated errors. We can further test this by performing Extended Collocation and having ERA5 and TerraClimate be correlated (which we test in the [next notebook](3_EC_application.ipynb)).\n",
    "\n",
    "Overall though, the strong overlap between each distribution is a nice indication that the TC method can derive reasonable uncertainty estimates that are relatively consistent regardless of the chosen data set triplets. Therefore, including cross-correlated data sets, while having some effect on the estimates, does not cause multiple factor differences in the error estimates. In other words, as long as it is understood that the TC derived uncertainies are lower limits on the uncertainty, using cross-correlated data sets will give viable results when utilized in a TC analysis."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}