2019_climate4R_SU.8.1_A.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [SU.8.1 ADAPTECCA SCENARIOS: use of regional climate change scenarios for impact and adaptation studies](https://web.unican.es/cursosdeveranoyextension/cursos-de-verano/curso?c=2861)\n",
    "\n",
    "## Tools for accessing and processing climate data: Case study with R\n",
    "\n",
    "### Herramientas para el acceso y procesamiento de datos climáticos: Caso práctico con R\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****\n",
    "![c4R](shared/figs/climate4R_esquema.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This worked example contains the code that reproduces part of the examples shown in the paper [\"climate4R: An R-based Framework for Climate Data Access, Post-processing and Bias Correction\"](https://www.sciencedirect.com/science/article/pii/S1364815218303049)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(loadeR)\n",
    "library(transformeR)\n",
    "library(visualizeR)\n",
    "library(downscaleR)\n",
    "library(climate4R.climdex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# 1. Cliamte data loading from OPeNDAP server: E-OBS observational data \n",
    "## Loading, collocating and harmonizing data \n",
    "\n",
    "\n",
    "The domain of the study area is defined by the following bounding coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon <- c(-10, 5)\n",
    "lat <- c(36, 44)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The SU index (summer days) can be obtained on-the-fly by loading maximum temperature data with function `loadGridData` and by the following argument settings: `aggr.m = \"sum\"`, `condition = \"GT\"` and `threshold = 25`. First we load E-OBS observational data by pointing to a NetCDF file via OPeNDAP. Previous to loading, function `dataInventory` might be applied for an overview of the dataset, which returns an inventory (object `di`) of the available variables names, units, coordinates, etc.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eobs<-\"http://opendap.knmi.nl/knmi/thredds/dodsC/e-obs_0.25regular/tx_0.25deg_reg_v17.0.nc\"\n",
    "di <- dataInventory(eobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "str(di)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the NetCDF file contains maximum temperature data named as \"tx\", thus, we set `var = \"tx\"` when calling to `loadGridData`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SU <- loadGridData(eobs, var = \"tx\",\n",
    "                    season = 1:12, \n",
    "                    years = 1971:2000,\n",
    "                    lonLim = lon, \n",
    "                    latLim = lat,\n",
    "                    aggr.m = \"sum\", \n",
    "                    condition = \"GT\", \n",
    "                    threshold = 25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "?loadGridData"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using a dictionary\n",
    "\n",
    "In order to load and work with harmonized data we can repeat the above operation using a dictionary file, that defines the necessary name and unit transformations to the standard parameters. Function `C4R.vocabulary` displays the climate4R standard variable naming and units:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C4R.vocabulary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the only non-standard parameter in the E-OBS dataset is the variable name (\"tx\"), however, we could perform further loading requests using the standard name if a dictionary file is crated previously (see the [`loadeR` wiki](https://github.com/SantanderMetGroup/loadeR/wiki/Harmonization)). This can be done easily, for instance, in the following manner:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file.create(\"dicEOBS.dic\")\n",
    "writeLines(c(\"identifier,short_name,time_step,lower_time_bound,upper_time_bound, cell_method,offset,scale,deaccum,derived,interface\",\n",
    "             \"tasmax,tx,24h,0,24,max,0,1,0,0,\"), \"dicEOBS.dic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next the loading operation is repeated but using the standard name for the maximum temperature (`var = \"tasmax\"`) and by passing the path to our *.dic file (\"dicEOBS.dic\") in argument `dictionary`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SU <- loadGridData(eobs,\n",
    "                         var = \"tasmax\",\n",
    "                         season = 1:12,\n",
    "                         lonLim = lon,\n",
    "                         latLim = lat,\n",
    "                         years = 1971:2000,\n",
    "                         aggr.m = \"sum\", \n",
    "                         threshold = 25,\n",
    "                         condition = \"GT\",\n",
    "                         dictionary = \"dicEOBS.dic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Transformation and visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An useful plotting function is `temporalPlot` that displays temporal series of multiple datasets\n",
    "and periods on the same plot. Here we plot the series corresponding to a single grid box (`latLim = 37.89`, `lonLim = -4.78`). If several grid boxes are considered (e.g. the whole domain) `temporalPlot` performs the spatial (`lat` and `lon` dimensions) aggregation before plotting (the `mean` is computed by default, type `?temporalPlot`).\n",
    "\n",
    "Note that function `temporalPlot` is based on `lattice` and arguments from function `xyplot` are optionally passed to argument `xyplot.custom`, allowing for a fine tuning of multiple graphical parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cordoba\n",
    "temporalPlot(SU, lonLim = -4.78, latLim = 37.89, xyplot.custom = list(ylim = c(0, 40)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Macizo de Aralar (cerca de San Sebastian)\n",
    "temporalPlot(SU, lonLim = -2, latLim = 43, xyplot.custom = list(ylim = c(0, 40)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Note that `loadGridData` returns monthly summer days (SU). To compute the annual index we only need to apply function `aggregateGrid` that performs the aggregation of the desired data dimension (in this case `time`). We use argument `aggr.y` to perform annual aggregation with function `sum`:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SU.annual <- aggregateGrid(SU, aggr.y = list(FUN = \"sum\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Macizo de Aralar (cerca de San Sebastian)\n",
    "temporalPlot(SU.annual, lonLim = -2, latLim = 43, xyplot.custom = list(ylim = c(15, 75)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ETCCDI index calculation (SU) from daily data\n",
    "\n",
    "Daily data of the original variable is loaded as above but excluding arguments `aggr.m`, `threshold` y `condition`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TX.obs <- loadGridData(eobs,\n",
    "                         var = \"tasmax\",\n",
    "                         season = 1:12,\n",
    "                         lonLim = lon,\n",
    "                         latLim = lat,\n",
    "                         years = 1971:2000,\n",
    "                         #aggr.m = \"sum\", \n",
    "                         #threshold = 25,\n",
    "                         #condition = \"GT\",\n",
    "                         dictionary = \"dicEOBS.dic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For plotting maps function `spatialPlot` is used. For plotting the maximum temperature here we select the `\"Spectral\"` color theme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(TX.obs), backdrop.theme = \"countries\",\n",
    "            at = seq(5, 35, 1), color.theme = \"Spectral\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the ETCCDI indices: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "climdexShow()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next the raw SU index (summer days) is calculated (object `SU.obs`), in a single line with function `climdexGrid`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SU.obs <- climdexGrid(tx = TX.obs, index.code = \"SU\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Macizo de Aralar (cerca de San Sebastian)\n",
    "temporalPlot(SU.obs, lonLim = -2, latLim = 43, xyplot.custom = list(ylim = c(15, 75)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case we set the `\"RdYlBu\"` color theme to visualize the mean annual SU for the reference period (1971-2000).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(SU.annual), backdrop.theme = \"countries\", \n",
    "            color.theme = \"RdYlBu\", rev.colors = TRUE,\n",
    "            at = seq(0, 260, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# 2. Cliamte data loading from local files: CORDEX historical projections\n",
    "Next, projection data (for both the historical and the RCP8.5 scenarios) is loaded from local NetCDF files, which correspond to a particular RCM (Regional Climate Model ICHEC-EC-EARTH_r12i1p1_SMHI-RCA4_v1) from EURO-CORDEX. These files were downloaded from ESGF (see Appendix A in the manuscript) and stored locally. Next we list them in objects `dir` and `dirf`, the first corresponding to the historical scenario and the second to the future RCP8.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir <- \"shared/local_nc_data/historical/\"\n",
    "dirf <- \"shared/local_nc_data/rcp85/\"\n",
    "print(list.files(dir, recursive = T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each file in the list contains data for a 5-year period of the same variable (tasmax). Therefore, we use a \"catalog\" (*.ncml file) to load data for the required period without worrying about the different files that need to be read and bound. Next we create two catalogs (for each scenario) with function `makeAggregateDataset` (\"CDX_hist.ncml\" and \"CDX_rcp85.ncml\"):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "makeAggregatedDataset(source.dir = dir, recursive = T, ncml.file = \"data/CDX_hist.ncml\")\n",
    "makeAggregatedDataset(source.dir = dirf, recursive = T, ncml.file = \"data/CDX_rcp85.ncml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The created *.ncml files are then used as a single access point to load data and to do the data inventory as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "di <- dataInventory(\"data/CDX_hist.ncml\")\n",
    "str(di$tasmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Contrarily to the case of the E-OBS dataset, the variable name is standard, but not the units (K). Therefore we define the harmonization parameters in another dictionary file (\"dicCDX.dic\"), where the offset is -273.15 to convert the data to the standard units (ºC):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file.create(\"dicCDX.dic\")\n",
    "writeLines(c(\"identifier,short_name,time_step,lower_time_bound,upper_time_bound,cell_metod,offset,scale,deaccum,derived,interface\",\n",
    "             \"tasmax,tasmax,24h,0,24,max,-273.15,1,0,0,\"), \"dicCDX.dic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Historical data\n",
    "Next, harmonized data is loaded for a single CORDEX model, for the historical scenario and the same reference period used to load E-OBS observational data (1971-2000):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TX.hist <- loadGridData(dataset = \"data/CDX_hist.ncml\",\n",
    "                     var = \"tasmax\",\n",
    "                     season = 1:12,\n",
    "                     lonLim = lon,\n",
    "                     latLim = lat,\n",
    "                     years = 1971:2000,\n",
    "                     dictionary = \"dicCDX.dic\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(TX.hist), backdrop.theme = \"countries\",\n",
    "            color.theme = \"Spectral\", rev.colors = TRUE,\n",
    "            at = seq(5, 35, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be noted in the resulting figure, the spatial grid of CORDEX is different from E-OBS. We can use function `interpGrid` to interpolate CORDEX data to the E-OBS spatial grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TX.hist.i <- interpGrid(TX.hist, getGrid(TX.obs))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(TX.hist.i), backdrop.theme = \"countries\",\n",
    "            color.theme = \"Spectral\", rev.colors = TRUE,\n",
    "            at = seq(5, 35, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nest we calculate the  `SU` index for CORDEX historical data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SU.hist <- climdexGrid(tx = TX.hist.i, index.code = \"SU\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(SU.hist), backdrop.theme = \"countries\",             \n",
    "            color.theme = \"RdYlBu\", rev.colors = TRUE,\n",
    "           at = seq(0, 260, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite not being necessary, here we apply a land mask before calculating the bias in order to eliminate the values projected by the CORDEX model over the sea. To do so, `gridArithmetics` can be used, first to create the mask and second to apply it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eobs.mask <- gridArithmetics(SU.obs, 0, operator = \"*\")\n",
    "spatialPlot(climatology(eobs.mask))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SU.hist <- gridArithmetics(SU.hist, eobs.mask, operator = \"+\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "spatialPlot(climatology(SU.hist), backdrop.theme = \"countries\",             \n",
    "            color.theme = \"RdYlBu\", rev.colors = TRUE,\n",
    "           at = seq(0, 260, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bias is computed by subtracting the SU index of E-OBS (object `SU.obs`) to the SU index of historical CORDEX (object `SU.hist`), for which function `gridArithmetics` is again used. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bias <- gridArithmetics(SU.hist, SU.obs, operator = \"-\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we plot the bias (object `bias`), in this case we select the `\"PiYG\"` color theme:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(bias), backdrop.theme = \"countries\", \n",
    "            at = seq(-100, 100, 10), color.theme = \"PiYG\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Macizo de Aralar (cerca de San Sebastian)\n",
    "temporalPlot(SU.obs, SU.hist, xyplot.custom = list(ylim = c(50, 120)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets also visualize the quantile quantile plot for a single location (Cordoba in this example). Function `subsetGrid` performs spatial subsetting by setting arguments `lonLim` and `latLim`. We can use function `qqplot` for plotting the quantiles of the historical CORDEX data vs observed data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs1 <- subsetGrid(TX.obs, lonLim = -4.78, latLim = 37.89)\n",
    "hist1 <- subsetGrid(TX.hist.i, lonLim = -4.78, latLim = 37.89)\n",
    "qqplot(obs1$Data, hist1$Data)\n",
    "lines(c(0, 50), c(0, 50), col = \"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# 3. Cliamte data loading from UDG: CORDEX future climate change projections\n",
    "\n",
    "Future projection can be loaded following the same steps followed above to load historical data. However, in order to illustrate the third option for data loading, here we load data for the RCP8.5 scenario and future period 2071-2100 via the Santander User Data Gateway (UDG)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display available datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "UDG.datasets(pattern = \"EUR44.*rcp85.*RCA4\")$name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "login:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loginUDG(username = \"miturbide\", password = \"lukinvela9&9\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load harmonized historical CORDEX data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TX.fut <- loadGridData(dataset = \"CORDEX-EUR44_EC-EARTH_r12i1p1_rcp85_RCA4_v1\",\n",
    "                     var = \"tasmax\",\n",
    "                     season = 1:12,\n",
    "                     lonLim = lon,\n",
    "                     latLim = lat,\n",
    "                     years = 2071:2100,\n",
    "                     #dictionary = \"local_nc_data/dicCDX.dic\"\n",
    "                      )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We perfome the same operations of interpolation, index calculation and land/sea masking as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TX.fut.i <- interpGrid(TX.fut, getGrid(TX.obs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SU.fut <- climdexGrid(tx = TX.fut.i, index.code = \"SU\")\n",
    "SU.fut <- gridArithmetics(SU.fut, eobs.mask, operator = \"+\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in this case the application of `gridArithmetics` gives the projected climate change signal (object `CCsignal`) w.r.t the historical period (object `SUh.interp`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "CCsignal <- gridArithmetics(SU.fut, \n",
    "                            SU.hist,\n",
    "                            operator = \"-\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we plot the SU index and the change signal:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(SU.fut), backdrop.theme = \"countries\", \n",
    "            color.theme = \"RdYlBu\", rev.colors = TRUE,\n",
    "            at = seq(0, 260, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(CCsignal), backdrop.theme = \"countries\",\n",
    "            at = seq(0, 80, 5), color.theme = \"Reds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# 4. Post-processing: Bias Correction\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Media espacial\n",
    "temporalPlot(SU.obs, SU.hist, SU.fut, \n",
    "             cols = c(\"black\", \"red\", \"red\"),\n",
    "             xyplot.custom = list(ylim = c(50, 170)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Next the \"eqm\" method is applied to bias correct historical CORDEX data (object `TX.hist`) by means of function `biasCorrection`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TX.hist.bc <- biasCorrection(y = TX.obs, x = TX.hist, \n",
    "                         method = \"eqm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "qqplots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist.bc1 <- subsetGrid(TX.hist.bc, lonLim = -4.78, latLim = 37.89)\n",
    "par(mfrow = c(1, 2))\n",
    "qqplot(obs1$Data, hist1$Data)\n",
    "lines(c(0, 50), c(0, 50), col = \"red\")\n",
    "qqplot(obs1$Data, hist.bc1$Data)\n",
    "lines(c(0, 50), c(0, 50), col = \"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Next the \"eqm\" method is applied to bias correct future CORDEX data (object `TX.fut`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "TX.fut.bc <- biasCorrection(y = TX.obs, x = TX.hist, newdata = TX.fut, \n",
    "                         method = \"eqm\", extrapolation = \"constant\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "qqplots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fut1 <- subsetGrid(TX.fut.i, lonLim = -4.78, latLim = 37.89)\n",
    "fut.bc1 <- subsetGrid(TX.fut.bc, lonLim = -4.78, latLim = 37.89)\n",
    "par(mfrow = c(1, 2))\n",
    "qqplot(obs1$Data, fut1$Data)\n",
    "lines(c(0, 50), c(0, 50), col = \"red\")\n",
    "qqplot(obs1$Data, fut.bc1$Data)\n",
    "lines(c(0, 50), c(0, 50), col = \"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Index calculation is repeated but for bias adjusted data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SU.hist.bc <- climdexGrid(tx = TX.hist.bc, index.code = \"SU\")\n",
    "SU.fut.bc <- climdexGrid(tx = TX.fut.bc, index.code = \"SU\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Media espacial\n",
    "temporalPlot(SU.obs, SU.hist, SU.fut, SU.hist.bc, SU.fut.bc, \n",
    "             cols = c(\"black\", \"red\", \"red\", \"blue\", \"blue\"),\n",
    "             xyplot.custom = list(ylim = c(50, 170)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Climate change signal calculation is also repeated for the corrected data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CCsignal.bc <- gridArithmetics(SU.fut.bc, \n",
    "                            SU.obs,\n",
    "                            operator = \"-\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By plotting the resulting objects we obtain..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(SU.fut.bc), backdrop.theme = \"countries\", \n",
    "            at = seq(0, 260, 10), color.theme = \"RdYlBu\", rev.colors = TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatialPlot(climatology(CCsignal.bc), backdrop.theme = \"countries\",\n",
    "            at = seq(0, 80, 5), color.theme = \"Reds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Export figures using functions `pdf` and `dev.off`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdf(\"CCsignal.bc.pdf\")\n",
    "spatialPlot(climatology(CCsignal.bc), backdrop.theme = \"countries\",\n",
    "            at = seq(0, 80, 5), color.theme = \"Reds\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "****\n",
    "![results](shared/figs/final_figures.pdf.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Other available material\n",
    "\n",
    "* [2018_climate4R_example1.pdf](https://github.com/SantanderMetGroup/notebooks/blob/devel/2018_climate4R_example1.pdf) and [2018_climate4R_example2.pdf](https://github.com/SantanderMetGroup/notebooks/blob/devel/2018_climate4R_example2.pdf) contain the full code for **Examples 1 and 2** of the paper `climate4R: An Ecosystem of R packages for Climate Data Access, Post-processing and Bias Correction'.\n",
    "* Find more worked examples on the utilization of climate4R packages in their respective GitHub **wiki**-s at [https://github.com/SantanderMetGroup](https://github.com/SantanderMetGroup):\n",
    "    + [loadeR: https://github.com/SantanderMetGroup/loadeR/wiki](https://github.com/SantanderMetGroup/loadeR/wiki)\n",
    "    + [transformeR: https://github.com/SantanderMetGroup/transformeR/wiki](https://github.com/SantanderMetGroup/transformeR/wiki)\n",
    "    + [downscaleR: https://github.com/SantanderMetGroup/downscaleR/wiki](https://github.com/SantanderMetGroup/downscaleR/wiki)\n",
    "    + [visualizeR: https://github.com/SantanderMetGroup/visualizeR/wiki](https://github.com/SantanderMetGroup/visualizeR/wiki) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Links for software installation\n",
    "* Jupyter notebook: [https://test-jupyter.readthedocs.io/en/latest/install.html](https://test-jupyter.readthedocs.io/en/latest/install.html)\n",
    "* R: [https://cran.r-project.org/doc/manuals/r-release/R-admin.html] (https://cran.r-project.org/doc/manuals/r-release/R-admin.html](https://cran.r-project.org/doc/manuals/r-release/R-admin.html)\n",
    "* Rstudio: [https://www.rstudio.com/products/rstudio/download/](https://www.rstudio.com/products/rstudio/download/)\n",
    "\n",
    "# climate4R installation \n",
    "[https://github.com/SantanderMetGroup/climate4R](https://github.com/SantanderMetGroup/climate4R)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}