notebooks/accessing_data_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "583bc943-739f-401d-8727-00495979ed56",
   "metadata": {},
   "source": [
    "<img width='50' src='https://carbonplan-assets.s3.amazonaws.com/monogram/dark-small.png' style='margin-left:0px;margin-top:20px'/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ab3a930-5bcf-410f-8906-fb744f593982",
   "metadata": {},
   "source": [
    "# Accessing CarbonPlan CMIP6 downscaled climate datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4c7235a-0a96-4fe1-a983-5e30f34d3a3b",
   "metadata": {},
   "source": [
    "Authors: Oriana Chegwidden and Max Jones"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01f99de5-12d9-49a4-8eb8-bc55503a04ce",
   "metadata": {},
   "source": [
    "This notebook offers users examples of accessing and working with CarbonPlan's downscaled climate datasets. The dataset collection is further described in an [explainer article](https://carbonplan.org/research/cmip6-downscaling-explainer). Monthly and annual summaries of the data products are visible in an [interactive mapping tool](https://cmip6.carbonplan.org/). We recommend using Python to interact with the datasets. Below we show examples of reading the data, performing basic visualization, and downloading subsets in space and time. We welcome further requests for interaction and encourage [feedback via GitHub](/~https://github.com/carbonplan/cmip6-downscaling/issues)!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b837672-dea0-41e5-ae20-191596ef16c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import regionmask\n",
    "import cartopy.crs as ccrs\n",
    "from carbonplan import styles  # noqa: F401\n",
    "import intake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8100767",
   "metadata": {},
   "outputs": [],
   "source": [
    "from cmip6_downscaling.analysis.analysis import (\n",
    "    grab_big_city_data,\n",
    "    load_big_cities,\n",
    ")\n",
    "from cmip6_downscaling.analysis.plot import plot_city_data\n",
    "\n",
    "xr.set_options(keep_attrs=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "066d5b6d-c209-473d-b689-3d6c8e6bfcbe",
   "metadata": {},
   "source": [
    "## Loading the data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f0e99c9-83a5-4027-bda8-b6a90f56632f",
   "metadata": {},
   "source": [
    "Let's load in the catalog of datasets available in this release."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4655fb8-eb57-4218-b060-8f8b538db39b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat = intake.open_esm_datastore(\n",
    "    \"https://rice1.osn.mghpcc.org/carbonplan/cp-cmip/version1/catalog/osn-rechunked-global-downscaled-cmip6.json\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60b0d86b-d6c4-4890-aa62-32999d1b217b",
   "metadata": {},
   "source": [
    "We can inspect the contents according to different methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0bf7262-1276-4393-aba9-aacf91e3fb4b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cat_subset = cat.search(method=\"GARD-SV\")\n",
    "cat_subset.df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c99a53d6-3dfe-4c4d-98fe-a6cd3bde8289",
   "metadata": {},
   "source": [
    "Now let's specify what models we're interested in. We're going to select a daily maximum temperature run from the `MRI-ESM2-0` GCM and the `SSP2-4.5` future scenario, downscaled using the `GARD-MV` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24848c46-7afe-40f8-8aed-06c641eb8dd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat_subset = cat.search(\n",
    "    method=\"GARD-MV\",\n",
    "    source_id=\"MRI-ESM2-0\",\n",
    "    experiment_id=\"ssp245\",\n",
    "    variable_id=\"tasmax\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29189ef6-9ad9-45f7-92c5-db6c5c9acb6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dsets = cat_subset.to_dataset_dict()\n",
    "dsets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "992c93ae-547d-48c4-8255-c9cb61ded79c",
   "metadata": {},
   "source": [
    "And now let's load that dataset into the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4a180ba-3103-4f4a-bdfc-e1a176e11c57",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = dsets[\"ScenarioMIP.MRI.MRI-ESM2-0.ssp245.day.GARD-MV\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdd9adbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "china_region = {'lat': slice(18, 54), 'lon': slice(17, 135)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6953e42",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.tasmax.isel(time=0).sel(china_region).plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "579f82be",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.sel(china_region).tasmax.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7d6fc0c-0e3f-40aa-8cbf-08db4cdac81d",
   "metadata": {},
   "source": [
    "## Visualizing the data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c40f76a-6751-45b0-8ea9-fae883101c7f",
   "metadata": {},
   "source": [
    "The temperature data are in units of Kelvin. Let's convert to Celcius to make it easer to interpret."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1595177-5f87-4578-8a14-a322533e93b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds -= 273.15"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90954da7-3c97-4ebe-a55e-968d4f521004",
   "metadata": {},
   "source": [
    "### Plotting maps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b175e45-1734-46de-a4fe-9ba85b6c2738",
   "metadata": {},
   "source": [
    "Let's load in a land mask and a projection for mapping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b543059a-6f7d-4645-a591-b465bc9c5f0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "land = regionmask.defined_regions.natural_earth_v5_0_0.land_110\n",
    "projection = ccrs.PlateCarree()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "256d548b-0cae-4eff-a93e-55384cbf8c14",
   "metadata": {},
   "source": [
    "We'll mask out the ocean values and load data for a single timestep (August 1, 2089) and a region of interest. We'll start with the East Africa region we reference in Figure 1 of the companion web article."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "627aad61-7c91-4d87-9642-be0a671b4b22",
   "metadata": {},
   "outputs": [],
   "source": [
    "east_africa_region = {\"lat\": slice(-3, 17), \"lon\": slice(17, 57)}\n",
    "east_africa_tasmax = ds.tasmax.sel(time=\"2089-08-01\").sel(**east_africa_region).load()\n",
    "east_africa_tasmax = east_africa_tasmax.where(land.mask(east_africa_tasmax) == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d905905b-f095-4ac1-b859-b0945a77ecb0",
   "metadata": {},
   "source": [
    "And now let's plot that single timeslice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26e0292d-628d-4a82-9e01-698013a646e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(subplot_kw=dict(projection=projection), figsize=(10, 4))\n",
    "east_africa_tasmax.plot(\n",
    "    cbar_kwargs=dict(label=r\"Maximum temperature $^\\circ$C\"),\n",
    "    cmap=\"warm_dark\",\n",
    ")\n",
    "ax.coastlines()\n",
    "ax.set_xticks([20, 30, 40, 50], crs=projection)\n",
    "ax.set_xlabel(r\"Longitude ($^\\circ$E)\")\n",
    "ax.set_yticks([0, 5, 10, 15], crs=projection)\n",
    "ax.set_ylabel(r\"Latitude ($^\\circ$N)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05ee0856-7b04-4eb9-bbe6-47463c4e0a60",
   "metadata": {},
   "source": [
    "Now let's do the same thing but for the whole globe. And while we're at it, let's grab the annual mean temperatures of that daily data for a 30 year period at the end of the century (the '2080s'). Caution - this could take a while (i.e., minutes)- it's a lot of data! We'll define two timeslices (we'll use the second one later in the notebook)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82ed1928-42e0-449f-89f9-862f780ba592",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_slices = {\"2030s\": slice(\"2020\", \"2049\"), \"2080s\": slice(\"2070\", \"2099\")}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ffd94ad-73c4-4938-a6b5-4d12a19d889a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the 30-year mean bracketing the 2080s.\n",
    "tasmax_2080s = ds.tasmax.sel(time=time_slices[\"2080s\"]).mean(dim=\"time\").load()\n",
    "# mask as above\n",
    "tasmax_2080s = tasmax_2080s.where(land.mask(tasmax_2080s) == 0).load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92f8106b-6091-4485-a9ae-a6685345ddf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(subplot_kw=dict(projection=projection), figsize=(10, 4))\n",
    "tasmax_2080s.plot(\n",
    "    cbar_kwargs=dict(label=r\"Mean daily maximum temperature $^\\circ$C\"),\n",
    "    cmap=\"warm_dark\",\n",
    ")\n",
    "ax.coastlines()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b197b4b4-7a81-4223-8e17-d690dbc5c738",
   "metadata": {},
   "source": [
    "### Plotting timeseries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e790b6f-39a7-4f77-a409-0d5cdcd17580",
   "metadata": {},
   "source": [
    "Let's look at the data at some individual points! After all, one of the main goals of downscaling is to provide more local information. We'll grab timeseries from 20 big cities around the world to explore what climate change might look like for them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "884a6506-0b9d-4185-9b66-74a59ee8560b",
   "metadata": {},
   "outputs": [],
   "source": [
    "big_cities = load_big_cities(num_cities=20, add_additional_cities=False, plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b4a88c8-f9b4-40d6-a4f0-a44098d40185",
   "metadata": {},
   "outputs": [],
   "source": [
    "[downscaled_cities] = grab_big_city_data([ds], big_cities)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb70d3e5-58f8-427b-bd28-f5a2aa34fcf9",
   "metadata": {},
   "source": [
    "Let's plot a timeseries of the daily data at just one of those 20 cities. Let's look at Tokyo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfc73f0b-966c-4698-af59-f644126f2467",
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = downscaled_cities.sel(cities=\"Tokyo\").tasmax\n",
    "ts.plot()\n",
    "plt.ylabel(r\"Daily maximum temperature $^\\circ$C\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd8b25e5-953f-495e-8666-7a3acef80360",
   "metadata": {},
   "source": [
    "That's a lot of daily data though- we can make it clearer by summarizing it into a seasonal cycle. Let's comparing 30 year periods near the start and end of the 21st century to see this model's projection of changes in temperature at this location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "004863c6-29de-4114-9630-e83ac4acc3c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for label, time_slice in time_slices.items():\n",
    "    ts.sel(time=time_slice).groupby(\"time.month\").mean().plot(label=label)\n",
    "ax.set_xticks(np.arange(1, 13))\n",
    "ax.set_xticklabels([\"J\", \"F\", \"M\", \"A\", \"M\", \"J\", \"J\", \"A\", \"S\", \"O\", \"N\", \"D\"])\n",
    "plt.legend()\n",
    "plt.ylabel(r\"Mean daily maximum temperature $^\\circ$C\")\n",
    "plt.xlabel(\"\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aad6900-c51d-467a-9693-214fb1cacf25",
   "metadata": {},
   "source": [
    "Now let's do that for all of the cities we grabbed to get a sense of how popualtion centers around the world will fare. We'll first plot the seasonal cycle to show the projected change from the 2030s to the 2080s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d3f34e4-37bb-453d-8463-372a240711cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_city_data(\n",
    "    downscaled_cities.tasmax,\n",
    "    time_slices=time_slices,\n",
    "    aggregation=\"seasonal_cycle\",\n",
    "    ylabel=r\"Mean daily maximum temperature $^\\circ$C\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46228b38-c55f-4293-9dcb-6f45afa7a24a",
   "metadata": {},
   "source": [
    "Now let's look at the annual means for the entire 21st century."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce2b027d-73e8-4dbd-846b-ce03b6777b27",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_city_data(\n",
    "    downscaled_cities.tasmax,\n",
    "    aggregation=\"annual\",\n",
    "    ylabel=r\"Mean daily maximum temperature $^\\circ$C\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9111a5f0-8151-4539-9a52-492388a93217",
   "metadata": {},
   "source": [
    "## Downloading the data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30ed71f5-1b31-4332-a930-cd4a844f4e55",
   "metadata": {},
   "source": [
    "And now let's download one of these daily timeseries to work with on our own computer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2168ea63-533a-43dd-8e77-3379b2082fb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ts.to_dataframe().reset_index().drop(columns=[\"member_id\", \"cities\", \"lat\", \"lon\"]).set_index(\n",
    "    \"time\"\n",
    ").to_csv(\"tokyo.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b12834e-12cf-4eb0-b2a0-1c6402eb75e3",
   "metadata": {},
   "source": [
    "We can also download a small regional subset as a netcdf file locally. Note that this might be quite large. We'll first check out how big it is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06db6d73-b7b7-40e6-b426-55ddef04ffdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Dataset is {} GB\".format(ds.sel(**east_africa_region).nbytes * 1e-9))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17a34d91-712f-47a3-934f-b57d01a09d61",
   "metadata": {},
   "source": [
    "If you want to save it locally switch the flag to `True`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10131b01-578d-4ca9-b1e0-1c79786a2414",
   "metadata": {},
   "outputs": [],
   "source": [
    "save_subset = False\n",
    "if save_subset:\n",
    "    ds.sel(**east_africa_region).to_netcdf(\"region.nc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e519f0f-565b-4bba-b292-37481552ba3d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}