ECCO_data_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ECCOv4 Data Analysis Example\n",
    "\n",
    "Largely inspired by the notebook of the same name from the [Pangeo Example Gallery](http://gallery.pangeo.io/repos/pangeo-gallery/physical-oceanography/04_eccov4.html).\n",
    "\n",
    "This Jupyter notebook demonstrates how to use [xarray](http://xarray.pydata.org/en/latest/) and [xgcm](http://xgcm.readthedocs.org) to analyze data from the [ECCO v4r3](https://ecco.jpl.nasa.gov/products/latest/) ocean state estimate.\n",
    "\n",
    "First we import our standard python packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data\n",
    "\n",
    "The ECCOv4r3 data was converted from its raw MDS (.data / .meta file) format to zarr format, using the [xmitgcm](http://xmitgcm.readthedocs.io) package. [Zarr](http://zarr.readthedocs.io) is a powerful data storage format that can be thought of as an alternative to HDF. In contrast to HDF, zarr works very well with cloud object storage. Zarr is currently useable in python, java, C++, and julia. It is likely that zarr will form the basis of the next major version of the netCDF library.\n",
    "\n",
    "If you're curious, here are some resources to learn more about zarr:\n",
    "  - https://zarr.readthedocs.io/en/stable/tutorial.html\n",
    "  - https://speakerdeck.com/rabernat/pangeo-zarr-cloud-data-storage\n",
    "  - https://mrocklin.github.com/blog/work/2018/02/06/hdf-in-the-cloud\n",
    "\n",
    "The ECCO zarr data currently lives in [Google Cloud Storage](https://cloud.google.com/storage/) as part of the [Pangeo Data Catalog](http://catalog.pangeo.io/). This means we can open the whole dataset using one line of code.\n",
    "\n",
    "This takes a bit of time to run because the metadata must be downloaded and parsed. The type of object returned is an [Xarray dataset](http://xarray.pydata.org/en/latest/data-structures.html#dataset)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import intake\n",
    "cat = intake.open_catalog(\"https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/ocean.yaml\")\n",
    "ds = cat.ECCOv4r3.to_dask()\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that no data has actually been downloaded yet. Xarray uses the approach of _lazy evaluation_, in which loading of data and execution of computations is delayed as long as possible (i.e. until data is actually needed for a plot). The data are represented symbolically as [dask arrays](http://docs.dask.org/en/latest/array.html). For example:\n",
    "\n",
    "    SALT       (time, k, face, j, i) float32 dask.array<shape=(288, 50, 13, 90, 90), chunksize=(1, 50, 13, 90, 90)>\n",
    "    \n",
    "The full shape of the array is (288, 50, 13, 90, 90), quite large. But the _chunksize_ is (1, 50, 13, 90, 90). Here the chunks correspond to the individual granuales of data (objects) in cloud storage. The chunk is the minimum amount of data we can read at one time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the size of the full dataset by typing the following, and dividing by $10^9$ to get the result in Gigabytes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.nbytes/(1e9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this tells us that the dataset is ~153 GB - that's pretty big! In fact, it's far larger than the memory of most personal computers, and so would be challenging to analyze without the help of something like cloud computing, combined with libraries like `xarray` and `dask`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a trick to make things work a bit faster\n",
    "coords = ds.coords.to_dataset().reset_coords()\n",
    "coords"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing Data\n",
    "\n",
    "### A Direct Plot\n",
    "\n",
    "Let's try to visualize something simple: the `Depth` variable. Here is how the data are stored:\n",
    "\n",
    "    Depth      (face, j, i) float32 dask.array<shape=(13, 90, 90), chunksize=(13, 90, 90)>\n",
    "\n",
    "Although depth is a 2D field, there is an extra, dimension (`face`) corresponding to the LLC face number. Let's use xarray's built in plotting functions to plot each face individually."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's check what shape `Depth` has."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords.Depth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that `Depth` is 3-dimensional with 13 faces. So if we simply use `Xarray`'s plotting function, we get the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords.Depth.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's a histogram! `Xarray` didn't know what kind of plot we wanted, so it binned values to make a histogram. But what we actually want is to see plots in terms of `i` and `j`, so we will set the columns argument as such: `col='face'`. The `col_wrap` call sets the number of columns to display in each row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords.Depth.plot(col='face', col_wrap=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This view is not the most useful. It reflects how the data is arranged logically, rather than geographically.\n",
    "\n",
    "### A Pretty Map\n",
    "\n",
    "To make plotting easier, we can define a quick function to plot the data in a more geographically friendly way. Eventually these plotting functions may be provided by the gcmplots package: <https://github.com/xecco/gcmplots>. For now, it is easy enough to roll our own.\n",
    "\n",
    "These are a bit complicated, so for now you can just take these as givens! Notice that we are using `cartopy` to make our plots here!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import cartopy as cart\n",
    "import pyresample\n",
    "\n",
    "class LLCMapper:\n",
    "\n",
    "    def __init__(self, ds, dx=0.25, dy=0.25):\n",
    "        # Extract LLC 2D coordinates\n",
    "        lons_1d = ds.XC.values.ravel()\n",
    "        lats_1d = ds.YC.values.ravel()\n",
    "\n",
    "        # Define original grid\n",
    "        self.orig_grid = pyresample.geometry.SwathDefinition(lons=lons_1d, lats=lats_1d)\n",
    "\n",
    "        # Longitudes latitudes to which we will we interpolate\n",
    "        lon_tmp = np.arange(-180, 180, dx) + dx/2\n",
    "        lat_tmp = np.arange(-90, 90, dy) + dy/2\n",
    "\n",
    "        # Define the lat lon points of the two parts.\n",
    "        self.new_grid_lon, self.new_grid_lat = np.meshgrid(lon_tmp, lat_tmp)\n",
    "        self.new_grid  = pyresample.geometry.GridDefinition(lons=self.new_grid_lon,\n",
    "                                                            lats=self.new_grid_lat)\n",
    "\n",
    "    def __call__(self, da, ax=None, projection=cart.crs.Robinson(), lon_0=-60, **plt_kwargs):\n",
    "\n",
    "        assert set(da.dims) == set(['face', 'j', 'i']), \"da must have dimensions ['face', 'j', 'i']\"\n",
    "\n",
    "        if ax is None:\n",
    "            fig, ax = plt.subplots(figsize=(12, 6), subplot_kw={'projection': projection})\n",
    "        else:\n",
    "            m = plt.axes(projection=projection)\n",
    "            \n",
    "        field = pyresample.kd_tree.resample_nearest(self.orig_grid, da.values,\n",
    "                                                    self.new_grid,\n",
    "                                                    radius_of_influence=100000,\n",
    "                                                    fill_value=None)\n",
    "\n",
    "        vmax = plt_kwargs.pop('vmax', field.max())\n",
    "        vmin = plt_kwargs.pop('vmin', field.min())\n",
    "\n",
    "        \n",
    "        x,y = self.new_grid_lon, self.new_grid_lat\n",
    "\n",
    "        # Find index where data is splitted for mapping\n",
    "        split_lon_idx = round(x.shape[1]/(360/(lon_0 if lon_0>0 else lon_0+360)))\n",
    "\n",
    "\n",
    "        p = ax.pcolormesh(x[:,:split_lon_idx], y[:,:split_lon_idx], field[:,:split_lon_idx],\n",
    "                         vmax=vmax, vmin=vmin, transform=cart.crs.PlateCarree(), zorder=1, **plt_kwargs)\n",
    "        p = ax.pcolormesh(x[:,split_lon_idx:], y[:,split_lon_idx:], field[:,split_lon_idx:],\n",
    "                         vmax=vmax, vmin=vmin, transform=cart.crs.PlateCarree(), zorder=2, **plt_kwargs)\n",
    "\n",
    "        ax.add_feature(cart.feature.LAND, facecolor='0.5', zorder=3)\n",
    "        label = ''\n",
    "        if da.name is not None:\n",
    "            label = da.name\n",
    "        if 'units' in da.attrs:\n",
    "            label += ' [%s]' % da.attrs['units']\n",
    "        cb = plt.colorbar(p, shrink=0.4, label=label)\n",
    "        return ax\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapper = LLCMapper(coords)\n",
    "mapper(coords.Depth);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this with any 2D cell-centered LLC variable.\n",
    "\n",
    "## Selecting data\n",
    "\n",
    "The entire ECCOv4e3 dataset is contained in a single `Xarray.Dataset` object. How do we view specific pieces of data? This is handled by Xarray's [indexing and selecting functions](http://xarray.pydata.org/en/latest/indexing.html). To get the SST from January 2000, we use the `.sel()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sst = ds.THETA.sel(time='2000-01-15', k=0) # select k=0 for surface\n",
    "sst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Still no data has actually been downloaded. That doesn't happen until we call `.load()` explicitly or try to make a plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapper(sst, cmap='RdBu_r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Do some Calculations\n",
    "\n",
    "Now let's start doing something besides just plotting the existing data. For example, let's calculate the time-mean SST."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_sst = ds.THETA.sel(k=0).mean(dim='time')\n",
    "mean_sst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, no data was loaded. Instead, `mean_sst` is a symbolic representation of the data that needs to be pulled and the computations that need to be executed to produce the desired result. In this case, the 288 original chunks all need to be read from cloud storage. Dask coordinates this automatically for us. But it does take some time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time mean_sst.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapper(mean_sst, cmap='RdBu_r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Speeding things up with a Dask Cluster\n",
    "\n",
    "How can we speed things up? In general, the main bottleneck for this type of data analysis is the speed with which we can read the data. With cloud storage, the access is highly parallelizeable.\n",
    "\n",
    "From a Pangeo environment, we can create a [Dask cluster](https://distributed.dask.org/en/latest/) to spread the work out amongst many compute nodes. This works on both HPC and cloud. In the cloud, the compute nodes are provisioned on the fly and can be shut down as soon as we are done with our analysis.\n",
    "\n",
    "The code below will create a cluster with five compute nodes. It can take a few minutes to provision our nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask_gateway import GatewayCluster\n",
    "from dask.distributed import Client\n",
    "cluster = GatewayCluster()\n",
    "cluster.scale(5)\n",
    "client = Client(cluster)\n",
    "cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Click the `Dashboard` link above to view the Dask Dashboard."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we re-run the mean calculation. Note how the dashboard helps us visualize what the cluster is doing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time ds.THETA.isel(k=0).mean(dim='time').load()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatially-Integrated Heat Content Anomaly\n",
    "\n",
    "Now let's do something harder. We will calculate the horizontally integrated heat content anomaly for the full 3D model domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the monthly climatology\n",
    "theta_clim = ds.THETA.groupby('time.month').mean(dim='time')\n",
    "# the anomaly\n",
    "theta_anom = ds.THETA.groupby('time.month') - theta_clim\n",
    "rho0 = 1029\n",
    "cp = 3994\n",
    "ohc = rho0 * cp * (theta_anom *\n",
    "                   coords.rA *\n",
    "                   coords.hFacC).sum(dim=['face', 'j', 'i'])\n",
    "ohc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# actually load the data\n",
    "ohc.load()\n",
    "# put the depth coordinate back for plotting purposes\n",
    "ohc.coords['Z'] = coords.Z\n",
    "ohc.swap_dims({'k': 'Z'}).transpose().plot(vmax=1e20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial Derivatives: Heat Budget\n",
    "\n",
    "As our final exercise, we will do something much more complicated. We will compute the time-mean convergence of vertically-integrated heat fluxes. This is hard for several reasons.\n",
    "\n",
    "The first reason it is hard is because it involves variables located at different grid points.\n",
    "Following MITgcm conventions, xmitgcm (which produced this dataset) labels the center point with the coordinates `j, i`, the u-velocity point as `j, i_g`, and the v-velocity point as `j_g, i`. \n",
    "The horizontal advective heat flux variables are\n",
    "\n",
    "    ADVx_TH    (time, k, face, j, i_g) float32 dask.array<shape=(288, 50, 13, 90, 90), chunksize=(1, 50, 13, 90, 90)>\n",
    "    ADVy_TH    (time, k, face, j_g, i) float32 dask.array<shape=(288, 50, 13, 90, 90), chunksize=(1, 50, 13, 90, 90)>\n",
    "    \n",
    "Xarray won't allow us to add or multiply variables that have different dimensions, and xarray by itself doesn't understand how to transform from one grid position to another.\n",
    "\n",
    "**That's why [xgcm](https://xgcm.readthedocs.io/en/latest/) was created.**\n",
    "\n",
    "Xgcm allows us to create a `Grid` object, which understands how to interpolate and take differences in a way that is compatible with finite volume models such at MITgcm. Xgcm also works with many other models, including ROMS, POP, MOM5/6, NEMO, etc.\n",
    "\n",
    "A second reason this is hard is because of the complex topology connecting the different MITgcm faces. Fortunately xgcm also [supports this](https://xgcm.readthedocs.io/en/latest/grid_topology.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xgcm\n",
    "\n",
    "# define the connectivity between faces\n",
    "face_connections = {'face':\n",
    "                    {0: {'X':  ((12, 'Y', False), (3, 'X', False)),\n",
    "                         'Y':  (None,             (1, 'Y', False))},\n",
    "                     1: {'X':  ((11, 'Y', False), (4, 'X', False)),\n",
    "                         'Y':  ((0, 'Y', False),  (2, 'Y', False))},\n",
    "                     2: {'X':  ((10, 'Y', False), (5, 'X', False)),\n",
    "                         'Y':  ((1, 'Y', False),  (6, 'X', False))},\n",
    "                     3: {'X':  ((0, 'X', False),  (9, 'Y', False)),\n",
    "                         'Y':  (None,             (4, 'Y', False))},\n",
    "                     4: {'X':  ((1, 'X', False),  (8, 'Y', False)),\n",
    "                         'Y':  ((3, 'Y', False),  (5, 'Y', False))},\n",
    "                     5: {'X':  ((2, 'X', False),  (7, 'Y', False)),\n",
    "                         'Y':  ((4, 'Y', False),  (6, 'Y', False))},\n",
    "                     6: {'X':  ((2, 'Y', False),  (7, 'X', False)),\n",
    "                         'Y':  ((5, 'Y', False),  (10, 'X', False))},\n",
    "                     7: {'X':  ((6, 'X', False),  (8, 'X', False)),\n",
    "                         'Y':  ((5, 'X', False),  (10, 'Y', False))},\n",
    "                     8: {'X':  ((7, 'X', False),  (9, 'X', False)),\n",
    "                         'Y':  ((4, 'X', False),  (11, 'Y', False))},\n",
    "                     9: {'X':  ((8, 'X', False),  None),\n",
    "                         'Y':  ((3, 'X', False),  (12, 'Y', False))},\n",
    "                     10: {'X': ((6, 'Y', False),  (11, 'X', False)),\n",
    "                          'Y': ((7, 'Y', False),  (2, 'X', False))},\n",
    "                     11: {'X': ((10, 'X', False), (12, 'X', False)),\n",
    "                          'Y': ((8, 'Y', False),  (1, 'X', False))},\n",
    "                     12: {'X': ((11, 'X', False), None),\n",
    "                          'Y': ((9, 'Y', False),  (0, 'X', False))}}}\n",
    "\n",
    "# create the grid object\n",
    "grid = xgcm.Grid(ds, periodic=False, face_connections=face_connections)\n",
    "grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the `grid` object we created to take the divergence of a 2D vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "advx_th_vint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "advy_th_vint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# vertical integral and time mean of horizontal diffusive heat flux\n",
    "advx_th_vint = ds.ADVx_TH.sum(dim='k').mean(dim='time')\n",
    "advy_th_vint = ds.ADVy_TH.sum(dim='k').mean(dim='time')\n",
    "\n",
    "# difference in the x and y directions\n",
    "diff_ADV_th = grid.diff_2d_vector({'X': advx_th_vint.load(), 'Y': advy_th_vint.load()}, boundary='fill')\n",
    "# convergence\n",
    "conv_ADV_th = -diff_ADV_th['X'] - diff_ADV_th['Y']\n",
    "conv_ADV_th"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xgcm.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vertical integral and time mean of horizontal diffusive heat flux\n",
    "difx_th_vint = ds.DFxE_TH.sum(dim='k').mean(dim='time')\n",
    "dify_th_vint = ds.DFyE_TH.sum(dim='k').mean(dim='time')\n",
    "\n",
    "# difference in the x and y directions\n",
    "diff_DIF_th = grid.diff_2d_vector({'X': difx_th_vint.load(), 'Y': dify_th_vint.load()}, boundary='fill')\n",
    "# convergence\n",
    "conv_DIF_th = -diff_DIF_th['X'] - diff_DIF_th['Y']\n",
    "conv_DIF_th"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert to Watts / m^2 and load\n",
    "mean_adv_conv = rho0 * cp * (conv_ADV_th/coords.rA).fillna(0.).load()\n",
    "mean_dif_conv = rho0 * cp * (conv_DIF_th/coords.rA).fillna(0.).load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = mapper(mean_adv_conv, cmap='RdBu_r', vmax=300, vmin=-300);\n",
    "ax.set_title(r'Convergence of Advective Flux (W/m$^2$)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = mapper(mean_dif_conv, cmap='RdBu_r', vmax=300, vmin=-300)\n",
    "ax.set_title(r'Convergence of Diffusive Flux (W/m$^2$)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = mapper(mean_dif_conv + mean_adv_conv, cmap='RdBu_r', vmax=300, vmin=-300)\n",
    "ax.set_title(r'Convergence of Net Horizontal Flux (W/m$^2$)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = mapper(ds.TFLUX.mean(dim='time').load(), cmap='RdBu_r', vmax=300, vmin=-300);\n",
    "ax.set_title(r'Surface Heat Flux (W/m$^2$)');"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}