data_example_R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "89781879-bb66-4786-a084-a4388bb2627c",
   "metadata": {},
   "source": [
    "# Climate Model and Climate data workshop, hands-on session\n",
    "\n",
    "This jupyter notebook is part of the climate model and climate data workshop. It showcases how to use climate model data in R. However, it is not a R course!\n",
    "\n",
    "Use kernel Teaching (R SLAC 2021).\n",
    "\n",
    "Here we will use the ncdf4 (https://search.r-project.org/CRAN/refmans/ncdf4/html/00Index.html) library."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec4b8b13-9980-468d-b749-26a4e0933e3b",
   "metadata": {},
   "source": [
    "## 1 Loading and inspecting data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adb1297b-bb66-4d6f-ade6-23f8923c33f2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#-----------------------------------------------------------------------------#\n",
    "# Load libraries\n",
    "#-----------------------------------------------------------------------------#\n",
    "library(ncdf4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2e09982-551d-4bdf-bc6b-6b988bbbce60",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get working directory\n",
    "WD <- getwd()\n",
    "\n",
    "# define directory where data is\n",
    "indir <- \"/net/co2/c2sm-data/rlorenz/climate_model_data_workshop/\"\n",
    "# define output directory in case you want to save any data\n",
    "outdir <- paste(WD, \"/data/\", sep=\"\")\n",
    "print(outdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92d135ee-248d-4431-b240-d5aaf06467b7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#-----------------------------------------------------------------------------#\n",
    "# Read data\n",
    "#-----------------------------------------------------------------------------#\n",
    "file <- paste(indir, \"tas_Amon_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_201001-201412.nc\", sep=\"\")\n",
    "ex.nc <- nc_open(file)\n",
    "print(paste(\"The file has\",ex.nc$nvars,\"variables\"))\n",
    "print(ex.nc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "378543e5-20b0-4428-9eb3-228a49d12c94",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This illustrates how to read all the data from a variable (the 5th) if you do not know the variable name\n",
    "v5 <- ex.nc$var[[5]]\n",
    "data1 <- ncvar_get( ex.nc, v5 )\t# by default, reads ALL the data\n",
    "\n",
    "print(paste(\"Var 5 has name\",v5$name,\"and\",v5$ndim, \"dimensions.\",\n",
    "    \"Dimensions are:\"))\n",
    "print(v5$varsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "378be25c-438f-45a0-8c91-0cd97f39e257",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read variable \"tas\" directly if variable name is known\n",
    "tas <- ncvar_get(ex.nc, \"tas\")\n",
    "print(dim(tas))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfb88c70-1a81-48f8-a2c6-9f97bc74ff41",
   "metadata": {},
   "source": [
    "❓ **Questions**\n",
    "\n",
    "1. How many years are covered in the dataset?\n",
    "2. What is the data frequency?\n",
    "3. How many variables does the dataset contain?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47ca7864-2f01-4eec-97ad-a00f4170942e",
   "metadata": {},
   "source": [
    "### 1.1 Exercise\n",
    "Calculate the grid spacing in latitudinal and longitudinal direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bdcde4a-c920-49b7-8eee-6c546965e1a4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3f674634-60fa-44eb-946f-e68fb9c839d1",
   "metadata": {},
   "source": [
    "## 2. Subsetting data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8289bff1-c9ca-494f-8f25-02bd50e6240d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Read time dimension\n",
    "time <- ncvar_get(ex.nc, \"time\")\n",
    "time_units <- ncatt_get(ex.nc, \"time\", \"units\") #time unit, for how defined\n",
    "# extract base date, third bit of time_units \"days since 1950-01-01 00:00\"\n",
    "base_date <- as.character(lapply(strsplit(as.character(time_units$value),\n",
    "                                          split = \" \"), \"[\", 3))\n",
    "time_d <- as.Date(time, format = \"%j\",\n",
    "                  origin = as.Date(base_date)) # define time as date\n",
    "time.y <- format(time_d, \"%Y\") #extract years only from time_d\n",
    "\n",
    "Y.Start <- time.y[1] # first year in timeseries\n",
    "                     # -> can be used as info for output file\n",
    "Y.End <- time.y[length(time.y)] # last year in timeseries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95f9b0e3-ea26-4238-a32c-0029236ee331",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# print start and end years, convert to integers\n",
    "print(Y.Start)\n",
    "print(Y.End)\n",
    "first_year <- strtoi(Y.Start, base = 0L)\n",
    "last_year <- strtoi(Y.End, base = 0L)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43f072d6-79ae-442d-959d-48c8c5046b8b",
   "metadata": {},
   "source": [
    "### Subset data for a specific year and print the dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8be7b520-ddff-4fc6-aff5-507dddbc4aae",
   "metadata": {},
   "source": [
    "#### Hint\n",
    "R's which method can be used to find indeces where a specific condition is met, see:\n",
    "\n",
    "- https://stat.ethz.ch/R-manual/R-devel/library/base/html/which.html\n",
    "\n",
    "check ncvar_get and what kind of parameters can be used:\n",
    "- https://www.rdocumentation.org/packages/ncdf4/versions/1.21/topics/ncvar_get"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08002105-f71c-4f84-a2f2-7fa7147ecc0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the time slice you are interested in (e.g., time = 2012)\n",
    "time_index <- which(time.y == 2012)\n",
    "print(time_index)\n",
    "\n",
    "# Read the variable data for the selected time slice\n",
    "variable_name <- \"tas\"\n",
    "data <- ncvar_get(ex.nc, variable_name, start = c(1, 1, time_index[1]), count = c(-1, -1, length(time_index)))\n",
    "print(dim(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7f71d98-325c-49c9-81b4-1daada9d7d4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is an example how a netcdf can be saved from R\n",
    "print(\"save monthly data to netcdf\")\n",
    "\n",
    "#define dimensions for output file\n",
    "dim1 <- ncdim_def(\"longitude\", \"degrees_east\", as.double(lon))\n",
    "dim2 <- ncdim_def(\"latitude\", \"degrees_north\", as.double(lat))\n",
    "dimT <- ncdim_def(\"time\", \"months\", time_index, unlim = FALSE)\n",
    "missval <- -9999\n",
    "\n",
    "# define the EMPTY netcdf variable \"tas\"\n",
    "var1 <- ncvar_def(\"tas\", \"K\", list(dim1, dim2, dimT), missval,\n",
    "                  longname = \"temperature\")\n",
    "\n",
    "# associate the netcdf variable with a netcdf file  \n",
    "file.out <- paste(outdir, \"T_2012.nc\", sep = \"\")\n",
    "print(file.out)\n",
    "# in case the output directory does not exist yet, create it\n",
    "if (!file.exists(outdir)){\n",
    "   dir.create(file.path(outdir), showWarnings = FALSE)\n",
    "}\n",
    "\n",
    "# create the netcdf file\n",
    "nc.ex = nc_create(file.out, list(var1))\n",
    "\n",
    "# put data into file\n",
    "ncvar_put(nc.ex, var1, data) # one could write a subset of the data using\n",
    "                           # start = c(1, 1, 1),  count = c(nlon, nlat, t2))\n",
    "ncatt_put(nc.ex , var1, '_FillValue', missval )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ca5987a-e841-4ee2-8546-b3e31ba3b9dc",
   "metadata": {},
   "source": [
    "#### Exercise 2.1\n",
    "\n",
    "Calculate the average temperature for the year 2012 and print the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "926929af-a65e-415e-8a1e-cc0ba6c4ea97",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "bc07b15f-25f9-48f3-859f-2fc4b71562e3",
   "metadata": {},
   "source": [
    "### Subset data for a specific region using latitude and longitude bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0af0de75-eb4a-4a54-a14a-3abc9be1e862",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read data but only for polygon\n",
    "min_lon <-  5. \n",
    "max_lon <- 16. \n",
    "# latitude is -90 to 90\n",
    "min_lat <- 45.\n",
    "max_lat <- 55.\n",
    "\n",
    "# find start indices for these longitudes and latitudes\n",
    "lon_index <- which(lon > min_lon & lon < max_lon)\n",
    "print(lon_index)\n",
    "\n",
    "lat_index <- which(lat > min_lat & lat < max_lat)\n",
    "print(lat_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "810f7fc0-6442-4d21-a16e-8b66d155a4f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the variable data for the selected latitudes and longitudes\n",
    "data_reg <- ncvar_get(ex.nc, variable_name, start = c(lon_index[1], lat_index[1], 1), count = c(length(lon_index), length(lat_index), -1))\n",
    "print(data_reg[,,1])\n",
    "\n",
    "# Close the NetCDF file\n",
    "nc_close(ex.nc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce594d94-6276-4379-8a9e-5b52e4b4a50a",
   "metadata": {},
   "source": [
    "#### Exercise 2.2\n",
    "\n",
    "Calculate and print the mean temperature over space and time for the selected region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a9ac6ae-95ff-4f7e-8ec9-e8f2a8203ded",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7238c2ca-c65b-4e84-b599-26b87e54b6a7",
   "metadata": {},
   "source": [
    "#### Exercise 2.3\n",
    "\n",
    "Calculate and print the mean temperature for the selected region for the full timeseries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b2cf70b-a63e-400f-87df-6b9b5b5d4b0c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "099df19f-fcd4-4653-8077-758a7810b53f",
   "metadata": {},
   "source": [
    "#### Exercise 2.4\n",
    "\n",
    "Calculate and print the mean temperature over space and time for another region of choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d310605-1f35-4a3f-80dd-2d6b0b59a3e4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "16c31e88-4a1e-46d4-b6a5-ad685c48a439",
   "metadata": {},
   "source": [
    "#### You can also use packages like \"raster\" to work with netcdf data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5bd6f93-0753-4e93-9fea-537f3e358eec",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# load necessary packages\n",
    "library(raster)\n",
    "\n",
    "# get working directory and load function eurocentric.r\n",
    "getwd()\n",
    "source(\"eurocentric.r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d658c25-fd30-41ec-9798-5b7ccc75fe0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "nc_file <- paste(indir, \"tas_Amon_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_201001-201412.nc\", sep=\"\")\n",
    "nc_data <- brick(nc_file) + 0\n",
    "print(nc_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9aaf4892-6221-4f07-bfc9-00ee4245d727",
   "metadata": {},
   "outputs": [],
   "source": [
    "nc_data = convert.to.eurocentric(nc_data)\n",
    "print(nc_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94b4675f-457e-473b-a44a-35aef5c0215f",
   "metadata": {},
   "source": [
    "❓ **Question**\n",
    "\n",
    "What does the function convert.to.eurocentric do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c601b5c4-3ff3-4342-be08-c1411d9665a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the time dimension\n",
    "nc <- nc_open(nc_file)\n",
    "\n",
    "time <- ncvar_get(nc, \"time\")\n",
    "time_units <- ncatt_get(nc, \"time\", \"units\")$value\n",
    "print(time_units)\n",
    "nc_close(nc)\n",
    "\n",
    "# extract base date, third bit of time_units \"days since 1950-01-01 00:00\"\n",
    "base_date <- as.character(lapply(strsplit(as.character(time_units),\n",
    "                                          split = \" \"), \"[\", 3))\n",
    "print(base_date)\n",
    "time_dates <- as.Date(time, format = \"%j\",\n",
    "                  origin = as.Date(base_date)) # define time as date\n",
    "print(time_dates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fffe53a-b41d-47dc-8233-1393285f16b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add time as date to rasterbrick\n",
    "raster_ts <- setZ(nc_data, time_dates)\n",
    "print(raster_ts)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Teaching (R SLAC 2021)",
   "language": "R",
   "name": "teach_r_slac_2021"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}