Tracker_Strains.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "    <h2><center>Global Glacier Velocity point Tracker</center></h2>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fetch the path for essential python scripts\n",
    "\n",
    "import os\n",
    "import sys\n",
    "path =  os.getcwd()\n",
    "# Avoid Windows users to have issues with how paths are written\n",
    "path = path.replace('\\\\','/')\n",
    "\n",
    "# Import python scripts from notebooks folder\n",
    "sys.path.append(path + '/scripts')\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import all the necessary packages\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.path as mpath\n",
    "import markdown \n",
    "\n",
    "matplotlib.rcParams['figure.figsize'] = [9, 5]\n",
    "matplotlib.rcParams[\"figure.autolayout\"] = True\n",
    "\n",
    "from velocity_widget import ITSLIVE\n",
    "velocity_widget = ITSLIVE()\n",
    "plt.close()\n",
    "from ipywidgets import widgets, HTML, Output\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import time\n",
    "from pyproj import Transformer\n",
    "from shapely.geometry import Point\n",
    "from shapely.geometry.polygon import Polygon\n",
    "import geopandas as gpd\n",
    "import ipyleaflet as ipl\n",
    "import math\n",
    "import pandas as pd\n",
    "import glob\n",
    "import seaborn\n",
    "import warnings\n",
    "import rasterio as rio\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instructions\n",
    "\n",
    "<br>\n",
    "\n",
    "This notebook will allow you to track pixels at the surface of your AOI, depending on their displacement from the velocity dataset ITS-LIVE.\n",
    "\n",
    "<br>\n",
    "\n",
    "1. Choose the type of dataset you want to work with. Just make sure you already downloaded some data.\n",
    "\n",
    "2. Select a date range. Make sure this range is at least the one you selected when you downloaded your data. It can be smaller if you want.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the type of dataset you want to download: 'Yearly' or 'Subyearly'\n",
    "type_dataset = \"Yearly\"\n",
    "\n",
    "# Select your date range (careful the datacubes get heavy quickly)\n",
    "sdate = '1998-01-01'\n",
    "edate = '2018-01-01'\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Select approximately the area you want to work on (if you already have it downloaded)\n",
    "\n",
    "<br>\n",
    "\n",
    "This section should select automatically the data closest to your AOI.\n",
    "Run the following cell.\n",
    "Instructions: same as for data downloading for Yearly dataset:\n",
    "- Left click once to put a marker\n",
    "- Left click a second time to plot the AOI's boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "dates_range = widgets.SelectionRangeSlider(\n",
    "    options=[i for i in range(546)],\n",
    "    index=(1, 120),\n",
    "    continuous_update=False,\n",
    "    description='Interval (days): ',\n",
    "    orientation='horizontal',\n",
    "    layout={'width': '90%',\n",
    "            'display': 'flex'},\n",
    "    style={'description_width': 'initial'})\n",
    "\n",
    "variables =  widgets.Dropdown(\n",
    "    options=['v', 'v_error', 'vx', 'vy'],\n",
    "    description='Variable: ',\n",
    "    disabled=False,\n",
    "    value='v',\n",
    "    layout={'width': '20%',\n",
    "            'display': 'flex'},\n",
    "    style={'description_width': 'initial'})\n",
    "\n",
    "plot_type =  widgets.Dropdown(\n",
    "    options=['location', 'satellite'],\n",
    "    description='Plot By: ',\n",
    "    disabled=False,\n",
    "    value='location',\n",
    "    layout={'width': '20%',\n",
    "            'display': 'flex'},\n",
    "    style={'description_width': 'initial'})\n",
    "\n",
    "plot_button =  widgets.Button(\n",
    "    description='Plot',\n",
    "    button_style='primary',\n",
    "    icon='line-chart',\n",
    "    style={'description_width': 'initial'})\n",
    "\n",
    "get_points =  widgets.Button(\n",
    "    description='Get points',\n",
    "    button_style='primary',\n",
    "    icon='line-chart',\n",
    "    style={'description_width': 'initial'})\n",
    "\n",
    "clear_button =  widgets.Button(\n",
    "    description='Clear Points',\n",
    "    # button_style='warning',\n",
    "    icon=\"trash\",\n",
    "    style={'description_width': 'initial'})\n",
    "\n",
    "latitude = widgets.BoundedFloatText(\n",
    "    value=0.0,\n",
    "    min=-90.0,\n",
    "    max=90.0,\n",
    "    step=0.1,\n",
    "    description='Lat: ',\n",
    "    disabled=False,\n",
    "    style={'description_width': 'initial'},\n",
    "    layout={'width': '20%',\n",
    "            'display': 'flex'},\n",
    ")\n",
    "\n",
    "longitude = widgets.BoundedFloatText(\n",
    "    value=0.0,\n",
    "    min=-180.0,\n",
    "    max=180.0,\n",
    "    step=0.1,\n",
    "    description='Lon: ',\n",
    "    disabled=False,\n",
    "    style={'description_width': 'initial'},\n",
    "    layout={'width': '20%',\n",
    "            'display': 'flex'},\n",
    ")\n",
    "\n",
    "add_button =  widgets.Button(\n",
    "    description='Add Point',\n",
    "    # button_style='info',\n",
    "    icon=\"map-marker\",\n",
    "    style={'description_width': 'initial'})\n",
    "\n",
    "include_running_mean =  widgets.Checkbox(\n",
    "            value=False,\n",
    "            description=\"Include Running Mean\",\n",
    "            style={'description_width': 'initial'},\n",
    "            disabled=False,\n",
    "            indent=False,\n",
    "            tooltip=\"Plot running mean through each time series\",\n",
    "            layout=widgets.Layout(width=\"25%\"),\n",
    "        )\n",
    "\n",
    "export_button = widgets.Button(\n",
    "    description='Export Data',\n",
    "    # button_style='info',\n",
    "    icon=\"file-export\",\n",
    "    style={'description_width': 'initial'})\n",
    "\n",
    "data_link = widgets.HTML(\n",
    "    value=\"<br>\"\n",
    ")\n",
    "\n",
    "# If this congiguration changes we need to rerun the cell.\n",
    "config = { \n",
    "    \"plot\": \"v\", # or other ITS_LIVE variables: vx, vy ...\n",
    "    \"min_separation_days\": 1,\n",
    "    \"max_separation_days\": 90,\n",
    "    \"color_by\": \"location\", # valid values: satellite, points\n",
    "    \"verbose\": True, # print operations\n",
    "    \"runnig_mean\": True,\n",
    "    \"coords\": {\n",
    "        \"latitude\": latitude,\n",
    "        \"longitude\": longitude\n",
    "    },\n",
    "    \"data_link\": data_link\n",
    "}\n",
    "\n",
    "\n",
    "def downloader(whatever):\n",
    "    print('Downloading...')\n",
    "    global pathsave\n",
    "    ######### YEARLY DATASET DOWNLOAD ########\n",
    "    if type_dataset == 'Yearly':\n",
    "\n",
    "            # Create list of years for the date range chosen earlier\n",
    "            list_years = np.arange(int(sdate.split('-')[0]), int(edate.split('-')[0])+1)\n",
    "\n",
    "            # Create path to the files\n",
    "            pathsave = velocity_widget.path_yearly_datacubes\n",
    "            os.makedirs(pathsave, exist_ok = True)\n",
    "            \n",
    "\n",
    "\n",
    "            for Y in range(len(list_years)):\n",
    "\n",
    "                    # Generate URL for the nc file\n",
    "                    url = f'{velocity_widget.url_region[0]}{int(list_years[Y])}.nc#mode=bytes'\n",
    "\n",
    "                    # Load datacube according to prerequisites (time, space and variables)\n",
    "                    start = time.time()\n",
    "                    xrds = xr.open_dataset(url\n",
    "                                            ).sel(x=slice(velocity_widget.xmin_proj, velocity_widget.xmax_proj),\n",
    "                                                y=slice(velocity_widget.ymax_proj, velocity_widget.ymin_proj)).load()\n",
    "\n",
    "                    print(f\"downloaded {list_years[Y]} spatial slice {time.time()-start:8.1f} seconds\")\n",
    "                    xrds.to_netcdf(f\"{pathsave}{list_years[Y]}.nc\")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    ######## SUBYEARLY DATASET DOWNLOAD ########\n",
    "    else:\n",
    "\n",
    "            # Create path to the files\n",
    "            pathsave = velocity_widget.path_subyearly_datacubes\n",
    "            os.makedirs(pathsave, exist_ok = True)\n",
    "\n",
    "            # Get the cube address\n",
    "            cubes = velocity_widget.dct.addresses\n",
    "            cubes = [*set(cubes)]\n",
    "\n",
    "            # List of variables to drop for the download (we drop everything but the variables written below)\n",
    "            variables_drop = [ele for ele in list(\n",
    "                    xr.open_dataset(cubes[0], engine='zarr').variables\n",
    "                    ) if ele not in ['mid_date','x','y','acquisition_date_img1', 'acquisition_date_img2', 'date_center', 'date_dt', 'satellite_img1','satellite_img2', 'v','vx','vy','roi_valid_percentage']\n",
    "            ]\n",
    "\n",
    "            \n",
    "            for n in range(len(cubes)):\n",
    "\n",
    "                    # Get the cube's URL\n",
    "                    url = cubes[n]\n",
    "\n",
    "                    # Load indices of slices above the quality threshold\n",
    "                    valid = xr.open_dataset(cubes[n], engine='zarr').roi_valid_percentage.values\n",
    "\n",
    "                    # Grab the time values\n",
    "                    t = xr.open_dataset(cubes[n], engine='zarr').mid_date.values\n",
    "\n",
    "                    # Create a time mask, based on the validity of layers and the custom date-range\n",
    "                    t_mask = np.logical_and(valid>threshold, np.logical_and(t>np.datetime64(sdate), t<np.datetime64(edate)))\n",
    "                    \n",
    "                    # Load datacube according to prerequisites (time, space and variables)\n",
    "                    start = time.time()\n",
    "                    xrds = xr.open_dataset(url,\n",
    "                                            engine='zarr',\n",
    "                                            drop_variables=variables_drop\n",
    "                                            ).sel(mid_date=t_mask,\n",
    "                                                x=slice(velocity_widget.xmin_proj, velocity_widget.xmax_proj),\n",
    "                                                y=slice(velocity_widget.ymax_proj, velocity_widget.ymin_proj)).load()\n",
    "\n",
    "                    print(f\"downloaded {cubes[n].split('/')[-1].split('.')[0]} spatial slice {time.time()-start:8.1f} seconds\")\n",
    "                    xrds.to_netcdf(f\"{pathsave}{cubes[n].split('/')[-1].split('.')[0]}_{sdate}_{edate}.nc\")\n",
    "    print('Done ! You can hit \"plot\" now')\n",
    "\n",
    "def plotter(whatever):\n",
    "    list_files = glob.glob(f'{pathsave}*.nc')\n",
    "\n",
    "    fig, ax = plt.subplots(len(list_files), figsize=(10,10))\n",
    "\n",
    "    if len(list_files) == 1:\n",
    "        if type_dataset == \"Yearly\":\n",
    "            ax.pcolormesh(xr.open_dataset(list_files[0]).x.values,xr.open_dataset(list_files[0]).y.values,xr.open_dataset(list_files[0]).v.values)\n",
    "        else:\n",
    "            ax.pcolormesh(xr.open_dataset(list_files[0]).x.values,xr.open_dataset(list_files[0]).y.values,np.nanmean(xr.open_dataset(list_files[0]).v.values,axis = 0))\n",
    "    else:\n",
    "        if type_dataset == \"Yearly\":\n",
    "            for i in range(len(list_files)):\n",
    "                ax[i].pcolormesh(xr.open_dataset(list_files[i]).x.values,xr.open_dataset(list_files[i]).y.values,xr.open_dataset(list_files[i]).v.values)\n",
    "        else:\n",
    "            for i in range(len(list_files)):\n",
    "                ax[i].pcolormesh(xr.open_dataset(list_files[i]).x.values,xr.open_dataset(list_files[i]).y.values,np.nanmean(xr.open_dataset(list_files[i]).v.values,axis = 0))\n",
    "\n",
    "\n",
    "def update_variable(change):\n",
    "        if change['type'] == 'change' and change['name'] == 'value':\n",
    "            config[\"plot\"] = variables.value\n",
    "            velocity_widget.set_config(config)\n",
    "            velocity_widget.plot_time_series()\n",
    "            \n",
    "def update_range(change):\n",
    "        if change['type'] == 'change' and change['name'] == 'value':\n",
    "            start, end = change['new']\n",
    "            config[\"min_separation_days\"] = start\n",
    "            config[\"max_separation_days\"] = end\n",
    "            velocity_widget.set_config(config)\n",
    "            velocity_widget.plot_time_series()\n",
    "            \n",
    "def update_plottype(change):\n",
    "        if change['type'] == 'change' and change['name'] == 'value':\n",
    "            config[\"color_by\"] = plot_type.value\n",
    "            velocity_widget.set_config(config)\n",
    "            velocity_widget.plot_time_series()\n",
    "            \n",
    "def update_mean(change):\n",
    "        if change['type'] == 'change' and change['name'] == 'value':\n",
    "            config[\"running_mean\"] = include_running_mean.value\n",
    "            velocity_widget.set_config(config)\n",
    "            velocity_widget.plot_time_series()\n",
    "            \n",
    "def add_point(event):\n",
    "    coords = (latitude.value, longitude.value)\n",
    "    velocity_widget.add_point(coords)\n",
    "    \n",
    "def export_ts(event):\n",
    "    velocity_widget.export_data()\n",
    "\n",
    "get_points.on_click(velocity_widget.plot_time_series)\n",
    "plot_button.on_click(plotter)\n",
    "clear_button.on_click(velocity_widget.clear_points)\n",
    "\n",
    "export_button.on_click(downloader)\n",
    "\n",
    "\n",
    "add_button.on_click(add_point)\n",
    "dates_range.observe(update_range, 'value')\n",
    "plot_type.observe(update_plottype, 'value')\n",
    "variables.observe(update_variable, 'value')\n",
    "include_running_mean.observe(update_mean, 'value')\n",
    "\n",
    "layout = widgets.Layout(align_items='stretch',\n",
    "                        display='flex',\n",
    "                        flex_flow='row wrap',\n",
    "                        border='none',\n",
    "                        grid_template_columns=\"repeat(auto-fit, minmax(720px, 1fr))\",\n",
    "                        # grid_template_columns='48% 48%',\n",
    "                        width='99%',\n",
    "                        height='100%')\n",
    "\n",
    "velocity_widget.set_config(config)\n",
    "\n",
    "velocity_widget.fig.canvas.capture_scroll = True\n",
    "\n",
    "# We render the widget\n",
    "widgets.GridBox([\n",
    "                widgets.VBox([velocity_widget.map,\n",
    "                            widgets.HBox([latitude, longitude, add_button, clear_button], layout=widgets.Layout(align_items=\"flex-start\",\n",
    "                                                                                                                flex_flow='row wrap'))],\n",
    "                            layout=widgets.Layout(min_width=\"100%\",\n",
    "                                                    display=\"flex\",\n",
    "                                                    # height=\"100%\",\n",
    "                                                    # max_height=\"100%\",\n",
    "                                                    max_width=\"100%\"))],\n",
    "\n",
    "                layout=layout)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Function to get the file's names in the folder corresponding to your AOI**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def folder_finder(name_region):\n",
    "        # Gather list of folders in the region of the AOI\n",
    "        list_folders = glob.glob(f'{path}/Datacubes/{type_dataset}/{name_region[0]}/*')\n",
    "        # Make sure that Windows users won't encounter troubles with the '\\\\' character in the path\n",
    "        list_folders = [(list_folders[i].replace('\\\\','/')) for i in range(len(list_folders))]\n",
    "        # Calculate the center-point of the folders based on their respective AOI\n",
    "        center_folders = [((float(list_folders[i].split('/')[-1].split('_')[2])-float(list_folders[i].split('/')[-1].split('_')[0]))/2+float(list_folders[i].split('/')[-1].split('_')[0]),\n",
    "                        (float(list_folders[i].split('/')[-1].split('_')[3])-float(list_folders[i].split('/')[-1].split('_')[1]))/2+float(list_folders[i].split('/')[-1].split('_')[1])) for i in range(len(list_folders))]\n",
    "        # Calculate the distance between the center of your AOI and the folders' AOI's centers\n",
    "        dist = [np.sqrt((center_folders[i][0]-velocity_widget.point_center.coords[0][0])**2 + (\n",
    "                center_folders[i][1]-velocity_widget.point_center.coords[0][1])**2) for i in range(len(center_folders))]\n",
    "        # Get the index of the folder the closest to your AOI\n",
    "        list_files = glob.glob(f'{list_folders[np.where(dist==min(dist))[0][0]]}/*.nc')\n",
    "        list_files = [(list_files[i].replace('\\\\','/')) for i in range(len(list_files))]\n",
    "        # Update the path to save files\n",
    "        if type_dataset == 'Yearly':\n",
    "                pathsave = velocity_widget.path_yearly_datacubes.replace(velocity_widget.path_yearly_datacubes.split('/')[-2], list_files[0].split('/')[-2] ).replace('Datacubes','Output')\n",
    "\n",
    "        else:\n",
    "                pathsave = velocity_widget.path_subyearly_datacubes.replace(velocity_widget.path_subyearly_datacubes.split('/')[-2], list_files[0].split('/')[-2] ).replace('Datacubes','Output')\n",
    "\n",
    "        # Create output folder\n",
    "        os.makedirs(pathsave, exist_ok = True)\n",
    "\n",
    "        return list_files, pathsave\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the list of files in the folder corresponding to your AOI\n",
    "list_files, pathsave = folder_finder(velocity_widget.name_region)\n",
    "\n",
    "# Get the coordinates of your AOI's boundaries\n",
    "coords = [(float(list_files[0].split('/')[-2].split('_')[1]), (float(list_files[0].split('/')[-2].split('_')[0]))), (\n",
    "           float(list_files[0].split('/')[-2].split('_')[3]), (float(list_files[0].split('/')[-2].split('_')[2]))) ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define your analysis Frequency:\n",
    "\n",
    "<br>\n",
    "\n",
    "In the next cell, you will choose on which time interval to calculate variable's means.\n",
    "\n",
    "The variable \"TYPE\" allows you to choose several options:\n",
    "- \"Y\" for year\n",
    "- \"M\" for end-of-month\n",
    "- \"MS\" for start-of-month\n",
    "- 'D\" for day\n",
    "\n",
    "<br>\n",
    "\n",
    "The variable \"amount\" allows you to choose over how many of your \"TYPE\" variable you want to average.\n",
    "\n",
    "For example, a mean every 5 months starting at the beginning of each month would have the variables:\n",
    "- TYPE = \"MS\"\n",
    "- amount = 5\n",
    "\n",
    "<br>\n",
    "\n",
    "Finally, \"delta\" calculates the time difference between each step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the variable type\n",
    "TYPE = 'Y'\n",
    "\n",
    "# Select the amount of units you want to average for (eg: if TYPE = \"Y\" and amount = 1, average every 1 year)\n",
    "# I personnally like to have a 5 days variable for subyearly (TYPE = 'D', amount = 5)\n",
    "amount = 1\n",
    "\n",
    "# Get timestep\n",
    "delta = np.timedelta64(amount, TYPE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Functions to import the datasets**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def importer(list_files, datatype, TYPE, sdate, edate):\n",
    "    \n",
    "    if datatype == 'Yearly':   \n",
    "        \n",
    "        # Initialize an array that will host the dates \n",
    "        years = np.zeros((len(list_files)), dtype='datetime64[s]')\n",
    "        list_years = np.zeros((len(list_files)), dtype='int')\n",
    "        \n",
    "        # Get the date of each file if they are in the boundaries in the starting date and ending date, and update the list of files to process\n",
    "        for i in range(0,len(list_files)):\n",
    "            years = [pd.to_datetime(list_files[i].split('/')[-1].split('.')[0],format='%Y') for i in range(len(list_files))]\n",
    "            list_years = [list_files[i].split('/')[-1].split('.')[0] for i in range(len(list_files)) if np.logical_and(np.datetime64((list_files[i].split('/')[-1].split('.')[0]))>=np.datetime64(sdate), np.datetime64((list_files[i].split('/')[-1].split('.')[0]))<=np.datetime64(edate))]\n",
    "            list_files = [list_files[i] for i in range(len(list_files)) if np.logical_and(np.datetime64((list_files[i].split('/')[-1].split('.')[0]))>=np.datetime64(sdate), np.datetime64((list_files[i].split('/')[-1].split('.')[0]))<=np.datetime64(edate))]\n",
    "        \n",
    "        # Initialize a dataset on which we will append the others\n",
    "        # List of variables to drop for the download (we drop everything but the variables written below)\n",
    "        drop = [ele for ele in list(\n",
    "                xr.open_dataset(list_files[0]).variables\n",
    "                ) if ele not in ['x','y','vx','vy','v']\n",
    "        ]\n",
    "\n",
    "                    \n",
    "        ds = xr.open_dataset(list_files[0], \n",
    "                             drop_variables=drop\n",
    "                             )\n",
    "        \n",
    "        # Loop through every file, which we will append to the existing dataset\n",
    "        for i in range(1,len(list_files)):\n",
    "\n",
    "            print(f'Opening {list_files[i]} ')\n",
    "            \n",
    "            drop = [ele for ele in list(\n",
    "                    xr.open_dataset(list_files[i]).variables\n",
    "                    ) if ele not in ['x','y','vx','vy','v']\n",
    "            ]\n",
    "              \n",
    "            dtemp = xr.open_dataset(list_files[i],\n",
    "                                    drop_variables=drop)  \n",
    "                                            \n",
    "            ds = xr.concat((ds,dtemp), 'time')\n",
    "        \n",
    "        # Assign the dates to the time dimension\n",
    "        ds = ds.assign_coords({'time':years})\n",
    "        \n",
    "    else:\n",
    "        \n",
    "        list_years = []\n",
    "        years = []\n",
    "        \n",
    "        # Grab the time values\n",
    "        t = xr.open_dataset(list_files[0]).mid_date.values\n",
    "\n",
    "        # Create a time mask, based on the validity of layers and the custom date-range\n",
    "        t_mask = np.logical_and(t>=np.datetime64(sdate), t<=np.datetime64(edate))\n",
    "\n",
    "        ds = xr.open_dataset(list_files[0]).sel(mid_date=t_mask)\n",
    "        ds = ds.sortby(ds.mid_date).resample(mid_date=f'{amount}{TYPE}').mean(dim=\"mid_date\",skipna = True)\n",
    "\n",
    "        for i in range(len(list_files)):\n",
    "            t = xr.open_dataset(list_files[i]).mid_date.values\n",
    "            t_mask = np.logical_and(t>=np.datetime64(sdate), t<=np.datetime64(edate))\n",
    "            temp = xr.open_dataset(list_files[i]).sel(mid_date=t_mask)\n",
    "            temp = temp.sortby(temp.mid_date).resample(mid_date=f'{amount}{TYPE}').mean(dim=\"mid_date\", skipna = True)\n",
    "            ds = xr.concat((ds,temp),dim='mid_date')\n",
    "        \n",
    "        # Sort the time index to allow time slicing\n",
    "        ds = ds.sortby(ds.mid_date).resample(mid_date=f'{amount}{TYPE}').mean(dim=\"mid_date\", skipna = True)\n",
    "        \n",
    "        # Rename the time dimension so we don't have to modify the entire code\n",
    "        ds = ds.rename({'mid_date': 'time'})\n",
    "    \n",
    "    return ds, years, list_years"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the files as an xarray object, their dates and years\n",
    "ds, years, list_years = importer(list_files, type_dataset, TYPE, sdate, edate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Select the points you want to track\n",
    "\n",
    "<br>\n",
    "\n",
    "To do that:\n",
    "- Left click to select points to track\n",
    "- Right click to change the group of points\n",
    "\n",
    "By changing the group of points, you can track several transects at once.\n",
    "For example you can select several transects along flow to assess how the ice flows along the glacier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cell to plot the map and get coordinates on click and store them in an array**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def handle_click(**kwargs):\n",
    "    global i\n",
    "    global n\n",
    "    global c\n",
    "    # If left click, add marker, grab coordinates and store in the coordinates array\n",
    "    # Ideally, I would like 'click' to be 'leftclick' if that exists\n",
    "    if kwargs.get('type') == 'click':\n",
    "        \n",
    "        icon = ipl.AwesomeIcon(name='fa-cog', marker_color = colors[c])\n",
    "        m.add_layer(ipl.Marker(location=kwargs.get('coordinates'), icon=icon))\n",
    "        coordinates[n,i,:] = (kwargs['coordinates'])\n",
    "        i += 1\n",
    "\n",
    "    \n",
    "    # If right click, change color of marker, change column in coordinates, reinitialize the row counter\n",
    "    elif kwargs.get('type') == 'contextmenu':\n",
    "        i = 0\n",
    "        n += 1\n",
    "        c += 1\n",
    "        if c > len(colors)-1:\n",
    "            c = 0\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Initialize coordinates array\n",
    "coordinates=np.full((1000,1000,2), np.nan)\n",
    "\n",
    "# Initialize indices\n",
    "i = 0\n",
    "n = 0\n",
    "c = 0\n",
    "# Create a rotating array of colors to differenciate the markers\n",
    "colors = ['red', 'green', 'blue', 'purple', 'gray', 'orange', 'beige']\n",
    "\n",
    "# plot the map\n",
    "m = ipl.Map(\n",
    "    basemap=ipl.basemap_to_tiles(ipl.basemaps.Esri.WorldImagery),\n",
    "    center=((coords[1][0]-coords[0][0])/2 + coords[0][0],\n",
    "            (coords[1][1]-coords[0][1])/2 + coords[0][1]),\n",
    "    zoom=7,\n",
    "    scroll_wheel_zoom=True\n",
    "    )\n",
    "\n",
    "# plot the AOI\n",
    "rectangle = ipl.Rectangle(bounds=((coords[1][0], coords[1][1]), (coords[0][0], coords[0][1])))\n",
    "m.add_layer(rectangle)\n",
    "\n",
    "m.on_interaction(handle_click)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate the displacement of each point selected over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    \n",
    "def get_vels(pts, ind_nan):\n",
    "    \n",
    "    # Create time position array\n",
    "    pts_time = np.zeros((ds.time.shape[0], ind_nan, pts.shape[1]), dtype = int)\n",
    "    pts_time[0] = pts\n",
    "\n",
    "    \n",
    "    for i in range(1,ds.time.shape[0]):\n",
    "        \n",
    "        time_step = ds.time.values[i] - ds.time.values[i-1]\n",
    "        time_step = time_step.astype('timedelta64[D]')\n",
    "        time_step = time_step / np.timedelta64(365, 'D')\n",
    "        \n",
    "        \n",
    "        for j in range(0, ind_nan):\n",
    "            \n",
    "            if np.isnan(ds.vx.isel(time = i).sel(x=pts_time[i-1,j,0], y=pts_time[i-1,j,1], method=\"nearest\").values):\n",
    "                pts_time[i,j,0] = pts_time[i-1,j,0]\n",
    "            else:\n",
    "                pts_time[i,j,0] = pts_time[i-1,j,0] + ds.vx.isel(time = i).sel(x=pts_time[i-1,j,0], y=pts_time[i-1,j,1], method=\"nearest\").values*time_step\n",
    "                \n",
    "            if np.isnan(ds.vy.isel(time = i).sel(x=pts_time[i-1,j,0], y=pts_time[i-1,j,1], method=\"nearest\").values):\n",
    "                pts_time[i,j,1] = pts_time[i-1,j,1]\n",
    "            else:\n",
    "                pts_time[i,j,1] = pts_time[i-1,j,1] + ds.vy.isel(time = i).sel(x=pts_time[i-1,j,0], y=pts_time[i-1,j,1], method=\"nearest\").values*time_step\n",
    "          \n",
    "    return pts_time\n",
    "\n",
    "\n",
    "# Convert the coordinates to the dataset's projection\n",
    "transformer = Transformer.from_crs(4326, int(velocity_widget.proj_region[0]))\n",
    "coords_proj = np.array([[pt for pt in transformer.itransform(coordinates[i])] for i in range(coordinates.shape[0])])\n",
    "\n",
    "# Determine how many points there are per transect\n",
    "inds_nan = [np.where(np.isnan(coords_proj[i,:,0]))[0][0] for i in range(0,n+1)]\n",
    "\n",
    "# Fill the array with points taken from the map\n",
    "pts_total = np.full((coords_proj.shape[0], len(ds.time), coords_proj.shape[1], 2), np.nan)\n",
    "\n",
    "for i in range(0,n+1):\n",
    "    temp = get_vels(coords_proj[i,0:inds_nan[i]], inds_nan[i])\n",
    "    pts_total[i,:,0:inds_nan[i],:] = temp   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot the displacements\n",
    "- First position of points: orange\n",
    "- Last position of points: black\n",
    "- Intermediate positions: gradient from blue to white"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Used for the colors of the scatter plots\n",
    "time = np.arange(0,len(ds.time))\n",
    "# Background image\n",
    "prog = np.nanmean(ds[list(ds.keys())[0]].values, axis = 0)\n",
    "\n",
    "# plot the displacement of points\n",
    "plt.figure(num = 2, figsize=(10,10))\n",
    "plt.pcolormesh(ds.x, ds.y, prog, cmap = 'icefire')\n",
    "b = plt.colorbar()\n",
    "b.set_label('Velocity [m/yr]')\n",
    "\n",
    "# Plot each timestep\n",
    "for i in range(0,n+1):\n",
    "    for j in range(0, inds_nan[i]):\n",
    "        plt.scatter(pts_total[i,:,j,0],pts_total[i,:,j,1], c=time, cmap = 'Blues_r')\n",
    "    \n",
    "    # Plot the first and last positions: first is orange, last is black\n",
    "    plt.plot(pts_total[i,0,:,0], pts_total[i,0,:,1],'-o', color='orange')\n",
    "    plt.plot(pts_total[i,-1,:,0], pts_total[i,-1,:,1],'-o', color='black')\n",
    "\n",
    "os.makedirs(f'{pathsave}Figures/Tracker', exist_ok = True)\n",
    "plt.savefig(f'{pathsave}Figures/Tracker')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Strain Rates\n",
    "\n",
    "<br>\n",
    "\n",
    "This section will be ran as it is. It uses the same time-range as for the pixel tracker.\n",
    "However you have some options ! \n",
    "\n",
    "<br>\n",
    "\n",
    "**Instructions:**\n",
    "\n",
    "Run the next two cells. If you have an ice-thickness tif for you AOI, change \"flux_div\" to \"True\" (in the next cell) and indicate the path to:\n",
    "- the path to ice thickness tif file\n",
    "- the path to boundary file (supposed to be a 2 columns, csv file. See the template in the \"Example\" folder)\n",
    "\n",
    "<br>\n",
    "\n",
    "**WARNING:** Make sure the boundary and ice thickness are IN THE SAME PROJECTION as your downloaded datacubes.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_div = False\n",
    "\n",
    "path_thickness = f'{path}/Example/Thickness.tif'\n",
    "path_boundaries = f'{path}/Example/Boundaries.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if flux_div:\n",
    "\n",
    "        # Import thickness boundaries\n",
    "        boundary_points = np.genfromtxt(path_boundaries, delimiter=',')\n",
    "        boundary_points = boundary_points[1:]\n",
    "        boundary_points = boundary_points[:,2:]\n",
    "\n",
    "        # Create a mask for thickness\n",
    "        xx = ds.x\n",
    "        yy = ds.y\n",
    "        xxx = np.linspace(0, len(xx), len(xx)+1)\n",
    "        yyy = np.linspace(0, len(yy), len(yy)+1)\n",
    "        gridx, gridy = np.meshgrid(xx, yy)\n",
    "\n",
    "        X = np.vstack((gridx.ravel(),gridy.ravel())).T\n",
    "\n",
    "        p = mpath.Path(np.vstack((boundary_points[:,0], boundary_points[:,1])).T)\n",
    "        glacier_mask = p.contains_points(X)\n",
    "        glacier_mask = np.hstack([glacier_mask]).reshape(ds.vx.shape[1],ds.vx.shape[2])\n",
    "\n",
    "\n",
    "        # Import thickness tif which we will use to define our search window for the xarray\n",
    "        with rio.open(path_thickness) as src:\n",
    "                thickness = src.read()[0]\n",
    "                height_thickness = thickness.shape[0]\n",
    "                width_thickness = thickness.shape[1]\n",
    "                cols_thickness, rows_thickness = np.meshgrid(np.arange(width_thickness), np.arange(height_thickness))\n",
    "                xs_thickness, ys_thickness = rio.transform.xy(src.transform, rows_thickness, cols_thickness)\n",
    "                lons_thickness= np.array(xs_thickness)[0]\n",
    "                lats_thickness = np.array(ys_thickness)[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grabs the closest indices of an array based on the input\n",
    "def closest_value(pty, ptx, array):\n",
    "    \n",
    "    minx = np.abs(lons_thickness-ptx)\n",
    "    miny = np.abs(lats_thickness-pty)\n",
    "    \n",
    "    col = np.where(minx == minx.min())[0][0]\n",
    "    row = np.where(miny == miny.min())[0][0]\n",
    "    \n",
    "    return row, col\n",
    "\n",
    "\n",
    "\n",
    "# Calculate strains and dynamic thinning\n",
    "def Strains_Calculator(datacube, sdate, edate, upscale, freq, Plotting, flux_div):\n",
    "\n",
    "    xx = ds.x\n",
    "    yy = ds.y\n",
    "    \n",
    "    # We suppress the warning for means of empty slice otherwise it pops-up often\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.filterwarnings(action='ignore', message='Mean of empty slice')\n",
    "        vxtot = np.nanmean(datacube.vx.sel(time = slice(sdate, edate)), axis = 0)\n",
    "        vytot = np.nanmean(datacube.vy.sel(time = slice(sdate, edate)), axis = 0)\n",
    "\n",
    "    # Calculate spatial steps\n",
    "    dx = np.abs(xx[0] - xx[upscale])\n",
    "    dy = np.abs(yy[0] - yy[upscale]) \n",
    "    \n",
    "    # Calculate spatial steps\n",
    "    dx = np.abs(xx[0] - xx[upscale])\n",
    "    dy = np.abs(yy[0] - yy[upscale])\n",
    "    \n",
    "    # Rescale X and Y\n",
    "    YY = yy.values[::upscale]\n",
    "    YY = YY[1:-1]\n",
    "    XX = xx.values[::upscale]\n",
    "    XX = XX[1:-1]\n",
    "    gridX, gridY = np.meshgrid(XX,YY)\n",
    "\n",
    "    #### Strain Rates #### \n",
    "    \n",
    "    # Prepare arrays for the strain calculation\n",
    "    vx = vxtot[::upscale, ::upscale]\n",
    "    vy = vytot[::upscale, ::upscale]\n",
    "      \n",
    "    # Calculate the spatial step\n",
    "    dx = np.abs(xx[0] - xx[upscale])\n",
    "    dy = np.abs(yy[0] - yy[upscale])\n",
    "    xxx = np.linspace(0, len(xx), len(xx)+1)\n",
    "    yyy = np.linspace(0, len(yy), len(yy)+1)   \n",
    "    \n",
    "    # Rescale rows and cols\n",
    "    dyy = yyy[::upscale]\n",
    "    dyy = dyy[1:-1]\n",
    "    dxx = xxx[::upscale]\n",
    "    dxx = dxx[1:-1]\n",
    "    \n",
    "    # Calculate velocity derivatives components\n",
    "    duy = (vx[2:,1:-1] - vx[:-2,1:-1])/(2*dx.values)\n",
    "    dux = (vx[1:-1,2:] - vx[1:-1,:-2])/(2*dy.values)\n",
    "    dvy = (vy[2:,1:-1] - vy[:-2,1:-1])/(2*dx.values)\n",
    "    dvx = (vy[1:-1,2:] - vy[1:-1,:-2])/(2*dy.values) \n",
    "\n",
    "    # Calculate the strain components\n",
    "    exx = dux\n",
    "    eyy = dvy\n",
    "    # Based on the assumption of flux divergence = 0\n",
    "    ezz = -exx-eyy\n",
    "    exy = 0.5*(duy + dvx)\n",
    "    \n",
    "    # Calculate the eigenvalues of the strain rate tensor\n",
    "    e1 = 0.5*(exx+eyy) + 0.5*np.sqrt(4*exy**2 + (exx-eyy)**2)\n",
    "    e2 = 0.5*(exx+eyy) - 0.5*np.sqrt(4*exy**2 + (exx-eyy)**2)\n",
    "\n",
    "    if flux_div:\n",
    "    #### Flux Divergence ####\n",
    "        \n",
    "        # Scale the thickness to the dataset\n",
    "        thickness_scaled = np.zeros(vxtot.shape)\n",
    "        \n",
    "        # Filter out non-glacier cells\n",
    "        thickness_scaled[glacier_mask==False] = np.nan\n",
    "        \n",
    "        # Resize the thickness to the dataset\n",
    "        for j in range(vxtot.shape[0]):\n",
    "            for i in range(vxtot.shape[1]):\n",
    "                if ~np.isnan(thickness_scaled[j,i]):    \n",
    "                    row, col = closest_value(yy.values[j], xx.values[i], thickness)\n",
    "                    thickness_scaled[j,i] = thickness[row,col]   \n",
    "        \n",
    "        # Filter out erratic values\n",
    "        thickness_scaled[thickness_scaled==np.nanmin(thickness_scaled)]=np.nan\n",
    "        \n",
    "        # Calculate thickness components qx and qy\n",
    "        qx = vxtot[::upscale,::upscale]*thickness_scaled[::upscale,::upscale]\n",
    "        qy = vytot[::upscale,::upscale]*thickness_scaled[::upscale,::upscale]\n",
    "    \n",
    "        # Calculate dqx and dqy\n",
    "        dqx = (qx[2:,1:-1] - qx[:-2,1:-1])/(2*dx.values)\n",
    "        dqy = (qy[1:-1,2:] - qy[1:-1,:-2])/(2*dy.values)\n",
    "        \n",
    "        # Calculate flux divergence\n",
    "        divflux = dqx + dqy\n",
    "    \n",
    "    else:\n",
    "        # Ensure divflux at least exists\n",
    "        divflux = np.zeros(ezz.shape)\n",
    "        divflux[divflux==0] = np.nan\n",
    "    \n",
    "    \n",
    "    if Plotting:       \n",
    "\n",
    "        # Make a colorbar for the plots\n",
    "        co = seaborn.diverging_palette(250, 30, l=65, center=\"dark\", as_cmap=True)\n",
    "\n",
    "        # plot flux divergence\n",
    "        fig = plt.figure(num=1, clear=True, figsize=[10,10])\n",
    "        ax = fig.add_subplot()\n",
    "        f = ax.pcolormesh(gridX, gridY, divflux, cmap = co, vmin = -100, vmax = 100)\n",
    "        fig.colorbar(f, ax = ax)\n",
    "        ax.set_title('Flux Divergence [m] '+sdate + ' / ' +edate, fontsize=18)\n",
    "        fig.savefig(f'{pathsave}Figures/Flux_div/{amount}_{TYPE}/{upscale}_{sdate}_{edate}.png', dpi = 400)\n",
    "                    \n",
    " \n",
    "        # plot 1st eigenvalue\n",
    "        fig = plt.figure(num=1, clear=True, figsize=[10,10])\n",
    "        ax = fig.add_subplot()\n",
    "        f = ax.pcolormesh(gridX, gridY, e1, cmap = co, vmin = -1, vmax = 1)\n",
    "        fig.colorbar(f, ax = ax)\n",
    "        ax.set_title('e1 '+sdate + ' / ' +edate, fontsize=18)\n",
    "        fig.savefig(f'{pathsave}Figures/e1/{amount}_{TYPE}/{upscale}_{sdate}_{edate}.png', dpi = 400)\n",
    "\n",
    "        # plot 2nd eigenvalue\n",
    "        fig = plt.figure(num=1, clear=True, figsize=[10,10])\n",
    "        ax = fig.add_subplot()\n",
    "        f = ax.pcolormesh(gridX, gridY, e2, cmap = co, vmin = -1, vmax = 1)\n",
    "        fig.colorbar(f, ax = ax)\n",
    "        ax.set_title('e2 '+sdate + ' / ' +edate, fontsize=18)\n",
    "        fig.savefig(f'{pathsave}Figures/e2/{amount}_{TYPE}/{upscale}_{sdate}_{edate}.png', dpi = 400)\n",
    "\n",
    "        # plot ezz eigenvalue\n",
    "        fig = plt.figure(num=1, clear=True, figsize=[10,10])\n",
    "        ax = fig.add_subplot()\n",
    "        f = ax.pcolormesh(gridX, gridY, ezz, cmap = 'bwr', vmin = -0.1, vmax = 0.1)\n",
    "        fig.colorbar(f, ax = ax)\n",
    "        ax.set_title('\\u03B5zz '+sdate + ' / ' +edate, fontsize=18)\n",
    "        fig.savefig(f'{pathsave}Figures/ezz/{amount}_{TYPE}/{upscale}_{sdate}_{edate}.png', dpi = 400)\n",
    "\n",
    "\n",
    "        \n",
    "    return ezz, divflux, e1, e2\n",
    "\n",
    "\n",
    "def Runner(Plotting, upscale, sdate, edate, flux_div):\n",
    "\n",
    "    # Create host subfolders\n",
    "    list_subfolders = ['Strains', 'Figures/e1', 'Figures/e2', 'Figures/ezz', 'Figures/Flux_div']\n",
    "    [os.makedirs(f'{pathsave}{list_subfolders[i]}/{amount}_{TYPE}/', exist_ok = True) for i in range(len(list_subfolders))]\n",
    "    \n",
    "    # Create a template array for the strain rate datacube.\n",
    "    template = np.zeros(ds.vx.shape)\n",
    "    template = template[:,::upscale, ::upscale]\n",
    "    template = template[:,1:-1,1:-1]\n",
    "    template[template==0] = np.nan\n",
    "    \n",
    "    # Rearrange the coordinates so they fit how we calculate the strain rates\n",
    "    xx = ds.x\n",
    "    yy = ds.y\n",
    "\n",
    "    Y = yy.values[::upscale]\n",
    "    Y = Y[1:-1]\n",
    "      \n",
    "    X = xx.values[::upscale]\n",
    "    X = X[1:-1]\n",
    "\n",
    "    # Create empty datacube to host the generated strain rates\n",
    "    strain_components = xr.Dataset(\n",
    "        data_vars=dict(\n",
    "            ezz=([\"time\", \"y\", \"x\"], template.copy()),\n",
    "            divflux=([\"time\", \"y\", \"x\"], template.copy()),\n",
    "            e1=([\"time\", \"y\", \"x\"], template.copy()),\n",
    "            e2=([\"time\", \"y\", \"x\"], template.copy()),\n",
    "            \n",
    "        ),\n",
    "        coords=dict(\n",
    "            time=([\"time\"],ds.time.values),\n",
    "            y=([\"y\"], Y),\n",
    "            x=([\"x\"], X),\n",
    "        ),\n",
    "        attrs=dict(description=f\"Strain Rate Datacube with an upscale = {upscale}\"),\n",
    "    )\n",
    "    \n",
    "    \n",
    "    # Loop through the strain calculation\n",
    "    for i in range(len(ds.time.values)-amount):\n",
    "\n",
    "        # At every loop, increment the spatial step\n",
    "        start = str(ds.time.values[i].astype(f'M8[{TYPE}]'))\n",
    "        end = str(ds.time.values[i+amount].astype(f'M8[{TYPE}]'))\n",
    "\n",
    "        # Calculate the strains\n",
    "        ezz, divflux, e1, e2 = Strains_Calculator(ds, start, end, upscale, TYPE, Plotting, flux_div)\n",
    "\n",
    "        # Save the calculated strains in the datacube\n",
    "        strain_components.ezz[i] = ezz\n",
    "        strain_components.divflux[i] = divflux\n",
    "        strain_components.e1[i] = e1\n",
    "        strain_components.e2[i] = e2\n",
    "        print(start)\n",
    "    \n",
    "    # Export the datacube\n",
    "    strain_components.to_netcdf(f'{pathsave}{list_subfolders[0]}/{amount}_{TYPE}/{sdate}_{edate}_{upscale}.nc')\n",
    "\n",
    "    return strain_components\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instructions\n",
    "\n",
    "<br>\n",
    "\n",
    "The following function can:\n",
    "- Calculate strain rates\n",
    "- Calculate flux divergence\n",
    "- Plot the strain rates\n",
    "- Plot arrows of the strain rates (beware, this option can take a lot of time depending on your upscale)\n",
    "\n",
    "<br>\n",
    "\n",
    "You can choose your input as it is written in the following cell. Just change statements to \"True\" or \"False\". It will print the date of each step processed.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "datacube = Runner(Plotting = True, upscale = 5, sdate = sdate, edate = edate, flux_div = flux_div)"
   ]
  }
 ],
 "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.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "2b7814ffd776579197ebe28bd31711c0ecb3137f77c21042d9f4200b395566ae"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}