{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with REMO domains"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The domain module should give some tools to work with preconfigured or user defined domains. Domains are defined as xarray datasets that will contain dimensions and coodinates according to CF-conventions. The pyremo domain module actually uses the [py-cordex](https://py-cordex.readthedocs.io/en/latest/domains.html) domains module in the background just with another set of tables.\n",
    "\n",
    "**NOTE**: Please be aware that a remo domain is usually a little larger than the “official” cordex domain according to the [archive specification](https://is-enes-data.github.io/cordex_archive_specifications.pdf) since the regional model usually accounts for a “buffer” zone where the lateral boundary conditions are nudged."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with domain information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyremo as pr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The domain module contains some useful functions to work with cordex meta data, e.g., you can get some domain grid information using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pr.domain_info(\"EUR-11\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All available cordex domains are in this table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pr.domains.table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `EUR-11` example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The heart of the module are some functions that create a dataset from the grid information, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eur11 = pr.remo_domain(\"EUR-11\", dummy=\"topo\")\n",
    "eur11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `dummy='topo'` argument means, we want a dummy variable in the dataset to see how the domain looks like. For the dummy topography, we use the `cdo topo` operator in the background. So maybe you have to install `python-cdo`, e.g., `conda install -c conda-forge python-cdo`. Working with xarray datasets means, that we can use all the nice functions of xarray including plotting, e.g.,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eur11.topo.plot(cmap=\"terrain\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eur11.topo.plot(x=\"lon\", y=\"lat\", cmap=\"terrain\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define a slightly more sophisticated plotting function that uses cartopy for the right [projection](https://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html) with a rotated pole:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot(da, pole, vmin=None, vmax=None, borders=True, title=\"\"):\n",
    "    \"\"\"plot a domain using the right projection with cartopy\"\"\"\n",
    "    %matplotlib inline\n",
    "    import cartopy.crs as ccrs\n",
    "    import cartopy.feature as cf\n",
    "    import matplotlib.pyplot as plt\n",
    "\n",
    "    plt.figure(figsize=(20, 10))\n",
    "    projection = ccrs.PlateCarree()\n",
    "    transform = ccrs.RotatedPole(pole_latitude=pole[1], pole_longitude=pole[0])\n",
    "    # ax = plt.axes(projection=projection)\n",
    "    ax = plt.axes(projection=transform)\n",
    "    # ax.set_extent([ds_sub.rlon.min(), ds_sub.rlon.max(), ds_sub.rlat.min(), ds_sub.rlat.max()], crs=transform)\n",
    "    ax.gridlines(\n",
    "        draw_labels=True,\n",
    "        linewidth=0.5,\n",
    "        color=\"gray\",\n",
    "        xlocs=range(-180, 180, 10),\n",
    "        ylocs=range(-90, 90, 5),\n",
    "    )\n",
    "    da.plot(\n",
    "        ax=ax,\n",
    "        cmap=\"terrain\",\n",
    "        transform=transform,\n",
    "        vmin=vmin,\n",
    "        vmax=vmax,\n",
    "        x=\"rlon\",\n",
    "        y=\"rlat\",\n",
    "    )\n",
    "    ax.coastlines(resolution=\"50m\", color=\"black\", linewidth=1)\n",
    "    if borders:\n",
    "        ax.add_feature(cf.BORDERS)\n",
    "    ax.set_title(\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pole_eur = (\n",
    "    eur11.rotated_latitude_longitude.grid_north_pole_longitude,\n",
    "    eur11.rotated_latitude_longitude.grid_north_pole_latitude,\n",
    ")\n",
    "pole_eur"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(eur11.topo, pole_eur)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User defined domain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The domains are created using the `create_dataset` function from the [`py-cordex`](https://py-cordex.readthedocs.io) package, e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cordex import create_dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create the EUR-11 domain manually from the numbers in the table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eur11_user = create_dataset(\n",
    "    nlon=433,\n",
    "    nlat=433,\n",
    "    dlon=0.11,\n",
    "    dlat=0.11,\n",
    "    ll_lon=-28.925,\n",
    "    ll_lat=-23.925,\n",
    "    pollon=-162.00,\n",
    "    pollat=39.25,\n",
    "    dummy=\"topo\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check that this gives the same result as our preconfigured domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eur11_user.equals(eur11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can now use the `create_dataset` function to create any domain as an xarray dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "afr11 = pr.remo_domain(\"AFR-11\", dummy=\"topo\")\n",
    "afr11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pole_afr = (\n",
    "    afr11.rotated_latitude_longitude.grid_north_pole_longitude,\n",
    "    afr11.rotated_latitude_longitude.grid_north_pole_latitude,\n",
    ")\n",
    "pole_afr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(afr11.topo, pole_afr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cropping the REMO domain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes it might be neccessary to crop the REMO data to the official CORDEX grid size, e.g., for cmorization. This can now easily be done like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cordex import cordex_domain\n",
    "\n",
    "eur11_cordex = cordex_domain(\"EUR-11\", dummy=\"topo\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# crop = eur11.sel(rlon=slice(eur11_cordex.rlon.min(), eur11_cordex.rlon.max()), rlat=slice(eur11_cordex.rlat.min(), eur11_cordex.rlat.max()))\n",
    "crop = eur11.sel(rlon=eur11_cordex.rlon, rlat=eur11_cordex.rlat, method=\"nearest\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(crop.topo, pole_eur)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crop"
   ]
  }
 ],
 "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.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}