PCP_Sesión_13.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "1752ce061957fd1ce93d2e4db9494b7b29da7935"
   },
   "source": [
    "<center>\n",
    "<img src=\"https://habrastorage.org/webt/jq/gc/r5/jqgcr5azmatev6zvgufn8xftehy.png\">\n",
    "    \n",
    "##  II Escuela de Verano en Física Computacional\n",
    "\n",
    "Autor: [Gerardo Rivera](https://github.com/DangoMelon).\\\n",
    "Institución: Instituto Geofísico del Perú / UNMSM\\\n",
    "Correo: gerardo.rivera1@unmsm.edu.pe\n",
    "    \n",
    "Este material está sujeto a los términos y condiciones de la licencia [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/). Se permite el uso gratuito para cualquier propósito no comercial.\n",
    "\n",
    "*También puede consultar la última versión de este notebook en nuestro [repositorio](https://github.com/PCPUNMSM) y los videos de clase [en nuestro canal de Youtube](https://www.youtube.com/channel/UCze8vzQLcplutz0nWDNjFCA).*\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_uuid": "9b690f2c8129ad9445c544337375b93246c26668"
   },
   "source": [
    "# <center> Sesión 13. Manejo de datos n-dimensionales en python: *conociendo xarray*  \n",
    "\n",
    "<br>\n",
    "<div align=\"center\">\n",
    "    <img src=\"http://xarray.pydata.org/en/stable/_static/dataset-diagram-logo.png\" width=30%>\n",
    "</div>\n",
    "<br>\n",
    "\n",
    "## Outline\n",
    "1. [Conociendo xarray](#1.-Conociendo-xarray)\n",
    "2. [Funciones Básicas](#2.-Funciones-básicas)\n",
    "3. [Trabajando con datos](#3.-Trabajando-con-datos)\n",
    "4. [Gráficos](#4.-Gráficos)\n",
    "5. [Integración con dask](#5.-Integración-con-dask)\n",
    "6. [Aplicación con datos reales](#6.-Aplicación-con-datos-reales)\n",
    "7. [Aplicación Numérica](#7.-Aplicación-Numérica)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Conociendo xarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para trabajar con data cubes usando python haremos uso extensivo de los paquetes estrella de la plataforma Pangeo. Pangeo es un esfuerzo de la comunidad que promueve la reproductibilidad, escalamiento y transparencia en la ciencia.\n",
    "\n",
    "El paquete de cabecera para el manejo de dichas estructuras de datos será xarray. Xarray extiende las capacidades de numpy en el manejo de datos n-dimensionales introduciendo etiquetas en forma de dimensiones, coordenadas y atributos que facilitan las operaciones sobre los mismos. El motor para las operaciones usando las etiquetas es proporcionado por pandas; la capacidad de manejar datos n-dimensionales con facilidad es propocionado por numpy; y la escalabilidad de los datos es proporcionado por dask, convirtiendo a este paquete en una herramienta poderosa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quitar el comentario y correr esta línea si se esta usando google colab\n",
    "#!pip install -q -U git+https://github.com/pydata/xarray.git#egg=xarray[complete]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1. Estructuras de xarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xarray` cuenta con dos estructuras simples que son fundamentales comprender para un mejor uso del paquete. Esta es una brevisima descripción sobre algunos conceptos claves que se deben tener presentes, mayor informacion la pueden encontrar siempre en la [documentacion](http://xarray.pydata.org/en/stable/data-structures.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1.1. `DataArray`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un `DataArray` es virtualmente igual a un arreglo n-dimensional de numpy que cuenta con etiquetas en sus dimensiones, siendo este la estructura de datos base de xarray. Para denifinir un DataArray necesitaremos proporcionar, como mínimo, una arreglo de elementos, como los generados por `numpy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr.DataArray(np.random.randn(10, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El `repr` en html de xarray es bien practico cuando estamos trabajando en notebooks, sin embargo, el `repr` en puro texto contiene un poco de información adicional sobre el estado de las coordenadas en el `DataArray` recién creado"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with xr.set_options(display_style=\"text\"):\n",
    "    display(xr.DataArray(np.random.randn(10, 5)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De manera automática, nuestras dimensiones fueron asignadas con nombres genéricos ya que no agregamos ningun detalle adicional. Adicionalmente, en el `repr` de texto podemos observar un mensaje adicional sobre `Dimensions without coordinates`, esto se debe a que no tenemos ninguna coordenada definida para nuestras dimensiones. Las coordenadas pueden ser entendidas como etiquetas a cada elemento posicional en una dimensión.  \n",
    "\n",
    "Ahora crearemos un `DataArray` con la información que hace falta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creamos las coordenadas para nuestras dimensiones\n",
    "lat = np.arange(-90, 90, 0.25)\n",
    "lon = np.arange(0, 360, 0.25)\n",
    "\n",
    "# Creamos un poco de data falsa\n",
    "llon, llat = np.meshgrid(lon, lat)\n",
    "r = np.sqrt(llon ** 2 + llat ** 2)\n",
    "data = np.sin(r)\n",
    "\n",
    "# Creamos el DataArray\n",
    "xarr = xr.DataArray(data, coords=[lat, lon], dims=[\"lat\", \"lon\"], name=\"false_data\")\n",
    "xarr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nuestro `DataArray` ha sido creado dentro de la variable xarr con las dimensiones y coordenadas asignadas correctamente. Ahora podemos acceder a los métodos que xarray ofrece para sus objetos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xarr.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Existen varias formas convenientes de crear un `DataArray` que son practicas en algunos casos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr.zeros_like(xarr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Ejercicio**\n",
    "\n",
    "Crear un `DataArray` llamado \"temperatura\" a partir de datos aleatorios.\n",
    "\n",
    "- Debe tener tres coordenadas con nombres \"latitud\", \"longitud\" y \"profundidad\"\n",
    "- Debe tener atributos en cada coordenada asi como atributos generales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xr.DataArray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1.2. `Dataset`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un `Dataset` es la representación en la memoria del sistema de un archivo netcdf. Al igual que en pandas un `DataFrame` es un conjunto de `Series`, en xarray un `Dataset` es un conjunto de `DataArray`'s. Su declaración es un poco más extensa que un `DataArray` pero a cambio proporciona un contenedor a varias variables.\n",
    "\n",
    "<div align=\"center\">\n",
    "<img src=\"http://xarray.pydata.org/en/stable/_images/dataset-diagram.png\" width=80%>\n",
    "<br>\n",
    "<a href=\"http://xarray.pydata.org/en/stable/data-structures.html#dataset\">Documentación de Dataset</a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definimos las coordenadas para nuestras dimensiones\n",
    "# en este caso tendremos 3 dimensiones: latitud, longitud y tiempo\n",
    "lat = np.arange(-15, 15, 0.25)\n",
    "lon = np.arange(140, 280, 0.25)\n",
    "time = pd.date_range(\"2019-01-01\", freq=\"D\", periods=120)\n",
    "\n",
    "# Datos aleatorios\n",
    "temp = np.random.rand(120, len(lat), len(lon))\n",
    "precip = np.random.rand(120, len(lat), len(lon))\n",
    "elev = np.random.rand(len(lat), len(lon))\n",
    "serie = np.random.rand(len(time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Declaramos el Dataset\n",
    "ds = xr.Dataset(\n",
    "    {\n",
    "        \"temp\": ([\"time\", \"lat\", \"lon\"], temp),\n",
    "        \"precip\": ([\"time\", \"lat\", \"lon\"], precip),\n",
    "        \"elev\": ([\"lat\", \"lon\"], elev),\n",
    "        \"serie\": ([\"time\"], serie),\n",
    "    },\n",
    "    coords={\"time\": time, \"lat\": lat, \"lon\": lon},\n",
    ")\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si la terminología es un poco dificil de entender, siempre pueden consultar la [documentación](http://xarray.pydata.org/en/stable/terminology.html0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2. Seleccion de datos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xarray` cuenta con dos métodos para selección de datos: `sel` y `isel`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sel` funciona usando los nombres (o etiquetas) asignados a cada dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.tutorial.load_dataset(\"air_temperature\")\n",
    "air = ds.air"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.sel(time=\"2013-04-01\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usando `slice` podemos seleccionar un rango de valores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.sel(time=slice(\"2013-04-01\", \"2013-09-15\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tampoco es necesario ser precisos con el valor de las coordenadas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.sel(lat=70.5, lon=211, method=\"nearest\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`isel` funciona usando los indices posicionales, pero aun usando los nombres de las dimensiones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.isel(time=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.isel(lat=slice(10, 30))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.isel(lat=10, lon=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esta forma de seleccionar los datos, usando el nombre de la dimensión, es de lejos mucho mejor que hacer indexing de la manera clásica\n",
    "\n",
    "```python\n",
    "air[10:40,10:40,14]\n",
    "```\n",
    "\n",
    "A simple vista uno no puede saber que valores esta seleccionando"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Funciones básicas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Muchas de las funciones disponibles en `xarray` son tomadas de `numpy` y `pandas`, presentando el mismo comportamiento en muchos casos.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El `Dataset` que hemos abierto tiene una sola variable: `air`. Para extraerla, podemos usar algunas notaciones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.air"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds[\"air\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si bien ambas notaciones son equivalente, la notacion \"punto\", i.e. `ds.air`, no funcionara si la variable tiene el nombre de alguna de las funciones incorporadas en python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air = ds[\"air\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Operaciones básicas de agregación disponibles en numpy estan disponibles en `xarray`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sin embargo, esto no impide que usemos operaciones de `numpy` sobre estructuras de `xarray`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(air)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Trabajando con datos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. Broadcasting y Alignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uno de los conceptos importantes al trabajar con arreglos de `numpy` es la propagacion de dimensiones o *broadcasting*.\n",
    "\n",
    "\n",
    "<div align=\"center\">\n",
    "<img src=\"http://scipy-lectures.org/_images/numpy_broadcasting.png\" width=50%>\n",
    "</div>\n",
    "\n",
    "<a style=\"font-size: 0.8em\" href=\"http://scipy-lectures.org/intro/numpy/operations.html#broadcasting\">Scipy Lecture Notes</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = xr.DataArray(np.arange(3), dims=\"time\", coords={\"time\": np.arange(3)})\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = xr.DataArray(np.arange(4), dims=\"space\", coords={\"space\": np.arange(4)})\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Por otro lado, la lineacion de coordenadas o *alignment* es una operacion que xarray realiza cuando encuentra coordenadas en comun entre dos o más `DataArray`\n",
    "\n",
    "<div align=\"center\">\n",
    "<img src=\"https://xarray-contrib.github.io/xarray-tutorial/_images/broadcasting.png\" width=50%>\n",
    "</div>\n",
    "\n",
    "<a style=\"font-size: 0.8em\" href=\"https://docs.google.com/presentation/d/16CMY3g_OYr6fQplUZIDqVtG-SKZqsG8Ckwoj2oOqepU/edit#slide=id.g67406e8d7_0_135\">Imagen de Stephan Hoyer</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = air.sel(time=slice(\"2013-01-01\", \"2013-12-31\")).isel(lat=0)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = air.sel(time=slice(\"2013-10-01\", \"2014-12-31\")).isel(lat=10)\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a + b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Este tipo de operaciones permite mantener los datos en la coordenada que les corresponde, mantiendo todo siempre en la misma grilla."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Operaciones de alto nivel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tomado de [xarray in 45 minutes](https://xarray-contrib.github.io/xarray-tutorial/oceanhackweek-2020/xarray-oceanhackweek20.html#High-level-computation:-groupby,-resample,-rolling,-coarsen,-weighted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xarray` cuenta con algunos objetos de alto nivel que permiten realizar operaciones comunes\n",
    "\n",
    "- groupby : Agrupa y reduce datos\n",
    "- resample : Groupby especializado para ejes de tiempo. Puede aumentar o reducir la resolución temporal.\n",
    "- rolling : Aplica ventanas moviles sobre los datos e.g. media corrida\n",
    "- coarsen : Reduce la resolución de los datos\n",
    "- weighted : Aplica pesos previo a un operación de reducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`groupby` es altamente útil en la dimensión del tiempo, especialmente para calcular climatologías"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.groupby(\"time.day\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.groupby(\"time.month\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.groupby(\"time.month\").mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`resample` nos permite cambiar la frecuencia en la dimensión del tiempo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.resample(time=\"1D\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.resample(time=\"1D\").mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`rolling` permite construir operaciones sobre ventanas móviles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air_ts = air.isel(lat=10, lon=10)\n",
    "air_ts.rolling(time=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air_ts.plot(label=\"orig\")\n",
    "air_ts.rolling(time=30).mean().plot(label=\"rolling\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. Otras operaciones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xarray` provee muchas operaciones de conveniencia. Siempre es una buena idea revisar la documentación para ver todo lo que tiene por ofrecer. [Link](http://xarray.pydata.org/en/stable/api.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.differentiate(\"time\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.integrate(\"time\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.isel(lat=10, lon=10).polyfit(\"time\", deg=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Gráficos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para una descripción más detallada, revisar https://xarray.pydata.org/en/stable/plotting.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xarray` permite realizar plots simples de nuestros datos usando el método `.plot`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Es posible hacer una grilla de gráficos especificando algunos parámetros"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.isel(time=slice(100, 105)).plot(col=\"time\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air.isel(time=slice(100, 105)).plot.contourf(col=\"time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Integración con dask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`xarray` tiene la capacidad de leer multiples archivos netcdf dentro de un solo objeto mediante la función `open_mfdataset`. Ya que no contamos con una gran cantidad de datos por el momento, vamos a simular este procedimiento aplicando `chunks` a nuestros datos de temperatura"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air_dask = air.chunk(dict(time=4))\n",
    "air_dask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Las operaciones a realizar sobre esta nueva variable seran evaluadas de manera \"vaga\" (*lazy*), esto quiere decir que ninguna operación será realizada hasta que se llame `compute` o `load`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air_dask.mean(dim=\"lat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air_dask.groupby(\"time.month\").mean(dim=\"time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Aplicación con datos reales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para cargar datos en formato netcdf a una variable, se hacer uso de la función open_dataset (en el caso de tener 1 solo archivo) o open_mfdataset en caso de tener muchos archivos. Es este caso vamos a leer datos de temperatura superficial del mar de ERSSTv5 para hacer el cálculo del Oceanic Niño Index (ONI)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ersstv5 = xr.open_dataset(\n",
    "    \"https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/noaa.ersst.v5/sst.mnmean.nc\"\n",
    ")\n",
    "ersstv5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculo del ONI\n",
    "\n",
    "El ONI esta definido como la media corrida de 3 meses sobre el indice del Niño 3.4 (5°S-5°N / 170°W-120°W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nino34 = ersstv5.sst.sel(lat=slice(5, -5), lon=slice(190, 240))\n",
    "nino34"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nino34 = nino34.mean(dim=[\"lat\", \"lon\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cálculo de la climatología\n",
    "\n",
    "Para realizar este cálculo, vamos a tomar como periodo base 1986-2015 con fines ilustrativos. En el cálculo actual, el periodo de la climatología cambia cada 5 años segun los criterios considerados por el centro de predicción del clima de la NOAA ([ref](https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_change.shtml))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nino34_clim = (\n",
    "    nino34.sel(time=slice(\"1986-01-01\", \"2015-12-31\"))\n",
    "    .groupby(\"time.month\")\n",
    "    .mean(dim=\"time\")\n",
    ")\n",
    "nino34_clim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "nino34_clim.plot(ax=ax)\n",
    "ax.set_yticks(np.arange(26, 28.1, 0.2))\n",
    "ax.set_ylim(26, 28)\n",
    "ax.grid(ls=\"--\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cálculo de la anomalía\n",
    "\n",
    "Utilizaremos una sintaxis similar para restar las climatologías a los valores absolutos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nino34_anom = nino34.groupby(\"time.month\") - nino34_clim\n",
    "nino34_anom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ONI = nino34_anom.rolling(time=3, center=True).mean()\n",
    "\n",
    "# truncamos los valores al segundo decimal\n",
    "ONI = np.around(ONI, decimals=2)\n",
    "ONI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ONI.attrs[\"long_name\"] = \"SST Anomaly\"\n",
    "ONI.attrs[\"units\"] = \"degrees C\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ONI.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Con algo de tiempo, podemos hacer un gráfico un poco más atractivo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Colocamos el tiempo que queremos plotear\n",
    "plot_data = ONI.sel(time=slice(\"2008-01-01\", \"2019-12-31\"))\n",
    "\n",
    "# Declaramos la figura y los ejes\n",
    "fig, ax = plt.subplots(dpi=200)\n",
    "\n",
    "ax.fill_between(\n",
    "    plot_data.time.data,\n",
    "    plot_data.data,\n",
    "    0.5,\n",
    "    where=plot_data > 0.5,\n",
    "    color=\"red\",\n",
    "    interpolate=True,\n",
    "    lw=0,\n",
    "    label=\"ONI > 0.5\",\n",
    ")\n",
    "ax.fill_between(\n",
    "    plot_data.time.data,\n",
    "    plot_data.data,\n",
    "    -0.5,\n",
    "    where=plot_data < -0.5,\n",
    "    color=\"blue\",\n",
    "    interpolate=True,\n",
    "    lw=0,\n",
    "    label=\"ONI < 0.5\",\n",
    ")\n",
    "plot_data.plot(ax=ax, c=\"k\", lw=0.5)\n",
    "\n",
    "minx = plot_data.time.min().data\n",
    "maxx = plot_data.time.max().data\n",
    "ax.hlines(0.5, minx, maxx, lw=0.5, linestyles=\"--\")\n",
    "ax.hlines(-0.5, minx, maxx, lw=0.5, linestyles=\"--\")\n",
    "\n",
    "ax.set_ylim(-4, 4)\n",
    "ax.set_xlim(minx, maxx)\n",
    "\n",
    "ax.text(\n",
    "    0.99,\n",
    "    0.01,\n",
    "    \"Clim 1981-2010\",\n",
    "    fontsize=5,\n",
    "    horizontalalignment=\"right\",\n",
    "    transform=ax.transAxes,\n",
    ")\n",
    "\n",
    "ax.set_title(\"ERSSTv5 ONI\")\n",
    "ax.grid(ls=\"--\", lw=1, alpha=0.3)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Aplicación Numérica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Tomado de un ejercicio del curso \"Modelado Numérico de la Atmósfera\" organizado por el IGP [[Repo]](https://github.com/DangoMelon/cmna-igp)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dada la siguiente malla de datos de temperatura, pronosticar la evolución del campo de la variable para una hora. Aplicar el método de Euler para las derivadas zonal y meridional.\n",
    "\n",
    "### Datos\n",
    "* Velocidad del viento = CONST = 7m/s\n",
    "* Dirección del viento = CONST = 140°\n",
    "* Resolución del dominio dx = dy = 18km\n",
    "* Paso del tiempo -> por determinar\n",
    "\n",
    "$$\n",
    "TEMP = \\begin{bmatrix}\n",
    "21.2 & 22.8 & 23.5 & 24.8 & 26.5 \\\\\n",
    "22.7 & 24.4 & 25.2 & 26.5 & 28.4 \\\\\n",
    "24.3 & 26.1 & 27.0 & 28.4 & 30.4 \\\\\n",
    "26.0 & 27.9 & 28.9 & 30.4 & 32.5 \\\\\n",
    "27.8 & 29.9 & 30.9 & 32.5 & 34.8 \\\\\n",
    "29.7 & 32.0 & 33.0 & 34.8 & 37.2 \\\\\n",
    "31.8 & 34.2 & 35.3 & 37.2 & 39.8\n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solución"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La ecuación que estamos tratando de resolver es la siguiente:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial T }{\\partial t} + u \\frac{\\partial T}{\\partial x} + v \\frac{\\partial T}{\\partial y} = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En donde el paso de tiempo estara determinado por el numero de Courant\n",
    "\n",
    "$$\n",
    "\\frac{\\lvert u\\rvert*\\Delta t}{\\Delta x} + \\frac{\\lvert v\\rvert*\\Delta t}{\\Delta y} \\leq 0.06\n",
    "\\\\\n",
    "$$\n",
    "\n",
    "Que resulta en\n",
    "\n",
    "$$\n",
    "\\Delta t \\leq \\frac{0.06}{\\left(\\frac{\\lvert u\\rvert}{\\Delta x} + \\frac{\\lvert v\\rvert}{\\Delta y}\\right)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Graficos interactivos\n",
    "import holoviews as hv\n",
    "import panel as pn\n",
    "from panel.interact import fixed, interact\n",
    "\n",
    "pn.extension()\n",
    "hv.extension(\"bokeh\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Declaración de datos\n",
    "dx = dy = 18\n",
    "magw = 7\n",
    "angw = 140\n",
    "temp_arr = np.array(\n",
    "    [\n",
    "        [21.2, 22.8, 23.5, 24.8, 26.5],\n",
    "        [22.7, 24.4, 25.2, 26.5, 28.4],\n",
    "        [24.3, 26.1, 27.0, 28.4, 30.4],\n",
    "        [26.0, 27.9, 28.9, 30.4, 32.5],\n",
    "        [27.8, 29.9, 30.9, 32.5, 34.8],\n",
    "        [29.7, 32.0, 33.0, 34.8, 37.2],\n",
    "        [31.8, 34.2, 35.3, 37.2, 39.8],\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculo de las componentes u,v\n",
    "u = magw * np.cos(np.deg2rad(140))\n",
    "v = magw * np.sin(np.deg2rad(140))\n",
    "\n",
    "print(f\"u = {u}\\n v = {v}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Con los valores obtenidos para u y v, calculamos el valor límite que debe tener $\\Delta t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt_approx = 0.06 / ((np.abs(u) / (dx * 1000)) + (np.abs(v) / (dy * 1000)))\n",
    "print(f\"dt_approx = {dt_approx}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entonces tomamos\n",
    "\n",
    "$$\n",
    "\\Delta t \\approx 100s\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construimos una grilla que cumpla con $\\Delta x$ y $\\Delta y$ (no es necesario realizar esto, pero es util para tener coordenadas definidas dentro de xarray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon = np.arange(-90, -90 + temp_arr.shape[1] * dx, dx)\n",
    "lat = np.arange(-50, -50 + temp_arr.shape[0] * dy, dy)\n",
    "print(f\"lon = {lon}\\nlat = {lat}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construimos nuestra matriz de datos con `lat` y `lon` como coordenadas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp = xr.DataArray(temp_arr, coords=[(\"lat\", lat[::-1]), (\"lon\", lon)])\n",
    "\n",
    "# Metadata adicional para variables/coordenadas que sigue la convención CF\n",
    "temp.lon.attrs[\"long_name\"] = \"Longitude\"\n",
    "temp.lon.attrs[\"units\"] = \"degrees_east\"\n",
    "temp.lat.attrs[\"long_name\"] = \"Latitude\"\n",
    "temp.lat.attrs[\"units\"] = \"degrees_north\"\n",
    "temp.attrs[\"long_name\"] = \"Temperature\"\n",
    "temp.attrs[\"units\"] = \"degrees_C\"\n",
    "temp.name = \"temperature\"\n",
    "\n",
    "temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Arreglos para grafico de quivers\n",
    "u_arr = np.full_like(temp_arr, u)\n",
    "v_arr = np.full_like(temp_arr, v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se puede explorar nuestra malla de datos iniciales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(dpi=100)\n",
    "temp.plot(ax=ax, vmin=20, vmax=40)\n",
    "ax.quiver(lon, lat, u_arr, v_arr)\n",
    "ax.set_title(\"En grilla\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "levels = np.arange(20, 40.5, 0.5)\n",
    "\n",
    "fix, ax = plt.subplots(dpi=100)\n",
    "temp.plot.contourf(ax=ax, vmin=20, vmax=40, levels=levels)\n",
    "CS = temp.plot.contour(ax=ax, vmin=20, vmax=40, colors=\"k\", levels=levels[::4])\n",
    "ax.quiver(lon + dx / 2, lat + dy / 2, u_arr, v_arr)\n",
    "ax.set_ylim(lat.min(), lat.max())\n",
    "ax.set_xlim(lon.min(), lon.max())\n",
    "ax.clabel(CS, colors=\"k\", fmt=\"%i\")\n",
    "ax.set_title(\"En contornos\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La discretización de la ecuacion a resolver se hara con un esquema centrado en el espacio y adelantado en el tiempo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\\frac{T_{i,j}^{n+1} - T_{i,j}^{n}}{\\Delta t}\\right) + u \\left(\\frac{T_{i+1,j}^{n} - T_{i-1,j}^{n}}{2\\Delta x} \\right) + v \\left( \\frac{T_{i,j+1}^{n} - T_{i,j-1}^{n}}{2\\Delta y} \\right) = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De donde"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "T_{i,j}^{n+1} = \\Delta t \\ T_{i,j}^{n} - \\frac{u \\Delta t}{2 \\Delta x} \\left( T_{i+1,j}^{n} - T_{i-1,j}^{n} \\right) - \\frac{v \\Delta t}{2 \\Delta y} \\left( T_{i,j+1}^{n} - T_{i,j-1}^{n} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implementamos la discretizacion obtenida en una función que nos permitirá obtener la temperatura pronosticada"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve(\n",
    "    inarr,\n",
    "    dt,\n",
    "    u,\n",
    "    v,\n",
    "    nt,\n",
    "    boundary=\"1st_order\",\n",
    "    tstart=\"2021-04-03T00:00\",\n",
    "    plot=False,\n",
    "    dynamic=False,\n",
    "    **hv_kwargs,\n",
    "):\n",
    "    \"\"\"\n",
    "    Resuelve la ecuacion de momento en dos dimensiones con velocidades constantes\n",
    "\n",
    "    Parametros\n",
    "    ----------\n",
    "    inarr: DataArray\n",
    "        Objeto de xarray con grilla regular (lat, lon) que\n",
    "        representa el campo inicial\n",
    "    dt: int\n",
    "        Salto de tiempo\n",
    "    u: float\n",
    "        Valor de la velocidad zonal\n",
    "    v: float\n",
    "        Valor de la velocidad meridional\n",
    "    nt: int\n",
    "        Saltos de tiempo a resolver\n",
    "    boundary: {'1st_order', '2nd_order', 'constant'}\n",
    "        Condicion de frontera a usar.\n",
    "        Por defecto se usa '1st_order'\n",
    "        Puede ser una de las siguientes opciones:\n",
    "            -'1st_order': aproximacion de primer orden hacia\n",
    "                adelante/atras en los bordes\n",
    "            -'2nd_order': aproximacion de segundo orden hacia\n",
    "                adelante/atras en los bordes\n",
    "            -'constant': frontera constante, i.e dT/dr = 0\n",
    "    tstart: string, opcional\n",
    "        Fecha inicial de la simulación. Útil para la dimension\n",
    "        tiempo a agregar en el array de salida.\n",
    "        Por defecto toma el valor 2021-04-03T00:00\n",
    "\n",
    "    Retorna\n",
    "    -------\n",
    "    outarr: DataArray\n",
    "        Arreglo con las mismas coordenadas que inarr pero con\n",
    "        la dimension `time` agregada, de longitud `nt`+1 y\n",
    "        saltos `dt`.\n",
    "    \"\"\"\n",
    "\n",
    "    boundaries = {\"1st_order\": 1, \"2nd_order\": 2}\n",
    "    diff_kwargs = {}\n",
    "\n",
    "    tdelta = pd.Timedelta(dt, unit=\"s\")\n",
    "    tstart = pd.to_datetime(tstart)\n",
    "    outarr = xr.DataArray(\n",
    "        np.nan,\n",
    "        coords=[\n",
    "            (\"time\", pd.date_range(tstart, tstart + (tdelta * nt), freq=f\"{dt}s\")),\n",
    "            (\"lat\", inarr.lat),\n",
    "            (\"lon\", inarr.lon),\n",
    "        ],\n",
    "    )\n",
    "    outarr[0] = inarr\n",
    "\n",
    "    const = 1\n",
    "    if boundary in boundaries.keys():\n",
    "        diff_kwargs = dict(edge_order=boundaries[boundary])\n",
    "    elif boundary == \"constant\":\n",
    "        const = xr.zeros_like(inarr)\n",
    "        const[1:-1, 1:-1] = 1\n",
    "\n",
    "    for i in range(outarr.time.size - 1):\n",
    "        current_time = outarr.isel(time=i)\n",
    "        cdiffx = current_time.differentiate(\"lon\", **diff_kwargs) * const / 1000\n",
    "        cdiffy = current_time.differentiate(\"lat\", **diff_kwargs) * const / 1000\n",
    "        outarr[i + 1] = (current_time) - (u * dt * cdiffx) - (v * dt * cdiffy)\n",
    "\n",
    "    # Colocamos la metadata\n",
    "    outarr[\"boundary\"] = boundary\n",
    "    outarr = outarr.expand_dims(\"boundary\")\n",
    "    outarr.lat.attrs = inarr.lat.attrs\n",
    "    outarr.lon.attrs = inarr.lon.attrs\n",
    "    outarr.time.attrs[\"long_name\"] = \"Time\"\n",
    "    outarr.boundary.attrs[\"long_name\"] = \"Boundary Condition\"\n",
    "    outarr.attrs = inarr.attrs\n",
    "    outarr.name = \"forecast_temp\"\n",
    "\n",
    "    # En caso de querer retornar el gráfico\n",
    "    if plot:\n",
    "        hmap = hv.Dataset(outarr).to(hv.Image, kdims=[\"lon\", \"lat\"])\n",
    "        cont = hv.operation.contours(\n",
    "            hmap, filled=True, dynamic=dynamic, levels=np.arange(-5, 60, 1.0).tolist()\n",
    "        ).redim.range(forecast_temp=(20, 40))\n",
    "        xx, yy = np.meshgrid(\n",
    "            np.linspace(lon.min(), lon.max(), 30), np.linspace(lat.min(), lat.max(), 30)\n",
    "        )\n",
    "        vector = hv.VectorField(\n",
    "            (\n",
    "                xx,\n",
    "                yy,\n",
    "                np.full_like(xx, np.arctan2(v, u)),\n",
    "                np.full_like(xx, 1),\n",
    "            )\n",
    "        )\n",
    "        layout = cont * vector\n",
    "        layout.opts(\n",
    "            hv.opts.Polygons(\n",
    "                tools=[\"hover\", \"save\"],\n",
    "                default_tools=[],\n",
    "                ylim=(lat.min(), lat.max()),\n",
    "                xlim=(lon.min(), lon.max()),\n",
    "                frame_width=500,\n",
    "                aspect=\"equal\",\n",
    "                colorbar=True,\n",
    "                color_levels=22,\n",
    "            ),\n",
    "            hv.opts.Overlay(default_tools=[]),\n",
    "            hv.opts.VectorField(default_tools=[]),\n",
    "        )\n",
    "        if dynamic:\n",
    "            pan = pn.panel(layout, widget_type=\"scrubber\", widget_location=\"bottom\")\n",
    "            # Modificamos la velocidad de reproducción por defecto\n",
    "            player = pan[1][1][0]\n",
    "            player.param.interval.default = 150\n",
    "            return pan\n",
    "        return layout\n",
    "\n",
    "    return outarr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fore_temp = solve(temp, 100, u, v, 36, boundary=\"2nd_order\")\n",
    "fore_temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exploramos el resultado para condiciones de frontera con solución de segundo orden y constante"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcst_2nd = solve(temp, 100, u, v, 36, plot=True, boundary=\"2nd_order\")\n",
    "fcst_2nd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcst_const = solve(temp, 100, u, v, 36, plot=True, boundary=\"constant\")\n",
    "fcst_const"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Live Demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "live_panel = interact(\n",
    "    solve,\n",
    "    inarr=fixed(temp),\n",
    "    dt=(10, 110, 10, 100),\n",
    "    u=(-15, 15, 0.1, u),\n",
    "    v=(-15, 15, 0.1, v),\n",
    "    nt=(18, 36 * 5, 18, 36),\n",
    "    plot=fixed(True),\n",
    "    boundary=[\"1st_order\", \"2nd_order\", \"constant\"],\n",
    "    dynamic=fixed(True),\n",
    ")\n",
    "\n",
    "pn.Row(\n",
    "    pn.layout.HSpacer(),\n",
    "    pn.Column(\n",
    "        pn.pane.PNG(\n",
    "            \"https://repositorio.igp.gob.pe/themes/Mirage/images/logo_igp_sharecontent.png\",\n",
    "            align=\"center\",\n",
    "            width=250,\n",
    "            margin=-25,\n",
    "        ),\n",
    "        pn.pane.Markdown(\"## Difusión de la Temperatura\\n **Parámetros**\", margin=0),\n",
    "        live_panel[0],\n",
    "        pn.layout.VSpacer(),\n",
    "        pn.pane.Markdown(\"*Club de Python para Físicos*\"),\n",
    "        margin=20,\n",
    "    ),\n",
    "    live_panel[1],\n",
    "    pn.layout.HSpacer(),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 8. Recursos útiles\n",
    "\n",
    "* [xarray Tutorial](https://xarray-contrib.github.io/xarray-tutorial/index.html)\n",
    "* [Pangeo](https://pangeo.io/)\n",
    "* [Scipy Lecture Notes](http://scipy-lectures.org/index.html)\n",
    "* [\"Scaling Scientific Python\"](https://matthewrocklin.com/slides/ams-2018.html#/) por Matthew Rocklin\n",
    "* [Documentación](http://xarray.pydata.org/en/stable/index.html) oficial de xarray\n",
    "* [\"multi-dimensional data analysis in Python\"](https://fabienmaussion.info/acinn_xarray_workshop/#/) por Fabien Maussion\n",
    "* [An Introduction to Earth and Environmental Data Science](https://earth-env-data-science.github.io/intro)\n",
    "## Donaciones\n",
    "<br>\n",
    "<center>\n",
    "Puedes donar una vez en el siguiente enlace (Ko-Fi):\n",
    "\n",
    "<br>\n",
    "    \n",
    "*Click en la imagen.*\n",
    "\n",
    "<a href=\"https://ko-fi.com/rcrdphysics\">\n",
    "<img src=\"https://habrastorage.org/webt/8r/ml/xf/8rmlxfpdzukegpxa62cxlfvgkqe.png\" width=20% />\n",
    "\n",
    "</center>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}