{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Grass growth anomalies\n", "\n", "- Weighted means take into account the number of days in each month" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "metadata": {} }, "outputs": [], "source": [ "import glob\n", "import importlib\n", "import itertools\n", "import os\n", "import sys\n", "from datetime import datetime, timezone\n", "import geopandas as gpd\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import xarray as xr\n", "import climag.climag as cplt\n", "from climag import climag_plot\n", "import seaborn as sns\n", "import fiona\n", "import pandas as pd\n", "import datetime" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "metadata": {} }, "outputs": [], "source": [ "exp_list = [\"historical\", \"rcp45\", \"rcp85\"]\n", "model_list = [\"CNRM-CM5\", \"EC-EARTH\", \"HadGEM2-ES\", \"MPI-ESM-LR\"]\n", "dataset_list = [\"HiResIreland\"]" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "metadata": {} }, "outputs": [], "source": [ "def keep_minimal_vars(data):\n", " \"\"\"\n", " Drop variables that are not needed\n", " \"\"\"\n", "\n", " data = data.drop_vars(\n", " [\n", " \"bm_gv\",\n", " \"bm_gr\",\n", " \"bm_dv\",\n", " \"bm_dr\",\n", " \"age_gv\",\n", " \"age_gr\",\n", " \"age_dv\",\n", " \"age_dr\",\n", " \"omd_gv\",\n", " \"omd_gr\",\n", " \"lai\",\n", " \"env\",\n", " \"wr\",\n", " \"aet\",\n", " \"sen_gv\",\n", " \"sen_gr\",\n", " \"abs_dv\",\n", " \"abs_dr\",\n", " \"gro\",\n", " \"c_bm\",\n", " # \"bm\",\n", " \"pgro\",\n", " \"i_bm\",\n", " \"h_bm\",\n", " ]\n", " )\n", "\n", " return data" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Corine land cover 2018\n", "# pasture only - vectorised (done in QGIS)\n", "pasture = gpd.read_file(\n", " os.path.join(\"data\", \"landcover\", \"clc-2018-pasture.gpkg\"),\n", " layer=\"dissolved\",\n", ").to_crs(cplt.ITM_EPSG)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "metadata": {} }, "outputs": [], "source": [ "def combine_datasets(dataset_dict, dataset_crs):\n", " dataset = xr.combine_by_coords(\n", " dataset_dict.values(), combine_attrs=\"override\"\n", " )\n", " dataset.rio.write_crs(dataset_crs, inplace=True)\n", "\n", " return dataset" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "metadata": {} }, "outputs": [], "source": [ "def generate_stats(dataset):\n", " ds = {}\n", " # ds_mean ={}\n", " # ds_max = {}\n", "\n", " for exp, model in itertools.product(exp_list, model_list):\n", " # auto-rechunking may cause NotImplementedError with object dtype\n", " # where it will not be able to estimate the size in bytes of object\n", " # data\n", " if model == \"HadGEM2-ES\":\n", " CHUNKS = 300\n", " else:\n", " CHUNKS = \"auto\"\n", "\n", " ds[f\"{model}_{exp}\"] = xr.open_mfdataset(\n", " glob.glob(\n", " os.path.join(\n", " \"data\",\n", " \"ModVege\",\n", " dataset,\n", " exp,\n", " model,\n", " f\"*{dataset}*{model}*{exp}*.nc\",\n", " )\n", " ),\n", " chunks=CHUNKS,\n", " decode_coords=\"all\",\n", " )\n", "\n", " # copy CRS\n", " crs_ds = ds[f\"{model}_{exp}\"].rio.crs\n", "\n", " # remove spin-up year\n", " if exp == \"historical\":\n", " ds[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].sel(\n", " time=slice(\"1976\", \"2005\")\n", " )\n", " else:\n", " ds[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].sel(\n", " time=slice(\"2041\", \"2070\")\n", " )\n", "\n", " # convert HadGEM2-ES data back to 360-day calendar\n", " # this ensures that the correct weighting is applied when\n", " # calculating the weighted average\n", " if model == \"HadGEM2-ES\":\n", " ds[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].convert_calendar(\n", " \"360_day\", align_on=\"year\"\n", " )\n", "\n", " # December data\n", " ds[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].sel(\n", " time=ds[f\"{model}_{exp}\"][\"time\"].dt.month.isin([12])\n", " )\n", "\n", " # assign new coordinates and dimensions\n", " ds[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].assign_coords(exp=exp)\n", " ds[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].expand_dims(dim=\"exp\")\n", " ds[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].assign_coords(model=model)\n", " ds[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].expand_dims(dim=\"model\")\n", "\n", " # drop unnecessary variables\n", " ds[f\"{model}_{exp}\"] = keep_minimal_vars(data=ds[f\"{model}_{exp}\"])\n", "\n", " ds[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].rio.clip(\n", " pasture[\"geometry\"].to_crs(ds[f\"{model}_{exp}\"].rio.crs),\n", " all_touched=True,\n", " )\n", "\n", " # weighted mean\n", " weights = (\n", " ds[f\"{model}_{exp}\"][\"time\"].dt.days_in_month.groupby(\"time.year\")\n", " / ds[f\"{model}_{exp}\"][\"time\"]\n", " .dt.days_in_month.groupby(\"time.year\")\n", " .sum()\n", " )\n", "\n", " # test that the sum of weights for each season is one\n", " np.testing.assert_allclose(\n", " weights.groupby(\"time.year\").sum().values,\n", " np.ones(len(set(weights[\"year\"].values))),\n", " )\n", "\n", " # calculate the weighted average\n", " ds[f\"{model}_{exp}\"] = (\n", " (ds[f\"{model}_{exp}\"] * weights)\n", " .groupby(\"time.year\")\n", " .sum(dim=\"time\")\n", " )\n", "\n", " # # max\n", " # ds_max[f\"{model}_{exp}\"] = ds[f\"{model}_{exp}\"].groupby(\"time.year\").max(dim=\"time\")\n", "\n", " # combine data\n", " ds = combine_datasets(ds, crs_ds)\n", " # ds_max = combine_datasets(ds_max, crs_ds)\n", "\n", " # ensemble mean\n", " # ds = ds.mean(dim=\"model\", skipna=True)\n", "\n", " # long-term average\n", " # ds_lta = ds.mean(dim=\"year\", skipna=True)\n", "\n", " return ds # ds_mean, ds_max #, ds_lta" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "metadata": {} }, "outputs": [], "source": [ "hiresireland = generate_stats(\"HiResIreland\")" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.Dataset> Size: 293MB\n",
"Dimensions: (year: 60, model: 4, exp: 3, rlat: 113, rlon: 90)\n",
"Coordinates:\n",
" * rlon (rlon) float64 720B -1.68 -1.645 -1.61 ... 1.365 1.4 1.435\n",
" * rlat (rlat) float64 904B -1.945 -1.91 -1.875 ... 1.905 1.94 1.975\n",
" * model (model) <U10 160B 'CNRM-CM5' 'EC-EARTH' ... 'MPI-ESM-LR'\n",
" * year (year) int64 480B 1976 1977 1978 1979 ... 2067 2068 2069 2070\n",
" * exp (exp) <U10 120B 'historical' 'rcp45' 'rcp85'\n",
" lon (rlat, rlon) float32 41kB dask.array<chunksize=(113, 90), meta=np.ndarray>\n",
" lat (rlat, rlon) float32 41kB dask.array<chunksize=(113, 90), meta=np.ndarray>\n",
" spatial_ref int64 8B 0\n",
"Data variables:\n",
" bm (year, model, exp, rlat, rlon) float64 59MB dask.array<chunksize=(1, 1, 1, 113, 90), meta=np.ndarray>\n",
" gro (year, model, exp, rlat, rlon) float64 59MB dask.array<chunksize=(1, 1, 1, 113, 90), meta=np.ndarray>\n",
" i_bm (year, model, exp, rlat, rlon) float64 59MB dask.array<chunksize=(1, 1, 1, 113, 90), meta=np.ndarray>\n",
" h_bm (year, model, exp, rlat, rlon) float64 59MB dask.array<chunksize=(1, 1, 1, 113, 90), meta=np.ndarray>\n",
" c_bm (year, model, exp, rlat, rlon) float64 59MB dask.array<chunksize=(1, 1, 1, 113, 90), meta=np.ndarray>| \n", " | model | \n", "year | \n", "exp | \n", "bm | \n", "gro | \n", "i_bm | \n", "h_bm | \n", "c_bm | \n", "
|---|---|---|---|---|---|---|---|---|
| 0 | \n", "CNRM-CM5 | \n", "1976 | \n", "historical | \n", "592.320147 | \n", "0.404723 | \n", "3287.816394 | \n", "170.450066 | \n", "0.082846 | \n", "
| 1 | \n", "CNRM-CM5 | \n", "1977 | \n", "historical | \n", "696.484553 | \n", "0.079848 | \n", "2967.021082 | \n", "30.719180 | \n", "0.677006 | \n", "
| 2 | \n", "CNRM-CM5 | \n", "1978 | \n", "historical | \n", "685.613975 | \n", "0.151322 | \n", "3034.121935 | \n", "8.065548 | \n", "0.259412 | \n", "
| 3 | \n", "CNRM-CM5 | \n", "1979 | \n", "historical | \n", "823.039920 | \n", "0.150739 | \n", "3151.410164 | \n", "52.672375 | \n", "0.343000 | \n", "
| 4 | \n", "CNRM-CM5 | \n", "1980 | \n", "historical | \n", "521.447003 | \n", "0.522444 | \n", "2299.919977 | \n", "66.140569 | \n", "0.038727 | \n", "