diff --git a/analyse-bios-trendy-output.ipynb b/analyse-bios-trendy-output.ipynb new file mode 100644 index 0000000..c91ef0e --- /dev/null +++ b/analyse-bios-trendy-output.ipynb @@ -0,0 +1,4085 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import math\n", + "import pandas as pd" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Todo: Improve chunking for larger datasets\n", + "## from dask.distributed import Client, LocalCluster\n", + "## client = Client()\n", + "## client" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "bios_output_path = \"/g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_output/\"\n", + "trendy_output_path = \"/g/data/rp23/experiments/2024-04-17_BIOS3-merge/lw5085/BIOS_through_TRENDY\"\n", + "\n", + "landmask_path = \"/g/data/rp23/experiments/2024-04-17_BIOS3-merge/lw5085/landmasks\"\n", + "\n", + "bios_act9_output = \"ACT9/S2_mpi\"\n", + "trendy_act9_output = \"S2-act9/output\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "test_landmask = f\"{landmask_path}/Australia_BIOS_9pts_at_0p05_resolution_landmask.nc\"\n", + "landmask = xr.open_dataset(test_landmask)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "bios_cable_output = f\"{bios_output_path}/{bios_act9_output}/bios_out_cable_1900_2022.nc\"\n", + "trendy_cable_output = f\"{trendy_output_path}/{trendy_act9_output}/bios_out_cable_1901_2022.nc\"\n", + "\n", + "bios_casa_output = f\"{bios_output_path}/{bios_act9_output}/bios_out_casa_1900_2022.nc\"\n", + "trendy_casa_output = f\"{trendy_output_path}/{trendy_act9_output}/bios_out_casa_1901_2022.nc\"" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "def preprocess_dataset(d, filter_vars):\n", + " \"\"\"Renames to longtiude/latitude, and only keeps variables listed in cable/casa\"\"\"\n", + " p_dataset = d.drop_vars(list(filter(lambda x: x not in filter_vars, d.keys())))\n", + " if \"x\" in p_dataset.dims:\n", + " p_dataset = p_dataset.rename(name_dict={\"x\": \"longitude\", \"y\": \"latitude\"})\n", + " return p_dataset\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "def extract_points(d, landmask):\n", + " \"\"\"Filter dataset points, given a landmask\"\"\"\n", + " landmask_stack = landmask.stack(z=('latitude', 'longitude'))\n", + " fil_pts = landmask_stack.where(landmask_stack[\"land\"] != 0, drop=True)\n", + " extractors = fil_pts.z.values\n", + " extractors_lat = [i[0] for i in extractors]\n", + " extractors_lon = [i[1] for i in extractors]\n", + " return d.sel(longitude=extractors_lon, latitude=extractors_lat, method=\"nearest\").drop_duplicates(...)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def differentiate_sec_forest(d):\n", + " \"\"\" \n", + " Differentiate primary and secondary forest\n", + " Assumption: Secondary forest can only exist at index 1\n", + " \"\"\"\n", + " new_iveg = d.sel(patch=1)[\"iveg\"]\n", + " new_iveg = new_iveg.where(lambda x: x != 2, -2)\n", + " d[\"iveg\"].loc[dict(patch=1)] = new_iveg\n", + " return d" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# Only used for vector-based outputs\n", + "def land_vec_to_loc(d):\n", + " temp = d.set_index(land=[\"local_lat\", \"local_lon\"])\n", + " temp = temp.unstack(\"land\")\n", + " if \"latitude\" in d.dims:\n", + " temp = temp.drop_vars([\"latitude\", \"longitude\"])\n", + " temp = temp.rename(name_dict={\"local_lon\": \"longitude\", \"local_lat\": \"latitude\"})\n", + " return temp" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "def setup_additional_vars(is_land_vec=False):\n", + " additional_vars = [\"iveg\"]\n", + " additional_vars += [\"local_lat\", \"local_lon\"] if is_land_vec else []\n", + " return additional_vars\n", + "\n", + "def run_cable_workflow(p, l, filter_vars, is_land_vec=False):\n", + " d = xr.open_dataset(p)\n", + "\n", + " # Keep certain variables in dataset\n", + " additional_vars = setup_additional_vars(is_land_vec)\n", + " d = preprocess_dataset(d, filter_vars + additional_vars)\n", + "\n", + " # Select time range\n", + " d = d.sel(time=slice(\"1901-01-01\", \"2022-01-01\")) # TODO: Change to whatever has the smaller dim \n", + "\n", + " # TRENDY merged outputs - only pick points from landmask\n", + " if not is_land_vec:\n", + " d = extract_points(d, l)\n", + "\n", + " # Secondary Evergreen Broadloaf classification\n", + " d = differentiate_sec_forest(d)\n", + "\n", + " # 1D -> 2D dimensionality\n", + " if is_land_vec:\n", + " d = land_vec_to_loc(d)\n", + "\n", + " return d" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Result specific utilities\n", + "\n", + "veg_type_map = {\n", + " -2 : \"peb\",\n", + " 2 : \"seb\",\n", + " 6 : \"c3\",\n", + " 7 : \"c4\",\n", + "}\n", + "\n", + "def gen_vmin_max(dpv, ndpv):\n", + " \"\"\"Generate min/max values for a variable across 2 datasets\"\"\"\n", + " concacted_dataset = xr.concat([dpv, ndpv], dim=\"latitude\")\n", + " return concacted_dataset.min().values, concacted_dataset.max().values\n", + "\n", + "def segregate_iveg(d, veg_type):\n", + " \"\"\"Values across a specific vegetation type\"\"\"\n", + " return d.where(d[\"iveg\"] == veg_type).max(dim=\"patch\")" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "# Analysis 1: Heatmaps\n", + "cable_vars = [\n", + "\"Rnet\",\n", + "\"Qh\",\n", + "\"Qle\",\n", + "\"Qs\",\n", + "\"Qsb\",\n", + "\"GPP\",\n", + "\"NPP\",\n", + "\"NEE\",\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 1MB\n",
+       "Dimensions:    (latitude: 3, longitude: 3, time: 1452, patch: 3)\n",
+       "Coordinates:\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: patch\n",
+       "Data variables:\n",
+       "    Qle        (time, patch, latitude, longitude) float32 157kB 90.39 ... 113.4\n",
+       "    Qh         (time, patch, latitude, longitude) float32 157kB 94.82 ... 19.66\n",
+       "    Qs         (time, patch, latitude, longitude) float32 157kB 0.0 ... 5.451...\n",
+       "    Qsb        (time, patch, latitude, longitude) float32 157kB 0.0 ... 2.846...\n",
+       "    NEE        (time, patch, latitude, longitude) float32 157kB -0.08288 ... ...\n",
+       "    Rnet       (time, patch, latitude, longitude) float32 157kB 178.3 ... 137.2\n",
+       "    GPP        (time, patch, latitude, longitude) float32 157kB 5.63 ... 13.62\n",
+       "    NPP        (time, patch, latitude, longitude) float32 157kB 3.437 ... 9.048\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "Attributes:\n",
+       "    Production:        2024/05/09 at 11:56:19\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  bios\n",
+       "    Output_averaging:  monthly
" + ], + "text/plain": [ + " Size: 1MB\n", + "Dimensions: (latitude: 3, longitude: 3, time: 1452, patch: 3)\n", + "Coordinates:\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: patch\n", + "Data variables:\n", + " Qle (time, patch, latitude, longitude) float32 157kB 90.39 ... 113.4\n", + " Qh (time, patch, latitude, longitude) float32 157kB 94.82 ... 19.66\n", + " Qs (time, patch, latitude, longitude) float32 157kB 0.0 ... 5.451...\n", + " Qsb (time, patch, latitude, longitude) float32 157kB 0.0 ... 2.846...\n", + " NEE (time, patch, latitude, longitude) float32 157kB -0.08288 ... ...\n", + " Rnet (time, patch, latitude, longitude) float32 157kB 178.3 ... 137.2\n", + " GPP (time, patch, latitude, longitude) float32 157kB 5.63 ... 13.62\n", + " NPP (time, patch, latitude, longitude) float32 157kB 3.437 ... 9.048\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + "Attributes:\n", + " Production: 2024/05/09 at 11:56:19\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: bios\n", + " Output_averaging: monthly" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bios_cable_dataset = run_cable_workflow(bios_cable_output, landmask, cable_vars, is_land_vec=True)\n", + "bios_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 1MB\n",
+       "Dimensions:    (longitude: 3, latitude: 3, patch: 3, time: 1452)\n",
+       "Coordinates:\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: patch\n",
+       "Data variables:\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "    Qle        (time, patch, latitude, longitude) float32 157kB 106.7 ... 79.86\n",
+       "    Qh         (time, patch, latitude, longitude) float32 157kB 69.14 ... 50.75\n",
+       "    Qs         (time, patch, latitude, longitude) float32 157kB 0.0 0.0 ... 0.0\n",
+       "    Qsb        (time, patch, latitude, longitude) float32 157kB 0.0 ... 4.22e-05\n",
+       "    NEE        (time, patch, latitude, longitude) float32 157kB 0.276 ... -0....\n",
+       "    Rnet       (time, patch, latitude, longitude) float32 157kB 174.3 ... 132.2\n",
+       "    GPP        (time, patch, latitude, longitude) float32 157kB 4.761 ... 5.517\n",
+       "    NPP        (time, patch, latitude, longitude) float32 157kB 2.345 ... 3.713\n",
+       "Attributes:\n",
+       "    Production:        2025/04/09 at 17:45:55\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  \n",
+       "    Output_averaging:  monthly\n",
+       "    history:           Tue Apr 15 10:45:20 2025: /g/data/rp23/experiments/202...
" + ], + "text/plain": [ + " Size: 1MB\n", + "Dimensions: (longitude: 3, latitude: 3, patch: 3, time: 1452)\n", + "Coordinates:\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: patch\n", + "Data variables:\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + " Qle (time, patch, latitude, longitude) float32 157kB 106.7 ... 79.86\n", + " Qh (time, patch, latitude, longitude) float32 157kB 69.14 ... 50.75\n", + " Qs (time, patch, latitude, longitude) float32 157kB 0.0 0.0 ... 0.0\n", + " Qsb (time, patch, latitude, longitude) float32 157kB 0.0 ... 4.22e-05\n", + " NEE (time, patch, latitude, longitude) float32 157kB 0.276 ... -0....\n", + " Rnet (time, patch, latitude, longitude) float32 157kB 174.3 ... 132.2\n", + " GPP (time, patch, latitude, longitude) float32 157kB 4.761 ... 5.517\n", + " NPP (time, patch, latitude, longitude) float32 157kB 2.345 ... 3.713\n", + "Attributes:\n", + " Production: 2025/04/09 at 17:45:55\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: \n", + " Output_averaging: monthly\n", + " history: Tue Apr 15 10:45:20 2025: /g/data/rp23/experiments/202..." + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "trendy_cable_dataset = run_cable_workflow(trendy_cable_output, landmask, cable_vars)\n", + "trendy_cable_dataset = trendy_cable_dataset.load()\n", + "trendy_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "casa_vars = [\n", + "\"clabile\",\n", + "\"cplant\",\n", + "\"csoil\",\n", + "\"clitter\",\n", + "]\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/namedarray/core.py:514: UserWarning: Duplicate dimension names present: dimensions {'msoil'} appear more than once in dims=('time', 'msoil', 'msoil', 'land'). We do not yet support duplicate dimension names, but we do allow initial construction of the object. We recommend you rename the dims immediately to become distinct, as most xarray functionality is likely to fail silently if you do not. To rename the dimensions you will need to set the ``.dims`` attribute of each variable, ``e.g. var.dims=('x0', 'x1')``.\n", + " warnings.warn(\n", + "/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/namedarray/core.py:514: UserWarning: Duplicate dimension names present: dimensions {'msoil'} appear more than once in dims=('time', 'msoil', 'msoil', 'land'). We do not yet support duplicate dimension names, but we do allow initial construction of the object. We recommend you rename the dims immediately to become distinct, as most xarray functionality is likely to fail silently if you do not. To rename the dimensions you will need to set the ``.dims`` attribute of each variable, ``e.g. var.dims=('x0', 'x1')``.\n", + " warnings.warn(\n", + "/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/namedarray/core.py:514: UserWarning: Duplicate dimension names present: dimensions {'msoil'} appear more than once in dims=('time', 'msoil', 'msoil', 'land'). We do not yet support duplicate dimension names, but we do allow initial construction of the object. We recommend you rename the dims immediately to become distinct, as most xarray functionality is likely to fail silently if you do not. To rename the dimensions you will need to set the ``.dims`` attribute of each variable, ``e.g. var.dims=('x0', 'x1')``.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 2MB\n",
+       "Dimensions:    (patch: 3, latitude: 3, longitude: 3, time: 1476, mplant: 3,\n",
+       "                mlitter: 3, msoil: 3)\n",
+       "Coordinates:\n",
+       "  * patch      (patch) int64 24B 0 1 2\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * time       (time) datetime64[ns] 12kB 1900-02-01 1900-03-01 ... 2023-01-01\n",
+       "Dimensions without coordinates: mplant, mlitter, msoil\n",
+       "Data variables:\n",
+       "    clabile    (time, patch, latitude, longitude) float32 159kB 0.0 ... 2.216\n",
+       "    cplant     (time, mplant, patch, latitude, longitude) float32 478kB 196.6...\n",
+       "    clitter    (time, mlitter, patch, latitude, longitude) float32 478kB 51.4...\n",
+       "    csoil      (time, msoil, patch, latitude, longitude) float32 478kB 46.55 ...\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "Attributes:\n",
+       "    icycle:     2\n",
+       "    startyear:  1900\n",
+       "    endyear:    2022\n",
+       "    runiden:    bios\n",
+       "    run-type:   cable-casa coupled run
" + ], + "text/plain": [ + " Size: 2MB\n", + "Dimensions: (patch: 3, latitude: 3, longitude: 3, time: 1476, mplant: 3,\n", + " mlitter: 3, msoil: 3)\n", + "Coordinates:\n", + " * patch (patch) int64 24B 0 1 2\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * time (time) datetime64[ns] 12kB 1900-02-01 1900-03-01 ... 2023-01-01\n", + "Dimensions without coordinates: mplant, mlitter, msoil\n", + "Data variables:\n", + " clabile (time, patch, latitude, longitude) float32 159kB 0.0 ... 2.216\n", + " cplant (time, mplant, patch, latitude, longitude) float32 478kB 196.6...\n", + " clitter (time, mlitter, patch, latitude, longitude) float32 478kB 51.4...\n", + " csoil (time, msoil, patch, latitude, longitude) float32 478kB 46.55 ...\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + "Attributes:\n", + " icycle: 2\n", + " startyear: 1900\n", + " endyear: 2022\n", + " runiden: bios\n", + " run-type: cable-casa coupled run" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# TODO: Work on after CABLE outputs are clarified\n", + "def set_casa_index(d):\n", + " pass\n", + "def run_casa_workflow(p, l, filter_vars, iveg, is_land_vec=False):\n", + " d = xr.open_dataset(p)\n", + " d[\"patch\"] = d[\"land\"] % 3\n", + " d = d.set_index(land=[\"patch\", \"latitude\", \"longitude\"])\n", + " d = preprocess_dataset(d, casa_vars)\n", + " d = d.unstack(\"land\")\n", + " d[\"iveg\"] = iveg\n", + " return d\n", + "\n", + "# xr.open_dataset(bios_casa_output)\n", + "bios_casa_dataset = run_casa_workflow(bios_casa_output, landmask, casa_vars, bios_cable_dataset[\"iveg\"], is_land_vec=True)\n", + "bios_casa_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "def calc_results(d):\n", + " results = {}\n", + " print(\"Calculating mean\")\n", + " results[\"mean\"] = d.mean(dim=\"time\").compute()\n", + " print(\"Calculating min\")\n", + " results[\"min\"] = d.min(dim=\"time\").compute()\n", + " print(\"Calculating max\")\n", + " results[\"max\"] = d.max(dim=\"time\").compute()\n", + " return results" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Calculating mean\n", + "Calculating min\n", + "Calculating max\n", + "Calculating mean\n", + "Calculating min\n", + "Calculating max\n" + ] + } + ], + "source": [ + "bios_cable_results = calc_results(bios_cable_dataset)\n", + "trendy_cable_results = calc_results(trendy_cable_dataset)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def draw_act9_heatmap(dpv, ax, title, **kwargs):\n", + " sns.heatmap(dpv, annot=True, ax=ax, **kwargs)\n", + " ax.set_yticklabels(dpv.longitude.values)\n", + " ax.set_xticklabels(dpv.latitude.values)\n", + " ax.set_title(title)\n", + " ax.set(xlabel='latitude', ylabel='longitude')\n", + "\n", + "def draw_act_dist(dpv, ndpv, ax):\n", + " diff_vals = abs(dpv - ndpv).values.ravel()\n", + " sns.histplot(diff_vals, ax=ax)\n", + " plt.axvline(diff_vals.mean(), color = 'green', linestyle = 'dashed', linewidth = 3, label='mean')\n", + " ax.set_title(f\"Spatial points differences distribution: mean={diff_vals.mean():.2f}, min={diff_vals.min():.2f}, max={diff_vals.max():.2f}\")\n", + "\n", + "def generate_heatmaps(output_dir, r, n_r, unique_veg_types=veg_type_map.keys()):\n", + " for veg_type in unique_veg_types:\n", + " for analysis in [\"mean\", \"min\", \"max\"]:\n", + " dp = segregate_iveg(r[analysis], veg_type)\n", + " n_dp = segregate_iveg(n_r[analysis], veg_type)\n", + " for v in cable_vars:\n", + " vmin, vmax = gen_vmin_max(dp[v], n_dp[v])\n", + " fig, ax = plt.subplot_mosaic(\"AB;CC\", figsize=(15, 10))\n", + " title = f\"Var: {v} Veg_type: {veg_type_map[veg_type]} Analysis: {analysis} \"\n", + " fig.suptitle(title)\n", + " print(title)\n", + " draw_act9_heatmap(dp[v], ax['A'], \"BIOS\", vmin=vmin, vmax=vmax)\n", + " draw_act9_heatmap(n_dp[v], ax['B'], \"TRENDY\", vmin=vmin, vmax=vmax)\n", + " draw_act_dist(dp[v], n_dp[v], ax['C'])\n", + " plt.savefig(f\"{output_dir}/{v}_{veg_type_map[veg_type]}_{analysis}.png\")\n", + " plt.close(\"all\")" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Var: Rnet Veg_type: peb Analysis: mean \n", + "Var: Qh Veg_type: peb Analysis: mean \n", + "Var: Qle Veg_type: peb Analysis: mean \n", + "Var: Qs Veg_type: peb Analysis: mean \n", + "Var: Qsb Veg_type: peb Analysis: mean \n", + "Var: GPP Veg_type: peb Analysis: mean \n", + "Var: NPP Veg_type: peb Analysis: mean \n", + "Var: NEE Veg_type: peb Analysis: mean \n", + "Var: Rnet Veg_type: peb Analysis: min \n", + "Var: Qh Veg_type: peb Analysis: min \n", + "Var: Qle Veg_type: peb Analysis: min \n", + "Var: Qs Veg_type: peb Analysis: min \n", + "Var: Qsb Veg_type: peb Analysis: min \n", + "Var: GPP Veg_type: peb Analysis: min \n", + "Var: NPP Veg_type: peb Analysis: min \n", + "Var: NEE Veg_type: peb Analysis: min \n", + "Var: Rnet Veg_type: peb Analysis: max \n", + "Var: Qh Veg_type: peb Analysis: max \n", + "Var: Qle Veg_type: peb Analysis: max \n", + "Var: Qs Veg_type: peb Analysis: max \n", + "Var: Qsb Veg_type: peb Analysis: max \n", + "Var: GPP Veg_type: peb Analysis: max \n", + "Var: NPP Veg_type: peb Analysis: max \n", + "Var: NEE Veg_type: peb Analysis: max \n", + "Var: Rnet Veg_type: seb Analysis: mean \n", + "Var: Qh Veg_type: seb Analysis: mean \n", + "Var: Qle Veg_type: seb Analysis: mean \n", + "Var: Qs Veg_type: seb Analysis: mean \n", + "Var: Qsb Veg_type: seb Analysis: mean \n", + "Var: GPP Veg_type: seb Analysis: mean \n", + "Var: NPP Veg_type: seb Analysis: mean \n", + "Var: NEE Veg_type: seb Analysis: mean \n", + "Var: Rnet Veg_type: seb Analysis: min \n", + "Var: Qh Veg_type: seb Analysis: min \n", + "Var: Qle Veg_type: seb Analysis: min \n", + "Var: Qs Veg_type: seb Analysis: min \n", + "Var: Qsb Veg_type: seb Analysis: min \n", + "Var: GPP Veg_type: seb Analysis: min \n", + "Var: NPP Veg_type: seb Analysis: min \n", + "Var: NEE Veg_type: seb Analysis: min \n", + "Var: Rnet Veg_type: seb Analysis: max \n", + "Var: Qh Veg_type: seb Analysis: max \n", + "Var: Qle Veg_type: seb Analysis: max \n", + "Var: Qs Veg_type: seb Analysis: max \n", + "Var: Qsb Veg_type: seb Analysis: max \n", + "Var: GPP Veg_type: seb Analysis: max \n", + "Var: NPP Veg_type: seb Analysis: max \n", + "Var: NEE Veg_type: seb Analysis: max \n", + "Var: Rnet Veg_type: c3 Analysis: mean \n", + "Var: Qh Veg_type: c3 Analysis: mean \n", + "Var: Qle Veg_type: c3 Analysis: mean \n", + "Var: Qs Veg_type: c3 Analysis: mean \n", + "Var: Qsb Veg_type: c3 Analysis: mean \n", + "Var: GPP Veg_type: c3 Analysis: mean \n", + "Var: NPP Veg_type: c3 Analysis: mean \n", + "Var: NEE Veg_type: c3 Analysis: mean \n", + "Var: Rnet Veg_type: c3 Analysis: min \n", + "Var: Qh Veg_type: c3 Analysis: min \n", + "Var: Qle Veg_type: c3 Analysis: min \n", + "Var: Qs Veg_type: c3 Analysis: min \n", + "Var: Qsb Veg_type: c3 Analysis: min \n", + "Var: GPP Veg_type: c3 Analysis: min \n", + "Var: NPP Veg_type: c3 Analysis: min \n", + "Var: NEE Veg_type: c3 Analysis: min \n", + "Var: Rnet Veg_type: c3 Analysis: max \n", + "Var: Qh Veg_type: c3 Analysis: max \n", + "Var: Qle Veg_type: c3 Analysis: max \n", + "Var: Qs Veg_type: c3 Analysis: max \n", + "Var: Qsb Veg_type: c3 Analysis: max \n", + "Var: GPP Veg_type: c3 Analysis: max \n", + "Var: NPP Veg_type: c3 Analysis: max \n", + "Var: NEE Veg_type: c3 Analysis: max \n" + ] + } + ], + "source": [ + "generate_heatmaps(\"output_heatmaps\", bios_cable_results, trendy_cable_results, unique_veg_types=np.unique(bios_cable_dataset[\"iveg\"]))" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "# Analysis 2\n", + "cable_vars_timeseries = [\n", + " \"GPP\",\n", + " \"Qh\",\n", + " \"Qle\",\n", + " \"TVeg\"\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 639kB\n",
+       "Dimensions:    (latitude: 3, longitude: 3, time: 1452, patch: 3)\n",
+       "Coordinates:\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: patch\n",
+       "Data variables:\n",
+       "    Qle        (time, patch, latitude, longitude) float32 157kB 90.39 ... 113.4\n",
+       "    Qh         (time, patch, latitude, longitude) float32 157kB 94.82 ... 19.66\n",
+       "    TVeg       (time, patch, latitude, longitude) float32 157kB 2.957e-05 ......\n",
+       "    GPP        (time, patch, latitude, longitude) float32 157kB 5.63 ... 13.62\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "Attributes:\n",
+       "    Production:        2024/05/09 at 11:56:19\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  bios\n",
+       "    Output_averaging:  monthly
" + ], + "text/plain": [ + " Size: 639kB\n", + "Dimensions: (latitude: 3, longitude: 3, time: 1452, patch: 3)\n", + "Coordinates:\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: patch\n", + "Data variables:\n", + " Qle (time, patch, latitude, longitude) float32 157kB 90.39 ... 113.4\n", + " Qh (time, patch, latitude, longitude) float32 157kB 94.82 ... 19.66\n", + " TVeg (time, patch, latitude, longitude) float32 157kB 2.957e-05 ......\n", + " GPP (time, patch, latitude, longitude) float32 157kB 5.63 ... 13.62\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + "Attributes:\n", + " Production: 2024/05/09 at 11:56:19\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: bios\n", + " Output_averaging: monthly" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bios_cable_dataset = run_cable_workflow(bios_cable_output, landmask, cable_vars_timeseries, is_land_vec=True)\n", + "bios_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "trendy_cable_dataset = run_cable_workflow(trendy_cable_output, landmask, cable_vars_timeseries)\n", + "trendy_cable_dataset = trendy_cable_dataset.load()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def draw_act9_timeseries(dpv, ax, title, **kwargs):\n", + " pass\n", + "\n", + "def generate_timeseries(output_dir, d, n_d, var_list, unique_veg_types=veg_type_map.keys()):\n", + " for veg_type in unique_veg_types:\n", + "\n", + " for lat in d[\"latitude\"]:\n", + " for lon in d[\"longitude\"]:\n", + " dp = segregate_iveg(d.sel(latitude=lat, longitude=lon), veg_type)\n", + " n_dp = segregate_iveg(d.sel(latitude=lat, longitude=lon), veg_type)\n", + "\n", + " # Suplot structure\n", + " fig, ax = plt.subplots(math.ceil(len(var_list) / 2), 2, figsize=(20, 10))\n", + " sup_title = f\"Veg_type: {veg_type_map[veg_type]} lon: {lon.values:.2f} lat: {lat.values:.2f}\"\n", + " fig.suptitle(sup_title)\n", + "\n", + " # Used later for common label\n", + " handles = None\n", + " labels = None\n", + "\n", + " for i, v in enumerate(var_list):\n", + " print(f\"{sup_title} Var: {v}\")\n", + " # NOTE: Difference of 1-2mins after around 1935 but same month for some values \n", + " combined_plot = xr.concat([dp[v], n_dp[v]], dim=pd.Index([\"BIOS\", \"TRENDY\"]), join=\"override\")\n", + " # Plot on the correct axis\n", + " cur_ax = ax[i // 2, i % 2]\n", + " combined_plot.plot.line(x=\"time\", ax=cur_ax, alpha=0.8)\n", + " # Remove Title\n", + " cur_ax.set_title(\"\")\n", + " # Remove legend (common for all)\n", + " handles, labels = cur_ax.get_legend_handles_labels()\n", + " # cur_ax.get_legend().remove()\n", + "\n", + " fig.legend(handles, labels, loc='upper center')\n", + " plt.tight_layout()\n", + " plt.savefig(f\"{output_dir}/lat-{lat.values:.2f}_lon-{lon.values:.2f}_{veg_type_map[veg_type]}.png\")\n", + " plt.close(\"all\")" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Veg_type: seb lon: 149.4499969482422 lat: -35.5 Var: GPP\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.5 Var: Qh\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.5 Var: Qle\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.5 Var: TVeg\n", + "Veg_type: seb lon: 149.5 lat: -35.5 Var: GPP\n", + "Veg_type: seb lon: 149.5 lat: -35.5 Var: Qh\n", + "Veg_type: seb lon: 149.5 lat: -35.5 Var: Qle\n", + "Veg_type: seb lon: 149.5 lat: -35.5 Var: TVeg\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.5 Var: GPP\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.5 Var: Qh\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.5 Var: Qle\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.5 Var: TVeg\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.45000076293945 Var: GPP\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.45000076293945 Var: Qh\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.45000076293945 Var: Qle\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.45000076293945 Var: TVeg\n", + "Veg_type: seb lon: 149.5 lat: -35.45000076293945 Var: GPP\n", + "Veg_type: seb lon: 149.5 lat: -35.45000076293945 Var: Qh\n", + "Veg_type: seb lon: 149.5 lat: -35.45000076293945 Var: Qle\n", + "Veg_type: seb lon: 149.5 lat: -35.45000076293945 Var: TVeg\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.45000076293945 Var: GPP\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.45000076293945 Var: Qh\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.45000076293945 Var: Qle\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.45000076293945 Var: TVeg\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.400001525878906 Var: GPP\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.400001525878906 Var: Qh\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.400001525878906 Var: Qle\n", + "Veg_type: seb lon: 149.4499969482422 lat: -35.400001525878906 Var: TVeg\n", + "Veg_type: seb lon: 149.5 lat: -35.400001525878906 Var: GPP\n", + "Veg_type: seb lon: 149.5 lat: -35.400001525878906 Var: Qh\n", + "Veg_type: seb lon: 149.5 lat: -35.400001525878906 Var: Qle\n", + "Veg_type: seb lon: 149.5 lat: -35.400001525878906 Var: TVeg\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.400001525878906 Var: GPP\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.400001525878906 Var: Qh\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.400001525878906 Var: Qle\n", + "Veg_type: seb lon: 149.5500030517578 lat: -35.400001525878906 Var: TVeg\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.5 Var: GPP\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.5 Var: Qh\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.5 Var: Qle\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.5 Var: TVeg\n", + "Veg_type: c3 lon: 149.5 lat: -35.5 Var: GPP\n", + "Veg_type: c3 lon: 149.5 lat: -35.5 Var: Qh\n", + "Veg_type: c3 lon: 149.5 lat: -35.5 Var: Qle\n", + "Veg_type: c3 lon: 149.5 lat: -35.5 Var: TVeg\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.5 Var: GPP\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.5 Var: Qh\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.5 Var: Qle\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.5 Var: TVeg\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.45000076293945 Var: GPP\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.45000076293945 Var: Qh\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.45000076293945 Var: Qle\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.45000076293945 Var: TVeg\n", + "Veg_type: c3 lon: 149.5 lat: -35.45000076293945 Var: GPP\n", + "Veg_type: c3 lon: 149.5 lat: -35.45000076293945 Var: Qh\n", + "Veg_type: c3 lon: 149.5 lat: -35.45000076293945 Var: Qle\n", + "Veg_type: c3 lon: 149.5 lat: -35.45000076293945 Var: TVeg\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.45000076293945 Var: GPP\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.45000076293945 Var: Qh\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.45000076293945 Var: Qle\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.45000076293945 Var: TVeg\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.400001525878906 Var: GPP\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.400001525878906 Var: Qh\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.400001525878906 Var: Qle\n", + "Veg_type: c3 lon: 149.4499969482422 lat: -35.400001525878906 Var: TVeg\n", + "Veg_type: c3 lon: 149.5 lat: -35.400001525878906 Var: GPP\n", + "Veg_type: c3 lon: 149.5 lat: -35.400001525878906 Var: Qh\n", + "Veg_type: c3 lon: 149.5 lat: -35.400001525878906 Var: Qle\n", + "Veg_type: c3 lon: 149.5 lat: -35.400001525878906 Var: TVeg\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.400001525878906 Var: GPP\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.400001525878906 Var: Qh\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.400001525878906 Var: Qle\n", + "Veg_type: c3 lon: 149.5500030517578 lat: -35.400001525878906 Var: TVeg\n" + ] + } + ], + "source": [ + "generate_timeseries(\"output_timeseries\", bios_cable_dataset, trendy_cable_dataset, cable_vars_timeseries, unique_veg_types=[2, 6])" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "cable_vars_moist = [\n", + " \"SoilMoist\",\n", + " \"swilt\",\n", + " \"patchfrac\"\n", + "]\n", + "\n", + "# NOTE: TRENDY runs have aggregated patchfrac for SoilMoist unlike BIOS\n", + "def aggregate_patchfrac(d, var: str):\n", + " return (d[var] * d[\"patchfrac\"]).sum(dim=\"patch\")" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 326kB\n",
+       "Dimensions:    (latitude: 3, longitude: 3, time: 1452, soil: 6, patch: 3)\n",
+       "Coordinates:\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: soil, patch\n",
+       "Data variables:\n",
+       "    SoilMoist  (time, soil, latitude, longitude) float32 314kB 0.1481 ... 0.2009\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "    patchfrac  (patch, latitude, longitude) float32 108B 0.4767 0.4767 ... 0.0\n",
+       "    swilt      (latitude, longitude) float32 36B 0.136 0.135 ... 0.135 0.136\n",
+       "Attributes:\n",
+       "    Production:        2024/05/09 at 11:56:19\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  bios\n",
+       "    Output_averaging:  monthly
" + ], + "text/plain": [ + " Size: 326kB\n", + "Dimensions: (latitude: 3, longitude: 3, time: 1452, soil: 6, patch: 3)\n", + "Coordinates:\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: soil, patch\n", + "Data variables:\n", + " SoilMoist (time, soil, latitude, longitude) float32 314kB 0.1481 ... 0.2009\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + " patchfrac (patch, latitude, longitude) float32 108B 0.4767 0.4767 ... 0.0\n", + " swilt (latitude, longitude) float32 36B 0.136 0.135 ... 0.135 0.136\n", + "Attributes:\n", + " Production: 2024/05/09 at 11:56:19\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: bios\n", + " Output_averaging: monthly" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bios_cable_dataset = run_cable_workflow(bios_cable_output, landmask, cable_vars_moist, is_land_vec=True)\n", + "bios_cable_dataset[\"SoilMoist\"] = aggregate_patchfrac(bios_cable_dataset, \"SoilMoist\")\n", + "bios_cable_dataset[\"swilt\"] = aggregate_patchfrac(bios_cable_dataset, \"swilt\")\n", + "bios_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [], + "source": [ + "trendy_cable_dataset = run_cable_workflow(trendy_cable_output, landmask, cable_vars_moist)\n", + "trendy_cable_dataset[\"swilt\"] = aggregate_patchfrac(trendy_cable_dataset, \"swilt\")\n", + "trendy_cable_dataset = trendy_cable_dataset.load()" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "def generate_timeseries_soilmoist(output_dir, d, n_d):\n", + " for lat in d[\"latitude\"]:\n", + " for lon in d[\"longitude\"]:\n", + " # Select point\n", + " dp = d.sel(latitude=lat, longitude=lon)\n", + " n_dp = n_d.sel(latitude=lat, longitude=lon)\n", + " # Create subplot\n", + " fig, ax = plt.subplots(math.ceil(d[\"soil\"].size / 2), 2, figsize=(20, 10))\n", + " sup_title = f\"lon: {lon.values:.2f} lat: {lat.values:.2f}\"\n", + " fig.suptitle(sup_title)\n", + " # Used later for common label\n", + " handles = None\n", + " labels = None\n", + " for i in range(d[\"soil\"].size):\n", + " vmin, vmax = gen_vmin_max(dp[\"SoilMoist\"], n_dp[\"SoilMoist\"])\n", + " print(f\"{sup_title} Soil: {i}\")\n", + " soilmoist_dp = dp[\"SoilMoist\"].sel(soil=i)\n", + " soilmoist_ndp = n_dp[\"SoilMoist\"].sel(soil=i)\n", + " cur_ax = ax[i // 2, i % 2]\n", + " combined_plot = xr.concat([soilmoist_dp, soilmoist_ndp],dim=pd.Index([\"BIOS\", \"TRENDY\"]), join=\"override\")\n", + " # Plot on the relevant axis\n", + " combined_plot.plot.line(x=\"time\", ax=cur_ax, alpha=0.8, ylim=(vmin, vmax), add_legend=False)\n", + " cur_ax.legend(labels=[\"BIOS\", \"TRENDY\"], loc=\"upper right\")\n", + " # Horizontal lines for swilt\n", + " cur_ax.axhline(dp[\"swilt\"], color = 'green', linestyle = 'dashed', linewidth = 3, label='swilt_BIOS', alpha=1)\n", + " cur_ax.axhline(n_dp[\"swilt\"], color = 'red', linestyle = 'dashed', linewidth = 3, label='swilt_TRENDY', alpha=1)\n", + " # Remove Title\n", + " cur_ax.set_title(\"\")\n", + " # Remove legend (common for all)\n", + " handles, labels = cur_ax.get_legend_handles_labels()\n", + " print(handles, labels)\n", + " # cur_ax.get_legend().remove()\n", + " fig.legend(handles, labels, loc='upper right')\n", + " plt.tight_layout()\n", + " plt.savefig(f\"{output_dir}/lat-{lat.values:.2f}_lon-{lon.values:.2f}.png\")\n", + " plt.close(\"all\")" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "lon: 149.45 lat: -35.50 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.50 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.50 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.50 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.50 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.50 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.50 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.50 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.50 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.50 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.50 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.50 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.50 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.50 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.50 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.50 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.50 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.50 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.45 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.45 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.45 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.45 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.45 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.45 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.45 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.45 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.45 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.45 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.45 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.45 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.45 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.45 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.45 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.45 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.45 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.45 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.40 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.40 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.40 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.40 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.40 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.45 lat: -35.40 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.40 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.40 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.40 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.40 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.40 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.50 lat: -35.40 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.40 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.40 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.40 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.40 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.40 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "lon: 149.55 lat: -35.40 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n" + ] + } + ], + "source": [ + "generate_timeseries_soilmoist(\"output_soilmoist\", bios_cable_dataset, trendy_cable_dataset)" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "TimedeltaIndex([ '0 days 00:00:00', '-1 days +23:58:56', '0 days 00:01:04',\n", + " '0 days 00:02:08', '-1 days +23:57:52'],\n", + " dtype='timedelta64[ns]', name='time', freq=None)" + ] + }, + "execution_count": 32, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "time_diffs = bios_cable_dataset.time.to_index() - trendy_cable_dataset.time.to_index()\n", + "time_diffs.unique()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "analysis3", + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/analyse-output.ipynb b/analyse-output.ipynb new file mode 100644 index 0000000..efbd04d --- /dev/null +++ b/analyse-output.ipynb @@ -0,0 +1,5327 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import math\n", + "import pandas as pd" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "# Todo: Improve chunking for larger datasets\n", + "## from dask.distributed import Client, LocalCluster\n", + "## client = Client()\n", + "## client" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "bios_output_path = \"/g/data/rp23/experiments/2024-04-17_BIOS3-merge/BIOS3_output/\"\n", + "trendy_output_path = \"/g/data/rp23/experiments/2024-04-17_BIOS3-merge/lw5085/BIOS_through_TRENDY\"\n", + "\n", + "landmask_path = \"/g/data/rp23/experiments/2024-04-17_BIOS3-merge/lw5085/landmasks\"\n", + "\n", + "bios_act9_output = \"ACT9/S2_mpi\"\n", + "trendy_act9_output = \"S2-act9/output\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "test_landmask = f\"{landmask_path}/Australia_BIOS_9pts_at_0p05_resolution_landmask.nc\"\n", + "landmask = xr.open_dataset(test_landmask)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "bios_cable_output = f\"{bios_output_path}/{bios_act9_output}/bios_out_cable_1900_2022.nc\"\n", + "trendy_cable_output = f\"{trendy_output_path}/{trendy_act9_output}/bios_out_cable_1901_2022.nc\"\n", + "\n", + "bios_casa_output = f\"{bios_output_path}/{bios_act9_output}/bios_out_casa_1900_2022.nc\"\n", + "trendy_casa_output = f\"{trendy_output_path}/{trendy_act9_output}/bios_out_casa_1901_2022.nc\"" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [], + "source": [ + "def preprocess_dataset(d, filter_vars):\n", + " \"\"\"Renames to longtiude/latitude, and only keeps variables listed in cable/casa\"\"\"\n", + " p_dataset = d.drop_vars(list(filter(lambda x: x not in filter_vars, d.keys())))\n", + " if \"x\" in p_dataset.dims:\n", + " p_dataset = p_dataset.rename(name_dict={\"x\": \"longitude\", \"y\": \"latitude\"})\n", + " return p_dataset\n" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "def extract_points(d, landmask):\n", + " \"\"\"Filter dataset points, given a landmask\"\"\"\n", + " landmask_stack = landmask.stack(z=('latitude', 'longitude'))\n", + " fil_pts = landmask_stack.where(landmask_stack[\"land\"] != 0, drop=True)\n", + " extractors = fil_pts.z.values\n", + " extractors_lat = [i[0] for i in extractors]\n", + " extractors_lon = [i[1] for i in extractors]\n", + " return d.sel(longitude=extractors_lon, latitude=extractors_lat, method=\"nearest\").drop_duplicates(...)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [], + "source": [ + "def differentiate_sec_forest(d):\n", + " \"\"\" \n", + " Differentiate primary and secondary forest\n", + " Assumption: Secondary forest can only exist at index 1\n", + " \"\"\"\n", + " new_iveg = d.sel(patch=1)[\"iveg\"]\n", + " new_iveg = new_iveg.where(lambda x: x != 2, -2)\n", + " d[\"iveg\"].loc[dict(patch=1)] = new_iveg\n", + " return d" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "# Only used for vector-based outputs\n", + "def land_vec_to_loc(d):\n", + " temp = d.set_index(land=[\"local_lat\", \"local_lon\"])\n", + " temp = temp.unstack(\"land\")\n", + " if \"latitude\" in d.dims:\n", + " temp = temp.drop_vars([\"latitude\", \"longitude\"])\n", + " temp = temp.rename(name_dict={\"local_lon\": \"longitude\", \"local_lat\": \"latitude\"})\n", + " return temp" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [], + "source": [ + "def setup_additional_vars(is_land_vec=False):\n", + " additional_vars = [\"iveg\"]\n", + " additional_vars += [\"local_lat\", \"local_lon\"] if is_land_vec else []\n", + " return additional_vars\n", + "\n", + "def run_cable_workflow(p, l, filter_vars, is_land_vec=False):\n", + " d = xr.open_dataset(p)\n", + "\n", + " # Keep certain variables in dataset\n", + " additional_vars = setup_additional_vars(is_land_vec)\n", + " d = preprocess_dataset(d, filter_vars + additional_vars)\n", + "\n", + " # Select time range\n", + " d = d.sel(time=slice(\"1901-01-01\", \"2022-01-01\")) # TODO: Change to whatever has the smaller dim \n", + "\n", + " # TRENDY merged outputs - only pick points from landmask\n", + " if not is_land_vec:\n", + " d = extract_points(d, l)\n", + "\n", + " # Secondary Evergreen Broadloaf classification\n", + " d = differentiate_sec_forest(d)\n", + "\n", + " # 1D -> 2D dimensionality\n", + " if is_land_vec:\n", + " d = land_vec_to_loc(d)\n", + "\n", + " return d" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [], + "source": [ + "# Result specific utilities\n", + "\n", + "veg_type_map = {\n", + " -2 : \"peb\",\n", + " 2 : \"seb\",\n", + " 6 : \"c3\",\n", + " 7 : \"c4\",\n", + "}\n", + "\n", + "def gen_vmin_max(dpv, ndpv):\n", + " \"\"\"Generate min/max values for a variable across latitude\"\"\"\n", + " concacted_dataset = xr.concat([dpv, ndpv], dim=\"latitude\")\n", + " return concacted_dataset.min().values, concacted_dataset.max().values\n", + "\n", + "def segregate_iveg(d, veg_type):\n", + " \"\"\"Values across a specific vegetation type\"\"\"\n", + " return d.where(d[\"iveg\"] == veg_type).max(dim=\"patch\")" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [], + "source": [ + "# Analysis 1: Heatmaps\n", + "cable_vars = [\n", + "\"Rnet\",\n", + "\"Qh\",\n", + "\"Qle\",\n", + "\"Qs\",\n", + "\"Qsb\",\n", + "\"GPP\",\n", + "\"NPP\",\n", + "\"NEE\",\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 1MB\n",
+       "Dimensions:    (latitude: 3, longitude: 3, time: 1452, patch: 3)\n",
+       "Coordinates:\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: patch\n",
+       "Data variables:\n",
+       "    Qle        (time, patch, latitude, longitude) float32 157kB 90.39 ... 113.4\n",
+       "    Qh         (time, patch, latitude, longitude) float32 157kB 94.82 ... 19.66\n",
+       "    Qs         (time, patch, latitude, longitude) float32 157kB 0.0 ... 5.451...\n",
+       "    Qsb        (time, patch, latitude, longitude) float32 157kB 0.0 ... 2.846...\n",
+       "    NEE        (time, patch, latitude, longitude) float32 157kB -0.08288 ... ...\n",
+       "    Rnet       (time, patch, latitude, longitude) float32 157kB 178.3 ... 137.2\n",
+       "    GPP        (time, patch, latitude, longitude) float32 157kB 5.63 ... 13.62\n",
+       "    NPP        (time, patch, latitude, longitude) float32 157kB 3.437 ... 9.048\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "Attributes:\n",
+       "    Production:        2024/05/09 at 11:56:19\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  bios\n",
+       "    Output_averaging:  monthly
" + ], + "text/plain": [ + " Size: 1MB\n", + "Dimensions: (latitude: 3, longitude: 3, time: 1452, patch: 3)\n", + "Coordinates:\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: patch\n", + "Data variables:\n", + " Qle (time, patch, latitude, longitude) float32 157kB 90.39 ... 113.4\n", + " Qh (time, patch, latitude, longitude) float32 157kB 94.82 ... 19.66\n", + " Qs (time, patch, latitude, longitude) float32 157kB 0.0 ... 5.451...\n", + " Qsb (time, patch, latitude, longitude) float32 157kB 0.0 ... 2.846...\n", + " NEE (time, patch, latitude, longitude) float32 157kB -0.08288 ... ...\n", + " Rnet (time, patch, latitude, longitude) float32 157kB 178.3 ... 137.2\n", + " GPP (time, patch, latitude, longitude) float32 157kB 5.63 ... 13.62\n", + " NPP (time, patch, latitude, longitude) float32 157kB 3.437 ... 9.048\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + "Attributes:\n", + " Production: 2024/05/09 at 11:56:19\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: bios\n", + " Output_averaging: monthly" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bios_cable_dataset = run_cable_workflow(bios_cable_output, landmask, cable_vars, is_land_vec=True)\n", + "bios_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 1MB\n",
+       "Dimensions:    (longitude: 3, latitude: 3, patch: 3, time: 1452)\n",
+       "Coordinates:\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: patch\n",
+       "Data variables:\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "    Qle        (time, patch, latitude, longitude) float32 157kB 106.7 ... 79.86\n",
+       "    Qh         (time, patch, latitude, longitude) float32 157kB 69.14 ... 50.75\n",
+       "    Qs         (time, patch, latitude, longitude) float32 157kB 0.0 0.0 ... 0.0\n",
+       "    Qsb        (time, patch, latitude, longitude) float32 157kB 0.0 ... 4.22e-05\n",
+       "    NEE        (time, patch, latitude, longitude) float32 157kB 0.276 ... -0....\n",
+       "    Rnet       (time, patch, latitude, longitude) float32 157kB 174.3 ... 132.2\n",
+       "    GPP        (time, patch, latitude, longitude) float32 157kB 4.761 ... 5.517\n",
+       "    NPP        (time, patch, latitude, longitude) float32 157kB 2.345 ... 3.713\n",
+       "Attributes:\n",
+       "    Production:        2025/04/09 at 17:45:55\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  \n",
+       "    Output_averaging:  monthly\n",
+       "    history:           Tue Apr 15 10:45:20 2025: /g/data/rp23/experiments/202...
" + ], + "text/plain": [ + " Size: 1MB\n", + "Dimensions: (longitude: 3, latitude: 3, patch: 3, time: 1452)\n", + "Coordinates:\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: patch\n", + "Data variables:\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + " Qle (time, patch, latitude, longitude) float32 157kB 106.7 ... 79.86\n", + " Qh (time, patch, latitude, longitude) float32 157kB 69.14 ... 50.75\n", + " Qs (time, patch, latitude, longitude) float32 157kB 0.0 0.0 ... 0.0\n", + " Qsb (time, patch, latitude, longitude) float32 157kB 0.0 ... 4.22e-05\n", + " NEE (time, patch, latitude, longitude) float32 157kB 0.276 ... -0....\n", + " Rnet (time, patch, latitude, longitude) float32 157kB 174.3 ... 132.2\n", + " GPP (time, patch, latitude, longitude) float32 157kB 4.761 ... 5.517\n", + " NPP (time, patch, latitude, longitude) float32 157kB 2.345 ... 3.713\n", + "Attributes:\n", + " Production: 2025/04/09 at 17:45:55\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: \n", + " Output_averaging: monthly\n", + " history: Tue Apr 15 10:45:20 2025: /g/data/rp23/experiments/202..." + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "trendy_cable_dataset = run_cable_workflow(trendy_cable_output, landmask, cable_vars)\n", + "trendy_cable_dataset = trendy_cable_dataset.load()\n", + "trendy_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "metadata": {}, + "outputs": [], + "source": [ + "casa_vars = [\n", + "\"clabile\",\n", + "\"cplant\",\n", + "\"csoil\",\n", + "\"clitter\",\n", + "]\n" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/namedarray/core.py:514: UserWarning: Duplicate dimension names present: dimensions {'msoil'} appear more than once in dims=('time', 'msoil', 'msoil', 'land'). We do not yet support duplicate dimension names, but we do allow initial construction of the object. We recommend you rename the dims immediately to become distinct, as most xarray functionality is likely to fail silently if you do not. To rename the dimensions you will need to set the ``.dims`` attribute of each variable, ``e.g. var.dims=('x0', 'x1')``.\n", + " warnings.warn(\n", + "/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/namedarray/core.py:514: UserWarning: Duplicate dimension names present: dimensions {'msoil'} appear more than once in dims=('time', 'msoil', 'msoil', 'land'). We do not yet support duplicate dimension names, but we do allow initial construction of the object. We recommend you rename the dims immediately to become distinct, as most xarray functionality is likely to fail silently if you do not. To rename the dimensions you will need to set the ``.dims`` attribute of each variable, ``e.g. var.dims=('x0', 'x1')``.\n", + " warnings.warn(\n", + "/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/namedarray/core.py:514: UserWarning: Duplicate dimension names present: dimensions {'msoil'} appear more than once in dims=('time', 'msoil', 'msoil', 'land'). We do not yet support duplicate dimension names, but we do allow initial construction of the object. We recommend you rename the dims immediately to become distinct, as most xarray functionality is likely to fail silently if you do not. To rename the dimensions you will need to set the ``.dims`` attribute of each variable, ``e.g. var.dims=('x0', 'x1')``.\n", + " warnings.warn(\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 2MB\n",
+       "Dimensions:    (patch: 3, latitude: 3, longitude: 3, time: 1476, mplant: 3,\n",
+       "                mlitter: 3, msoil: 3)\n",
+       "Coordinates:\n",
+       "  * patch      (patch) int64 24B 0 1 2\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * time       (time) datetime64[ns] 12kB 1900-02-01 1900-03-01 ... 2023-01-01\n",
+       "Dimensions without coordinates: mplant, mlitter, msoil\n",
+       "Data variables:\n",
+       "    clabile    (time, patch, latitude, longitude) float32 159kB 0.0 ... 2.216\n",
+       "    cplant     (time, mplant, patch, latitude, longitude) float32 478kB 196.6...\n",
+       "    clitter    (time, mlitter, patch, latitude, longitude) float32 478kB 51.4...\n",
+       "    csoil      (time, msoil, patch, latitude, longitude) float32 478kB 46.55 ...\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "Attributes:\n",
+       "    icycle:     2\n",
+       "    startyear:  1900\n",
+       "    endyear:    2022\n",
+       "    runiden:    bios\n",
+       "    run-type:   cable-casa coupled run
" + ], + "text/plain": [ + " Size: 2MB\n", + "Dimensions: (patch: 3, latitude: 3, longitude: 3, time: 1476, mplant: 3,\n", + " mlitter: 3, msoil: 3)\n", + "Coordinates:\n", + " * patch (patch) int64 24B 0 1 2\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * time (time) datetime64[ns] 12kB 1900-02-01 1900-03-01 ... 2023-01-01\n", + "Dimensions without coordinates: mplant, mlitter, msoil\n", + "Data variables:\n", + " clabile (time, patch, latitude, longitude) float32 159kB 0.0 ... 2.216\n", + " cplant (time, mplant, patch, latitude, longitude) float32 478kB 196.6...\n", + " clitter (time, mlitter, patch, latitude, longitude) float32 478kB 51.4...\n", + " csoil (time, msoil, patch, latitude, longitude) float32 478kB 46.55 ...\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + "Attributes:\n", + " icycle: 2\n", + " startyear: 1900\n", + " endyear: 2022\n", + " runiden: bios\n", + " run-type: cable-casa coupled run" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# TODO: Work on after CABLE outputs are clarified\n", + "def set_casa_index(d):\n", + " pass\n", + "def run_casa_workflow(p, l, filter_vars, iveg, is_land_vec=False):\n", + " d = xr.open_dataset(p)\n", + " d[\"patch\"] = d[\"land\"] % 3\n", + " d = d.set_index(land=[\"patch\", \"latitude\", \"longitude\"])\n", + " d = preprocess_dataset(d, casa_vars)\n", + " d = d.unstack(\"land\")\n", + " d[\"iveg\"] = iveg\n", + " return d\n", + "\n", + "# xr.open_dataset(bios_casa_output)\n", + "bios_casa_dataset = run_casa_workflow(bios_casa_output, landmask, casa_vars, bios_cable_dataset[\"iveg\"], is_land_vec=True)\n", + "bios_casa_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "metadata": {}, + "outputs": [], + "source": [ + "def calc_results(d):\n", + " results = {}\n", + " print(\"Calculating mean\")\n", + " results[\"mean\"] = d.mean(dim=\"time\").compute()\n", + " print(\"Calculating min\")\n", + " results[\"min\"] = d.min(dim=\"time\").compute()\n", + " print(\"Calculating max\")\n", + " results[\"max\"] = d.max(dim=\"time\").compute()\n", + " return results" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Calculating mean\n", + "Calculating min\n", + "Calculating max\n", + "Calculating mean\n", + "Calculating min\n", + "Calculating max\n" + ] + } + ], + "source": [ + "bios_cable_results = calc_results(bios_cable_dataset)\n", + "trendy_cable_results = calc_results(trendy_cable_dataset)" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "metadata": {}, + "outputs": [], + "source": [ + "def draw_act9_heatmap(dpv, ax, title, **kwargs):\n", + " sns.heatmap(dpv, annot=True, fmt=\".4f\", ax=ax, **kwargs)\n", + " ax.set_yticklabels(dpv.longitude.values)\n", + " ax.set_xticklabels(dpv.latitude.values)\n", + " ax.set_title(title)\n", + " ax.set(xlabel='latitude', ylabel='longitude')\n", + "\n", + "def draw_act_dist(dpv, ndpv, ax):\n", + " diff_vals = abs(dpv - ndpv).values.ravel()\n", + " sns.histplot(diff_vals, ax=ax)\n", + " plt.axvline(diff_vals.mean(), color = 'green', linestyle = 'dashed', linewidth = 3, label='mean')\n", + " ax.set_title(f\"Spatial points differences distribution: mean={diff_vals.mean():.2f}, min={diff_vals.min():.2f}, max={diff_vals.max():.2f}\")\n", + "\n", + "def generate_heatmaps(output_dir, r, n_r, unique_veg_types=veg_type_map.keys()):\n", + " for veg_type in unique_veg_types:\n", + " for analysis in [\"mean\", \"min\", \"max\"]:\n", + " dp = segregate_iveg(r[analysis], veg_type)\n", + " n_dp = segregate_iveg(n_r[analysis], veg_type)\n", + " for v in cable_vars:\n", + " vmin, vmax = gen_vmin_max(dp[v], n_dp[v])\n", + " fig, ax = plt.subplot_mosaic(\"AB;CC\", figsize=(15, 10))\n", + " title = f\"Var: {v} Veg_type: {veg_type_map[veg_type]} Analysis: {analysis} \"\n", + " fig.suptitle(title)\n", + " print(title)\n", + " draw_act9_heatmap(dp[v], ax['A'], \"BIOS\", vmin=vmin, vmax=vmax)\n", + " draw_act9_heatmap(n_dp[v], ax['B'], \"TRENDY\", vmin=vmin, vmax=vmax)\n", + " draw_act_dist(dp[v], n_dp[v], ax['C'])\n", + " plt.savefig(f\"{output_dir}/{v}_{veg_type_map[veg_type]}_{analysis}.png\")\n", + " plt.close(\"all\")" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Var: Rnet Veg_type: peb Analysis: mean \n", + "Var: Qh Veg_type: peb Analysis: mean \n", + "Var: Qle Veg_type: peb Analysis: mean \n", + "Var: Qs Veg_type: peb Analysis: mean \n", + "Var: Qsb Veg_type: peb Analysis: mean \n", + "Var: GPP Veg_type: peb Analysis: mean \n", + "Var: NPP Veg_type: peb Analysis: mean \n", + "Var: NEE Veg_type: peb Analysis: mean \n", + "Var: Rnet Veg_type: peb Analysis: min \n", + "Var: Qh Veg_type: peb Analysis: min \n", + "Var: Qle Veg_type: peb Analysis: min \n", + "Var: Qs Veg_type: peb Analysis: min \n", + "Var: Qsb Veg_type: peb Analysis: min \n", + "Var: GPP Veg_type: peb Analysis: min \n", + "Var: NPP Veg_type: peb Analysis: min \n", + "Var: NEE Veg_type: peb Analysis: min \n", + "Var: Rnet Veg_type: peb Analysis: max \n", + "Var: Qh Veg_type: peb Analysis: max \n", + "Var: Qle Veg_type: peb Analysis: max \n", + "Var: Qs Veg_type: peb Analysis: max \n", + "Var: Qsb Veg_type: peb Analysis: max \n", + "Var: GPP Veg_type: peb Analysis: max \n", + "Var: NPP Veg_type: peb Analysis: max \n", + "Var: NEE Veg_type: peb Analysis: max \n", + "Var: Rnet Veg_type: seb Analysis: mean \n", + "Var: Qh Veg_type: seb Analysis: mean \n", + "Var: Qle Veg_type: seb Analysis: mean \n", + "Var: Qs Veg_type: seb Analysis: mean \n", + "Var: Qsb Veg_type: seb Analysis: mean \n", + "Var: GPP Veg_type: seb Analysis: mean \n", + "Var: NPP Veg_type: seb Analysis: mean \n", + "Var: NEE Veg_type: seb Analysis: mean \n", + "Var: Rnet Veg_type: seb Analysis: min \n", + "Var: Qh Veg_type: seb Analysis: min \n", + "Var: Qle Veg_type: seb Analysis: min \n", + "Var: Qs Veg_type: seb Analysis: min \n", + "Var: Qsb Veg_type: seb Analysis: min \n", + "Var: GPP Veg_type: seb Analysis: min \n", + "Var: NPP Veg_type: seb Analysis: min \n", + "Var: NEE Veg_type: seb Analysis: min \n", + "Var: Rnet Veg_type: seb Analysis: max \n", + "Var: Qh Veg_type: seb Analysis: max \n", + "Var: Qle Veg_type: seb Analysis: max \n", + "Var: Qs Veg_type: seb Analysis: max \n", + "Var: Qsb Veg_type: seb Analysis: max \n", + "Var: GPP Veg_type: seb Analysis: max \n", + "Var: NPP Veg_type: seb Analysis: max \n", + "Var: NEE Veg_type: seb Analysis: max \n", + "Var: Rnet Veg_type: c3 Analysis: mean \n", + "Var: Qh Veg_type: c3 Analysis: mean \n", + "Var: Qle Veg_type: c3 Analysis: mean \n", + "Var: Qs Veg_type: c3 Analysis: mean \n", + "Var: Qsb Veg_type: c3 Analysis: mean \n", + "Var: GPP Veg_type: c3 Analysis: mean \n", + "Var: NPP Veg_type: c3 Analysis: mean \n", + "Var: NEE Veg_type: c3 Analysis: mean \n", + "Var: Rnet Veg_type: c3 Analysis: min \n", + "Var: Qh Veg_type: c3 Analysis: min \n", + "Var: Qle Veg_type: c3 Analysis: min \n", + "Var: Qs Veg_type: c3 Analysis: min \n", + "Var: Qsb Veg_type: c3 Analysis: min \n", + "Var: GPP Veg_type: c3 Analysis: min \n", + "Var: NPP Veg_type: c3 Analysis: min \n", + "Var: NEE Veg_type: c3 Analysis: min \n", + "Var: Rnet Veg_type: c3 Analysis: max \n", + "Var: Qh Veg_type: c3 Analysis: max \n", + "Var: Qle Veg_type: c3 Analysis: max \n", + "Var: Qs Veg_type: c3 Analysis: max \n", + "Var: Qsb Veg_type: c3 Analysis: max \n", + "Var: GPP Veg_type: c3 Analysis: max \n", + "Var: NPP Veg_type: c3 Analysis: max \n", + "Var: NEE Veg_type: c3 Analysis: max \n" + ] + } + ], + "source": [ + "generate_heatmaps(\"output_heatmaps\", bios_cable_results, trendy_cable_results, unique_veg_types=np.unique(bios_cable_dataset[\"iveg\"]))" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "# Analysis 2\n", + "cable_vars_timeseries = [\n", + " \"GPP\",\n", + " \"Qh\",\n", + " \"Qle\",\n", + " \"TVeg\"\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 639kB\n",
+       "Dimensions:    (latitude: 3, longitude: 3, time: 1452, patch: 3)\n",
+       "Coordinates:\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: patch\n",
+       "Data variables:\n",
+       "    Qle        (time, patch, latitude, longitude) float32 157kB 90.39 ... 113.4\n",
+       "    Qh         (time, patch, latitude, longitude) float32 157kB 94.82 ... 19.66\n",
+       "    TVeg       (time, patch, latitude, longitude) float32 157kB 2.957e-05 ......\n",
+       "    GPP        (time, patch, latitude, longitude) float32 157kB 5.63 ... 13.62\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "Attributes:\n",
+       "    Production:        2024/05/09 at 11:56:19\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  bios\n",
+       "    Output_averaging:  monthly
" + ], + "text/plain": [ + " Size: 639kB\n", + "Dimensions: (latitude: 3, longitude: 3, time: 1452, patch: 3)\n", + "Coordinates:\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: patch\n", + "Data variables:\n", + " Qle (time, patch, latitude, longitude) float32 157kB 90.39 ... 113.4\n", + " Qh (time, patch, latitude, longitude) float32 157kB 94.82 ... 19.66\n", + " TVeg (time, patch, latitude, longitude) float32 157kB 2.957e-05 ......\n", + " GPP (time, patch, latitude, longitude) float32 157kB 5.63 ... 13.62\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + "Attributes:\n", + " Production: 2024/05/09 at 11:56:19\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: bios\n", + " Output_averaging: monthly" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bios_cable_dataset = run_cable_workflow(bios_cable_output, landmask, cable_vars_timeseries, is_land_vec=True)\n", + "bios_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "trendy_cable_dataset = run_cable_workflow(trendy_cable_output, landmask, cable_vars_timeseries)\n", + "trendy_cable_dataset = trendy_cable_dataset.load()" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [], + "source": [ + "def draw_act9_timeseries(dpv, ax, title, **kwargs):\n", + " pass\n", + "\n", + "def generate_timeseries(output_dir, d, n_d, var_list, unique_veg_types=veg_type_map.keys()):\n", + " for veg_type in unique_veg_types:\n", + "\n", + " for lat in d[\"latitude\"]:\n", + " for lon in d[\"longitude\"]:\n", + " dp = segregate_iveg(d.sel(latitude=lat, longitude=lon), veg_type)\n", + " n_dp = segregate_iveg(n_d.sel(latitude=lat, longitude=lon), veg_type)\n", + "\n", + " # Suplot structure\n", + " fig, ax = plt.subplots(math.ceil(len(var_list) / 2), 2, figsize=(20, 10))\n", + " sup_title = f\"Veg_type: {veg_type_map[veg_type]} lon: {lon.values:.2f} lat: {lat.values:.2f}\"\n", + " fig.suptitle(sup_title)\n", + "\n", + " # Used later for common label\n", + " handles = None\n", + " labels = None\n", + "\n", + " for i, v in enumerate(var_list):\n", + " print(f\"{sup_title} Var: {v}\")\n", + " # NOTE: Difference of 1-2mins after around 1935 but same month for some values \n", + " combined_plot = xr.concat([dp[v], n_dp[v]], dim=pd.Index([\"BIOS\", \"TRENDY\"]), join=\"override\")\n", + " # Plot on the correct axis\n", + " cur_ax = ax[i // 2, i % 2]\n", + " combined_plot.plot.line(x=\"time\", ax=cur_ax, alpha=0.8)\n", + " # Remove Title\n", + " cur_ax.set_title(\"\")\n", + " # Remove legend (common for all)\n", + " handles, labels = cur_ax.get_legend_handles_labels()\n", + " # cur_ax.get_legend().remove()\n", + "\n", + " fig.legend(handles, labels, loc='upper center')\n", + " plt.tight_layout()\n", + " plt.savefig(f\"{output_dir}/lat-{lat.values:.2f}_lon-{lon.values:.2f}_{veg_type_map[veg_type]}.png\")\n", + " plt.close(\"all\")" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Veg_type: seb lon: 149.45 lat: -35.50 Var: GPP\n", + "Veg_type: seb lon: 149.45 lat: -35.50 Var: Qh\n", + "Veg_type: seb lon: 149.45 lat: -35.50 Var: Qle\n", + "Veg_type: seb lon: 149.45 lat: -35.50 Var: TVeg\n", + "Veg_type: seb lon: 149.50 lat: -35.50 Var: GPP\n", + "Veg_type: seb lon: 149.50 lat: -35.50 Var: Qh\n", + "Veg_type: seb lon: 149.50 lat: -35.50 Var: Qle\n", + "Veg_type: seb lon: 149.50 lat: -35.50 Var: TVeg\n", + "Veg_type: seb lon: 149.55 lat: -35.50 Var: GPP\n", + "Veg_type: seb lon: 149.55 lat: -35.50 Var: Qh\n", + "Veg_type: seb lon: 149.55 lat: -35.50 Var: Qle\n", + "Veg_type: seb lon: 149.55 lat: -35.50 Var: TVeg\n", + "Veg_type: seb lon: 149.45 lat: -35.45 Var: GPP\n", + "Veg_type: seb lon: 149.45 lat: -35.45 Var: Qh\n", + "Veg_type: seb lon: 149.45 lat: -35.45 Var: Qle\n", + "Veg_type: seb lon: 149.45 lat: -35.45 Var: TVeg\n", + "Veg_type: seb lon: 149.50 lat: -35.45 Var: GPP\n", + "Veg_type: seb lon: 149.50 lat: -35.45 Var: Qh\n", + "Veg_type: seb lon: 149.50 lat: -35.45 Var: Qle\n", + "Veg_type: seb lon: 149.50 lat: -35.45 Var: TVeg\n", + "Veg_type: seb lon: 149.55 lat: -35.45 Var: GPP\n", + "Veg_type: seb lon: 149.55 lat: -35.45 Var: Qh\n", + "Veg_type: seb lon: 149.55 lat: -35.45 Var: Qle\n", + "Veg_type: seb lon: 149.55 lat: -35.45 Var: TVeg\n", + "Veg_type: seb lon: 149.45 lat: -35.40 Var: GPP\n", + "Veg_type: seb lon: 149.45 lat: -35.40 Var: Qh\n", + "Veg_type: seb lon: 149.45 lat: -35.40 Var: Qle\n", + "Veg_type: seb lon: 149.45 lat: -35.40 Var: TVeg\n", + "Veg_type: seb lon: 149.50 lat: -35.40 Var: GPP\n", + "Veg_type: seb lon: 149.50 lat: -35.40 Var: Qh\n", + "Veg_type: seb lon: 149.50 lat: -35.40 Var: Qle\n", + "Veg_type: seb lon: 149.50 lat: -35.40 Var: TVeg\n", + "Veg_type: seb lon: 149.55 lat: -35.40 Var: GPP\n", + "Veg_type: seb lon: 149.55 lat: -35.40 Var: Qh\n", + "Veg_type: seb lon: 149.55 lat: -35.40 Var: Qle\n", + "Veg_type: seb lon: 149.55 lat: -35.40 Var: TVeg\n", + "Veg_type: c3 lon: 149.45 lat: -35.50 Var: GPP\n", + "Veg_type: c3 lon: 149.45 lat: -35.50 Var: Qh\n", + "Veg_type: c3 lon: 149.45 lat: -35.50 Var: Qle\n", + "Veg_type: c3 lon: 149.45 lat: -35.50 Var: TVeg\n", + "Veg_type: c3 lon: 149.50 lat: -35.50 Var: GPP\n", + "Veg_type: c3 lon: 149.50 lat: -35.50 Var: Qh\n", + "Veg_type: c3 lon: 149.50 lat: -35.50 Var: Qle\n", + "Veg_type: c3 lon: 149.50 lat: -35.50 Var: TVeg\n", + "Veg_type: c3 lon: 149.55 lat: -35.50 Var: GPP\n", + "Veg_type: c3 lon: 149.55 lat: -35.50 Var: Qh\n", + "Veg_type: c3 lon: 149.55 lat: -35.50 Var: Qle\n", + "Veg_type: c3 lon: 149.55 lat: -35.50 Var: TVeg\n", + "Veg_type: c3 lon: 149.45 lat: -35.45 Var: GPP\n", + "Veg_type: c3 lon: 149.45 lat: -35.45 Var: Qh\n", + "Veg_type: c3 lon: 149.45 lat: -35.45 Var: Qle\n", + "Veg_type: c3 lon: 149.45 lat: -35.45 Var: TVeg\n", + "Veg_type: c3 lon: 149.50 lat: -35.45 Var: GPP\n", + "Veg_type: c3 lon: 149.50 lat: -35.45 Var: Qh\n", + "Veg_type: c3 lon: 149.50 lat: -35.45 Var: Qle\n", + "Veg_type: c3 lon: 149.50 lat: -35.45 Var: TVeg\n", + "Veg_type: c3 lon: 149.55 lat: -35.45 Var: GPP\n", + "Veg_type: c3 lon: 149.55 lat: -35.45 Var: Qh\n", + "Veg_type: c3 lon: 149.55 lat: -35.45 Var: Qle\n", + "Veg_type: c3 lon: 149.55 lat: -35.45 Var: TVeg\n", + "Veg_type: c3 lon: 149.45 lat: -35.40 Var: GPP\n", + "Veg_type: c3 lon: 149.45 lat: -35.40 Var: Qh\n", + "Veg_type: c3 lon: 149.45 lat: -35.40 Var: Qle\n", + "Veg_type: c3 lon: 149.45 lat: -35.40 Var: TVeg\n", + "Veg_type: c3 lon: 149.50 lat: -35.40 Var: GPP\n", + "Veg_type: c3 lon: 149.50 lat: -35.40 Var: Qh\n", + "Veg_type: c3 lon: 149.50 lat: -35.40 Var: Qle\n", + "Veg_type: c3 lon: 149.50 lat: -35.40 Var: TVeg\n", + "Veg_type: c3 lon: 149.55 lat: -35.40 Var: GPP\n", + "Veg_type: c3 lon: 149.55 lat: -35.40 Var: Qh\n", + "Veg_type: c3 lon: 149.55 lat: -35.40 Var: Qle\n", + "Veg_type: c3 lon: 149.55 lat: -35.40 Var: TVeg\n" + ] + } + ], + "source": [ + "generate_timeseries(\"output_timeseries\", bios_cable_dataset, trendy_cable_dataset, cable_vars_timeseries, unique_veg_types=[2, 6])" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [], + "source": [ + "cable_vars_moist = [\n", + " \"SoilMoist\",\n", + " \"swilt\",\n", + " \"patchfrac\"\n", + "]\n", + "\n", + "# NOTE: TRENDY runs have aggregated patchfrac for SoilMoist unlike BIOS\n", + "def aggregate_patchfrac(d, var: str):\n", + " return (d[var] * d[\"patchfrac\"]).sum(dim=\"patch\")" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 326kB\n",
+       "Dimensions:    (latitude: 3, longitude: 3, time: 1452, soil: 6, patch: 3)\n",
+       "Coordinates:\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: soil, patch\n",
+       "Data variables:\n",
+       "    SoilMoist  (time, soil, latitude, longitude) float32 314kB 0.1481 ... 0.2009\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "    patchfrac  (patch, latitude, longitude) float32 108B 0.4767 0.4767 ... 0.0\n",
+       "    swilt      (latitude, longitude) float32 36B 0.136 0.135 ... 0.135 0.136\n",
+       "Attributes:\n",
+       "    Production:        2024/05/09 at 11:56:19\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  bios\n",
+       "    Output_averaging:  monthly
" + ], + "text/plain": [ + " Size: 326kB\n", + "Dimensions: (latitude: 3, longitude: 3, time: 1452, soil: 6, patch: 3)\n", + "Coordinates:\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: soil, patch\n", + "Data variables:\n", + " SoilMoist (time, soil, latitude, longitude) float32 314kB 0.1481 ... 0.2009\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + " patchfrac (patch, latitude, longitude) float32 108B 0.4767 0.4767 ... 0.0\n", + " swilt (latitude, longitude) float32 36B 0.136 0.135 ... 0.135 0.136\n", + "Attributes:\n", + " Production: 2024/05/09 at 11:56:19\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: bios\n", + " Output_averaging: monthly" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bios_cable_dataset = run_cable_workflow(bios_cable_output, landmask, cable_vars_moist, is_land_vec=True)\n", + "bios_cable_dataset[\"SoilMoist\"] = aggregate_patchfrac(bios_cable_dataset, \"SoilMoist\")\n", + "bios_cable_dataset[\"swilt\"] = aggregate_patchfrac(bios_cable_dataset, \"swilt\")\n", + "bios_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [], + "source": [ + "trendy_cable_dataset = run_cable_workflow(trendy_cable_output, landmask, cable_vars_moist)\n", + "trendy_cable_dataset[\"swilt\"] = aggregate_patchfrac(trendy_cable_dataset, \"swilt\")\n", + "trendy_cable_dataset = trendy_cable_dataset.load()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "def generate_timeseries_soilmoist(output_dir, d, n_d):\n", + " for lat in d[\"latitude\"]:\n", + " for lon in d[\"longitude\"]:\n", + " # Select point\n", + " dp = d.sel(latitude=lat, longitude=lon)\n", + " n_dp = n_d.sel(latitude=lat, longitude=lon)\n", + " # Create subplot\n", + " fig, ax = plt.subplots(math.ceil(d[\"soil\"].size / 2), 2, figsize=(20, 10))\n", + " sup_title = f\"lon: {lon.values:.2f} lat: {lat.values:.2f}\"\n", + " fig.suptitle(sup_title)\n", + " # Used later for common label\n", + " handles = None\n", + " labels = None\n", + " for i in range(d[\"soil\"].size):\n", + " vmin, vmax = gen_vmin_max(dp[\"SoilMoist\"], n_dp[\"SoilMoist\"])\n", + " print(f\"{sup_title} Soil: {i}\")\n", + " soilmoist_dp = dp[\"SoilMoist\"].sel(soil=i)\n", + " soilmoist_ndp = n_dp[\"SoilMoist\"].sel(soil=i)\n", + " cur_ax = ax[i // 2, i % 2]\n", + " combined_plot = xr.concat([soilmoist_dp, soilmoist_ndp],dim=pd.Index([\"BIOS\", \"TRENDY\"]), join=\"override\")\n", + " # Plot on the relevant axis\n", + " combined_plot.plot.line(x=\"time\", ax=cur_ax, alpha=0.8, ylim=(vmin, vmax), add_legend=False)\n", + " cur_ax.legend(labels=[\"BIOS\", \"TRENDY\"], loc=\"upper right\")\n", + " # Horizontal lines for swilt\n", + " cur_ax.axhline(dp[\"swilt\"], color = 'green', linestyle = 'dashed', linewidth = 3, label='swilt_BIOS', alpha=1)\n", + " cur_ax.axhline(n_dp[\"swilt\"], color = 'red', linestyle = 'dashed', linewidth = 3, label='swilt_TRENDY', alpha=1)\n", + " # Remove Title\n", + " cur_ax.set_title(\"\")\n", + " # Remove legend (common for all)\n", + " handles, labels = cur_ax.get_legend_handles_labels()\n", + " print(handles, labels)\n", + " # cur_ax.get_legend().remove()\n", + " fig.legend(handles, labels, loc='upper right')\n", + " plt.tight_layout()\n", + " plt.savefig(f\"{output_dir}/lat-{lat.values:.2f}_lon-{lon.values:.2f}.png\")\n", + " plt.close(\"all\")" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.136 0.136\n", + "lon: 149.45 lat: -35.50 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.45 lat: -35.50 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.45 lat: -35.50 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.45 lat: -35.50 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.45 lat: -35.50 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.45 lat: -35.50 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.50 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.50 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.50 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.50 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.50 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.50 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.50 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.50 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.50 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.50 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.50 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.50 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.45 lat: -35.45 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.45 lat: -35.45 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.45 lat: -35.45 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.45 lat: -35.45 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.45 lat: -35.45 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.45 lat: -35.45 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.45 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.45 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.45 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.45 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.45 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.45 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.45 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.45 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.45 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.45 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.45 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.45 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.134 0.134\n", + "lon: 149.45 lat: -35.40 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.134 0.134\n", + "lon: 149.45 lat: -35.40 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.134 0.134\n", + "lon: 149.45 lat: -35.40 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.134 0.134\n", + "lon: 149.45 lat: -35.40 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.134 0.134\n", + "lon: 149.45 lat: -35.40 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.134 0.134\n", + "lon: 149.45 lat: -35.40 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.40 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.40 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.40 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.40 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.40 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.135 0.135\n", + "lon: 149.50 lat: -35.40 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.40 Soil: 0\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.40 Soil: 1\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.40 Soil: 2\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.40 Soil: 3\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.40 Soil: 4\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n", + "0.136 0.136\n", + "lon: 149.55 lat: -35.40 Soil: 5\n", + "[, ] ['swilt_BIOS', 'swilt_TRENDY']\n" + ] + } + ], + "source": [ + "generate_timeseries_soilmoist(\"output_soilmoist\", bios_cable_dataset, trendy_cable_dataset)" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "TimedeltaIndex([ '0 days 00:00:00', '-1 days +23:58:56', '0 days 00:01:04',\n", + " '0 days 00:02:08', '-1 days +23:57:52'],\n", + " dtype='timedelta64[ns]', name='time', freq=None)" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "time_diffs = bios_cable_dataset.time.to_index() - trendy_cable_dataset.time.to_index()\n", + "time_diffs.unique()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Analysis 4: Constant var comparison\n", + "cable_vars = [\n", + " \"silt\",\n", + " \"clay\",\n", + " \"css\",\n", + " \"sfc\",\n", + " \"rhosoil\",\n", + " \"bch\",\n", + " \"hyds\",\n", + " \"ssat\",\n", + " \"swilt\",\n", + " \"sucs\",\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 13kB\n",
+       "Dimensions:    (latitude: 3, longitude: 3, time: 1452, patch: 3)\n",
+       "Coordinates:\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: patch\n",
+       "Data variables:\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "    bch        (patch, latitude, longitude) float32 108B 5.1 5.2 5.1 ... 5.2 5.1\n",
+       "    clay       (patch, latitude, longitude) float32 108B 0.15 0.15 ... 0.15 0.15\n",
+       "    silt       (patch, latitude, longitude) float32 108B 0.1 0.1 0.1 ... 0.1 0.1\n",
+       "    ssat       (patch, latitude, longitude) float32 108B 0.474 0.439 ... 0.474\n",
+       "    sfc        (patch, latitude, longitude) float32 108B 0.294 0.282 ... 0.294\n",
+       "    swilt      (patch, latitude, longitude) float32 108B 0.136 0.135 ... 0.136\n",
+       "    hyds       (patch, latitude, longitude) float32 108B 8.33e-05 ... 8.33e-05\n",
+       "    sucs       (patch, latitude, longitude) float32 108B -0.3795 ... -0.3795\n",
+       "    css        (patch, latitude, longitude) float32 108B 303.0 303.0 ... 303.0\n",
+       "    rhosoil    (patch, latitude, longitude) float32 108B 1.3e+03 ... 1.3e+03\n",
+       "Attributes:\n",
+       "    Production:        2024/05/09 at 11:56:19\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  bios\n",
+       "    Output_averaging:  monthly
" + ], + "text/plain": [ + " Size: 13kB\n", + "Dimensions: (latitude: 3, longitude: 3, time: 1452, patch: 3)\n", + "Coordinates:\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: patch\n", + "Data variables:\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + " bch (patch, latitude, longitude) float32 108B 5.1 5.2 5.1 ... 5.2 5.1\n", + " clay (patch, latitude, longitude) float32 108B 0.15 0.15 ... 0.15 0.15\n", + " silt (patch, latitude, longitude) float32 108B 0.1 0.1 0.1 ... 0.1 0.1\n", + " ssat (patch, latitude, longitude) float32 108B 0.474 0.439 ... 0.474\n", + " sfc (patch, latitude, longitude) float32 108B 0.294 0.282 ... 0.294\n", + " swilt (patch, latitude, longitude) float32 108B 0.136 0.135 ... 0.136\n", + " hyds (patch, latitude, longitude) float32 108B 8.33e-05 ... 8.33e-05\n", + " sucs (patch, latitude, longitude) float32 108B -0.3795 ... -0.3795\n", + " css (patch, latitude, longitude) float32 108B 303.0 303.0 ... 303.0\n", + " rhosoil (patch, latitude, longitude) float32 108B 1.3e+03 ... 1.3e+03\n", + "Attributes:\n", + " Production: 2024/05/09 at 11:56:19\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: bios\n", + " Output_averaging: monthly" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bios_cable_dataset = run_cable_workflow(bios_cable_output, landmask, cable_vars, is_land_vec=True)\n", + "bios_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 13kB\n",
+       "Dimensions:    (longitude: 3, latitude: 3, patch: 3, time: 1452)\n",
+       "Coordinates:\n",
+       "  * longitude  (longitude) float32 12B 149.4 149.5 149.6\n",
+       "  * latitude   (latitude) float32 12B -35.5 -35.45 -35.4\n",
+       "  * time       (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n",
+       "Dimensions without coordinates: patch\n",
+       "Data variables:\n",
+       "    iveg       (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n",
+       "    bch        (patch, latitude, longitude) float32 108B 5.1 5.2 5.1 ... 5.2 5.1\n",
+       "    clay       (patch, latitude, longitude) float32 108B 0.15 0.15 ... 0.15 0.15\n",
+       "    silt       (patch, latitude, longitude) float32 108B 0.1 0.1 0.1 ... 0.1 0.1\n",
+       "    ssat       (patch, latitude, longitude) float32 108B 0.474 0.439 ... 0.474\n",
+       "    sfc        (patch, latitude, longitude) float32 108B 0.294 0.282 ... 0.294\n",
+       "    swilt      (patch, latitude, longitude) float32 108B 0.136 0.135 ... 0.136\n",
+       "    hyds       (patch, latitude, longitude) float32 108B 8.33e-05 ... 8.33e-05\n",
+       "    sucs       (patch, latitude, longitude) float32 108B -0.3795 ... -0.3795\n",
+       "    css        (patch, latitude, longitude) float32 108B 303.0 303.0 ... 303.0\n",
+       "    rhosoil    (patch, latitude, longitude) float32 108B 1.3e+03 ... 1.3e+03\n",
+       "Attributes:\n",
+       "    Production:        2025/04/09 at 17:45:55\n",
+       "    Source:            CABLE LSM output file\n",
+       "    CABLE_input_file:  \n",
+       "    Output_averaging:  monthly\n",
+       "    history:           Tue Apr 15 10:45:20 2025: /g/data/rp23/experiments/202...
" + ], + "text/plain": [ + " Size: 13kB\n", + "Dimensions: (longitude: 3, latitude: 3, patch: 3, time: 1452)\n", + "Coordinates:\n", + " * longitude (longitude) float32 12B 149.4 149.5 149.6\n", + " * latitude (latitude) float32 12B -35.5 -35.45 -35.4\n", + " * time (time) datetime64[ns] 12kB 1901-01-16T12:00:00 ... 2021-12-16T...\n", + "Dimensions without coordinates: patch\n", + "Data variables:\n", + " iveg (patch, latitude, longitude) float64 216B 2.0 2.0 2.0 ... 6.0 6.0\n", + " bch (patch, latitude, longitude) float32 108B 5.1 5.2 5.1 ... 5.2 5.1\n", + " clay (patch, latitude, longitude) float32 108B 0.15 0.15 ... 0.15 0.15\n", + " silt (patch, latitude, longitude) float32 108B 0.1 0.1 0.1 ... 0.1 0.1\n", + " ssat (patch, latitude, longitude) float32 108B 0.474 0.439 ... 0.474\n", + " sfc (patch, latitude, longitude) float32 108B 0.294 0.282 ... 0.294\n", + " swilt (patch, latitude, longitude) float32 108B 0.136 0.135 ... 0.136\n", + " hyds (patch, latitude, longitude) float32 108B 8.33e-05 ... 8.33e-05\n", + " sucs (patch, latitude, longitude) float32 108B -0.3795 ... -0.3795\n", + " css (patch, latitude, longitude) float32 108B 303.0 303.0 ... 303.0\n", + " rhosoil (patch, latitude, longitude) float32 108B 1.3e+03 ... 1.3e+03\n", + "Attributes:\n", + " Production: 2025/04/09 at 17:45:55\n", + " Source: CABLE LSM output file\n", + " CABLE_input_file: \n", + " Output_averaging: monthly\n", + " history: Tue Apr 15 10:45:20 2025: /g/data/rp23/experiments/202..." + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "trendy_cable_dataset = run_cable_workflow(trendy_cable_output, landmask, cable_vars)\n", + "trendy_cable_dataset = trendy_cable_dataset.load()\n", + "trendy_cable_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "silt: [0.]\n", + "clay: [0.]\n", + "css: [0.]\n", + "sfc: [0.]\n", + "rhosoil: [0.]\n", + "bch: [0.]\n", + "hyds: [0.]\n", + "ssat: [0.]\n", + "swilt: [0.]\n", + "sucs: [0.]\n" + ] + }, + { + "ename": "KeyError", + "evalue": "\"No variable named 'mvg'. Variables on the dataset include ['longitude', 'latitude', 'iveg', 'bch', 'clay', ..., 'hyds', 'sucs', 'css', 'rhosoil', 'time']\"", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/core/dataset.py\u001b[0m in \u001b[0;36m?\u001b[0;34m(self, name)\u001b[0m\n\u001b[1;32m 1475\u001b[0m \u001b[0mvariable\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_variables\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1476\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1477\u001b[0;31m \u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvariable\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_get_virtual_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_variables\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msizes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1478\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mKeyError\u001b[0m: 'mvg'", + "\nDuring handling of the above exception, another exception occurred:\n", + "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/core/dataset.py\u001b[0m in \u001b[0;36m?\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 1574\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_construct_dataarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1575\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1576\u001b[0;31m raise KeyError(\n\u001b[0m\u001b[1;32m 1577\u001b[0m \u001b[0;34mf\"No variable named {key!r}. Variables on the dataset include {shorten_list_repr(list(self.variables.keys()), max_items=10)}\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/core/dataset.py\u001b[0m in \u001b[0;36m?\u001b[0;34m(self, name)\u001b[0m\n\u001b[1;32m 1475\u001b[0m \u001b[0mvariable\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_variables\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1476\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1477\u001b[0;31m \u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvariable\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_get_virtual_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_variables\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msizes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1478\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/core/dataset.py\u001b[0m in \u001b[0;36m?\u001b[0;34m(variables, key, dim_sizes)\u001b[0m\n\u001b[1;32m 208\u001b[0m \u001b[0msplit_key\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkey\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msplit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\".\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 209\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msplit_key\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 210\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 211\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mKeyError\u001b[0m: 'mvg'", + "\nThe above exception was the direct cause of the following exception:\n", + "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m/scratch/tm70/ag9761/tmp/ipykernel_839876/1252957591.py\u001b[0m in \u001b[0;36m?\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mcable_vars\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0muniq_val\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrendy_cable_dataset\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mbios_cable_dataset\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"{v}: {uniq_val}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/g/data/hh5/public/apps/miniconda3/envs/analysis3/lib/python3.10/site-packages/xarray/core/dataset.py\u001b[0m in \u001b[0;36m?\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 1572\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mutils\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhashable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1573\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1574\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_construct_dataarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1575\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1576\u001b[0;31m raise KeyError(\n\u001b[0m\u001b[1;32m 1577\u001b[0m \u001b[0;34mf\"No variable named {key!r}. Variables on the dataset include {shorten_list_repr(list(self.variables.keys()), max_items=10)}\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1578\u001b[0m ) from e\n\u001b[1;32m 1579\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mKeyError\u001b[0m: \"No variable named 'mvg'. Variables on the dataset include ['longitude', 'latitude', 'iveg', 'bch', 'clay', ..., 'hyds', 'sucs', 'css', 'rhosoil', 'time']\"" + ] + } + ], + "source": [ + "for v in cable_vars:\n", + " uniq_val = np.unique(trendy_cable_dataset[v] - bios_cable_dataset[v])\n", + " print(f\"{v}: {uniq_val}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "analysis3", + "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.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}