{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# readDiag\n",
    "\n",
    "In order to facilitate the access to the contents of the [Gridpoint Statistical Interpolation (GSI)](https://dtcenter.org/community-code/gridpoint-statistical-interpolation-gsi) diagnostic files, it has been written the [readDiag](https://github.com/GAD-DIMNT-CPTEC/readDiag) package wich is Fortran/Python interface that aims to read the GSI diagnostic binary files and create appropriate data structures to it. This Jupyter notebook is a simple walkthrough to its usage.\n",
    "\n",
    "## Usage\n",
    "\n",
    "In order to use the readDiag, first one need to load the `gsidiag` module from the readDiag package. This module can be loaded just like any other Python module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment next line to use in Google Colab\n",
    "#!pip install readDiag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import gsidiag as gd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the module loaded, open a file containing the GSI diagnostic information. In this example, the file `diag_conv_01.2020031112` (with ~64 Mb) is a diagnostic file provided by the GSI containing information from conventional observations used in the minimization of a 3DVar cost function for the date 2020031112. This file is provided to the user for testing purposes only.\n",
    "\n",
    "For this, we will create the object `file` from which we will use the builtin methods and functions of the readDiag.\n",
    "\n",
    "**Note:** the `%%time` is a Jupyter magic command that accounts for the time a command takes to run and it is not required for the use of the readDiag).\n",
    "\n",
    "To load the provided conventional diagnostic file, do as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "file = gd.read_diag('../data/diag_conv_01.2020031112')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## List of built-in methods and functions\n",
    "\n",
    "The readDiag package has several built-in methods and functions. To get a list of the available methods and functions, type `file.` and push the `<TAB>` button twice on the keyboard:\n",
    "\n",
    "```python\n",
    ">>> file.\n",
    "file.close(      file.obsInfo     file.pfileinfo(  file.pvmap(      file.zlevs       \n",
    "file.kxcount(    file.overview(   file.plot(       file.varNames    \n",
    "file.obs         file.pcount(     file.ptmap(      file.vcount(\n",
    "```\n",
    "\n",
    "The built-in methos and functions have documentation which can be accessed using Python's docstrings by issuing the following command:\n",
    "\n",
    "```python\n",
    "print(object.function_name.__doc__)\n",
    "```\n",
    "\n",
    "## Retrieve file information \n",
    "\n",
    "Now, let's get the basic information contained inside the GSI conventional diagnostic file. The information are related to the binary file itself as well as the observations and its related contents. Use the method `pfileinfo()` to retrieve what observations and types of observations (`kx`) are available inside the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(file.pfileinfo.__doc__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.pfileinfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to retrieve the file content information is through the method `overview()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(file.overview.__doc__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.overview()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another similar way to get a dictionary containing the basic information is by inspecting the attribute `obsInfo`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.obsInfo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the previous command returned a dictionary with all the information for all the available observations and types. To filter out these information by the observations, you may use the following commands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.obsInfo.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.obsInfo['ps']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the outcome of the command `file.obsInfo['ps']` returned a Pandas dataframe (acctually a multiindex dataframe) with all the information related with the `ps` (surface pressure) observation. More information on how to handle the dataframe and filter out some information are given on the next sections.\n",
    "\n",
    "Use the attribute `varNames` to get a list of the variable names from the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.varNames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to get the list of variables is by using some of the Pandas methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.obs.groupby(level=0).size().sort_values(ascending=False).keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a list of the available standard vertical levels, use the attribute `zlevs`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.zlevs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data exploration\n",
    "\n",
    "On the previous section it has been given an overview about the functions and methods of the readDiag package. In this section, it is given an overview of how to handle the Pandas dataframe in order to retrive and filter out some information related to the observations from the provided conventionan diagnostic file. \n",
    "\n",
    "To simply count the number of records of some kind of observation, use the function `pcount()`. This method returns a histogram with the count of kinds (i.e., `kx`) of a given observation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(gd.plot_diag.pcount.__doc__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gd.plot_diag.pcount(file,'ps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gd.plot_diag.pcount(file,'t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gd.plot_diag.pcount(file,'q')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gd.plot_diag.pcount(file,'uv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataframe exploration\n",
    "\n",
    "The Pandas dataframe is a tabular data structure with its own built in functions and methods. This section shows some generic operations related to it.\n",
    "\n",
    "Previously, we learned that the method `obsInfo` returns a dictionary with one DataFrame for each observation, e.g., for the `t` (temperature) observation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.obsInfo['t']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This dataframe is multiindex, i.e., a tabular data structure with multiple indices. The indices can be retrieve issuing the command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.obsInfo['t'].index.names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `kx` and `points` are the indices for whose all the columns of the dataframe are defined. This means that for a given observation, there are a number of `kx` (observation kinds) and `points` (the index of each point related to that observation kind).\n",
    "\n",
    "**Note:** although the `points` index each observation, its real world lat/lon coordinate is given by the column `geometry` (which is a tuple formed by the `lat` and `lon` columns).\n",
    "\n",
    "In order to simplify the use of some commands, we will create the variable `df_t` which will hold the temperature observation dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t = file.obsInfo['t']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The temperature observation has the following kinds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t.index.get_level_values('kx').unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to filter all the temperature observations by the kind 120, you may issue the command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t.iloc[df_t.index.get_level_values('kx') == 120]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the same way we did before, we will create another variable `df_t120` wich will hold the information of the temperature observation of the kind 120:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120 = df_t.iloc[df_t.index.get_level_values('kx') == 120]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to print the shape of the dataframe, you may issue:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the same way, to print the names of the columns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to print the head (i.e., the first few lines) of the dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's count the number of observations marked with the `iuse=1` and `iuse=-1` (i.e., observations marked as used and monitored, respectively). First, to show the unique values of the `iuse` mark, use the command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120['iuse'].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, the `dt_120_iuse` is a variable wich holds just the `iuse` column of the `df_t120` dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120_iuse = df_t120['iuse']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To count the number of temperature observations of the kind 120 marked with the flag `iuse=1`or `iused=-1`, use the pandas `count()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120_iuse.loc[df_t120['iuse'] == 1].count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120_iuse.loc[df_t120['iuse'] == -1].count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the `30572 + 184 = 30756` accounts for the total number of rows of the `df_t120` dataframe (which accounts for the total number of observations).\n",
    "\n",
    "If you want to plot a series of OMF (Observation Minus Forecast), which accounts for the innovation brought by the observations to the forecast during the assimilation process, its possible to filter out the `omf` column and retrieve just the observations marked with the flag `iuse = 1` (in this case the temperature observations of the kind 120):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120['omf'].loc[df_t120_iuse == 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a figure with the series of the values of the `omf` by index of the temperature observations of the kind 120, use the `plot()` method of the Pandas dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120['omf'].loc[df_t120_iuse == 1].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Altought in the figure we can see that some indices shows greater values of `omf` (e.g., greater than +/- 10), we cannot figure out where nor at what level those observations are. In the next sections, some built in plotting functions from the readDiag are introduced in order to facilitate this kind of analysis.\n",
    "\n",
    "An alternative way to use select the columns of the Pandas dataframe is given below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_t120.loc[:,'omf']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making spatial plots using the readDiag built-in functions\n",
    "\n",
    "Until here, all methods and functions from the readDiag package have been used to explore some information from the provided conventional GSI diagnostic file trough an appropriated data structure. The readDiag package ships with some built-in functions that leverage the its capability to explore data by creating spatial plots. \n",
    "\n",
    "These functions are:\n",
    "\n",
    "* `ptmap()`: this function plots the selected observation for the selected kinds;\n",
    "* `pvmap()`: this function plots the selected observations without specifying its kinds. It used the flag iuse instead;\n",
    "* `plot()`: this function makes a plot for the selected observation by using information of the following columns available within the dataframe.\n",
    "\n",
    "Use the built-in function `ptmap()` to plot the spatial distribution of some observation by type (`kx`). Remember to use the docstring to get an overview of its usage:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(gd.plot_diag.ptmap.__doc__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "gd.plot_diag.ptmap(file,'uv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the built-in function `ptmap()` uses the information about the observation and its kind, it is also possible to pass a list of `kx` related to that observation type (this also helps to unclutter the plot and takes significantly less time to get the plot the figure). Remember that you can use the built-in function `pfileinfo()` to retrieve the `kx` of each existing observation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "gd.plot_diag.ptmap(file,'uv', [253, 259, 243, 242])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another interesting built-in function is the function `pvmap()` which plots the spatial distribution of different observations considering all the kinds for a given observation. In this case, it is needed to specify the `iuse` flag of the observation to be plotted, i.e., `iuse = 1` for used and `iuse = -1` for monitored. The statuses of the observations must be passed as a list with at least one element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(gd.plot_diag.pvmap.__doc__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "gd.plot_diag.pvmap(file,['uv','ps','t','q'], mask='iuse==1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The built-in funtion `plot` plots the spatial distribution of variables using one of the columns of the dataframe. In the example below, we made a plot for the values (when the column `obs` is used) of the used surface pressure observations of the kind 187:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(gd.plot_diag.plot.__doc__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time \n",
    "\n",
    "gd.plot_diag.plot(file,'ps', 187, 'obs', mask='iuse==1',legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As final step, in order to release the memory used to hold the binary file, use the `close` method, which should return `0` when the file is correctly closed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}