s5p-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sentinel-5P Data Exploration\n",
    "\n",
    "In this notebook, we use [s5p-tools](https://github.com/bilelomrani1/s5p-tools) to automatically download and preprocess Sentinel-5P data. We provide starting code for preliminary analysis. The following notebook is shared under a MIT license. Fell free to reuse and adapt the code to your needs. If you find this code useful, please consider citing the following paper:\n",
    "```\n",
    "@article{OMRANI2020105089,\n",
    "title = \"Spatio-temporal data on the air pollutant nitrogen dioxide derived from Sentinel satellite for France\",\n",
    "journal = \"Data in Brief\",\n",
    "volume = \"28\",\n",
    "pages = \"105089\",\n",
    "year = \"2020\",\n",
    "issn = \"2352-3409\",\n",
    "doi = \"https://doi.org/10.1016/j.dib.2019.105089\",\n",
    "url = \"http://www.sciencedirect.com/science/article/pii/S2352340919314453\",\n",
    "author = \"Hichem Omrani and Bilel Omrani and Benoit Parmentier and Marco Helbich\",\n",
    "keywords = \"Air pollution, Remote sensing, Monitoring\"\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from itertools import product, cycle\n",
    "\n",
    "# Ignore Runtime warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore', category=RuntimeWarning)\n",
    "\n",
    "# Propagate attributes during computation\n",
    "xr.set_options(keep_attrs=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DS = xr.open_dataset('processed/processed__NO2___/NO2___1-6-2020__7-6-2020.nc')\n",
    "DS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preliminaries\n",
    "\n",
    "By default, each entry along the `time` dimension correspons to a single product. We resample the dataset so that we average together products belonging to the same day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DS = DS.resample(time='1D').mean(dim='time', skipna=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you wish to select a particular time slice, use the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DS = DS.sel(time=slice('2020-06-01', '2020-06-07'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable of interest is `tropospheric_NO2_column_number_density`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VARIABLE = 'tropospheric_NO2_column_number_density'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the column meta-data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DS[VARIABLE].attrs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We slightly reformat the labels and create a meta-data `long_name` which will be called during plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DS[VARIABLE].attrs['units'] = 'mol/m$^{2}$'\n",
    "DS[VARIABLE].attrs['long_name'] = 'Tropospheric vertical column of NO$_2$'\n",
    "\n",
    "LABEL = '{long_name}\\n [{unit}]'.format(long_name=DS[VARIABLE].attrs['long_name'],\n",
    "                                        unit=DS[VARIABLE].attrs['units'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a folder for the plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export directory\n",
    "from os import makedirs\n",
    "folder_plots = 'plots'\n",
    "makedirs(folder_plots, exist_ok=True)\n",
    "\n",
    "# Export format\n",
    "FORMAT_EXPORT = 'tiff'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some settings for plots..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "import matplotlib.dates as mdates\n",
    "import matplotlib.ticker as mticker\n",
    "import matplotlib.font_manager as fm\n",
    "from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar\n",
    "import seaborn as sns;sns.set();sns.set_style(\"white\")\n",
    "import cartopy.crs as ccrs\n",
    "%matplotlib inline\n",
    "\n",
    "# Properties\n",
    "COLOR_PALETTE = plt.cm.OrRd\n",
    "PROJECTION = ccrs.Mercator()\n",
    "\n",
    "# Scalebar\n",
    "SCALEBAR_LENGTH = 100000   # in meters\n",
    "SCALEBAR_UNIT = 'km'\n",
    "NUMBER_METERS_PER_UNIT = 1000\n",
    "SCALEBAR_LOCATION = 'lower left'\n",
    "SCALEBAR_SEP = 3\n",
    "SCALEBAR_TEXT_SIZE = 7\n",
    "\n",
    "# Font sizes\n",
    "SMALL_SIZE = 6\n",
    "MEDIUM_SIZE = 9\n",
    "BIGGER_SIZE = 18\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=MEDIUM_SIZE)    # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=MEDIUM_SIZE)   # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=MEDIUM_SIZE)   # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=MEDIUM_SIZE)   # legend fontsize\n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Global map averaged over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "FIGSIZE = 4\n",
    "\n",
    "# Temporal average\n",
    "data = DS.mean(dim='time')\n",
    "\n",
    "# Figure creation\n",
    "fig = plt.figure(figsize=(FIGSIZE*data.dims['longitude'] / data.dims['latitude'], FIGSIZE))\n",
    "ax = plt.axes(projection=PROJECTION)\n",
    "ax.set_facecolor('white')\n",
    "\n",
    "# Plot data\n",
    "p = data[VARIABLE].plot(\n",
    "        ax=ax,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        cmap=COLOR_PALETTE,\n",
    "        #norm=colors.LogNorm()\n",
    "    )\n",
    "\n",
    "# Add scale bar\n",
    "scalebar = AnchoredSizeBar(ax.transData,\n",
    "                           SCALEBAR_LENGTH, \n",
    "                           '{length} {unit}'.format(length=int(SCALEBAR_LENGTH/NUMBER_METERS_PER_UNIT), \n",
    "                                                    unit=SCALEBAR_UNIT),\n",
    "                           SCALEBAR_LOCATION, \n",
    "                           pad=0.1,\n",
    "                           fontproperties=fm.FontProperties(size=SCALEBAR_TEXT_SIZE),\n",
    "                           sep=SCALEBAR_SEP,\n",
    "                           color='black',\n",
    "                           label_top=True,\n",
    "                           frameon=False)\n",
    "ax.add_artist(scalebar)\n",
    "\n",
    "# Re-fit axis extent \n",
    "ax.set_extent([data.longitude.min(), data.longitude.max(), data.latitude.min(), data.latitude.max()],\n",
    "              ccrs.Geodetic())\n",
    "\n",
    "# Color bar\n",
    "p.colorbar.locator = mticker.AutoLocator()\n",
    "p.colorbar.formatter = mticker.ScalarFormatter(useMathText=True, useOffset=True)\n",
    "p.colorbar.formatter.set_scientific(True)\n",
    "p.colorbar.formatter.set_powerlimits((0,0))\n",
    "p.colorbar.ax.yaxis.set_offset_position('left')\n",
    "p.colorbar.ax.yaxis.labelpad = 20\n",
    "p.colorbar.ax.yaxis.set_tick_params(width=1)\n",
    "p.colorbar.update_ticks()\n",
    "\n",
    "# Export\n",
    "plt.savefig(\"{folder_plots}/france_global_NO2.{format}\".format(\n",
    "                folder_plots=folder_plots,\n",
    "                format=FORMAT_EXPORT),\n",
    "            bbox_inches='tight', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Administrative areas maps\n",
    "\n",
    "We generate regional maps with `fr-regions.shp` shapefile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.io.shapereader as shpreader\n",
    "import shapely.vectorized\n",
    "\n",
    "reader = shpreader.Reader('shp/fr-regions.shp')\n",
    "records = list(reader.records())\n",
    "\n",
    "pd.DataFrame([entry.attributes for entry in records])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To keep the legend consistent, we need to keep in memory the min and max value accross all the territory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = DS.mean(dim='time')[VARIABLE].values\n",
    "val_min = np.nanmin(values)\n",
    "val_max = np.nanmax(values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then iterate over the shapefile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "FIGSIZE = 1.85\n",
    "NROWS = 5\n",
    "NCOLS = 3\n",
    "\n",
    "SCALEBAR_LENGTH = 50000\n",
    "NAME_ATTR = 'nom'\n",
    "\n",
    "points = list(product(DS.latitude.values, DS.longitude.values))\n",
    "x = [i for j, i in points]\n",
    "y = [j for j, i in points]\n",
    "\n",
    "fig, axes = plt.subplots(\n",
    "    figsize=FIGSIZE * plt.figaspect(NROWS/NCOLS),\n",
    "    nrows=NROWS, ncols=NCOLS, subplot_kw={'projection': PROJECTION}\n",
    ")\n",
    "\n",
    "plt.subplots_adjust(hspace=0.5, wspace=0.2)\n",
    "\n",
    "for ax in axes.flat:\n",
    "    ax.set_facecolor('white')\n",
    "    ax.set(aspect=1, adjustable='datalim')\n",
    "\n",
    "index_ax = 0\n",
    "for index, region in enumerate(records):\n",
    "    ax = axes.flat[index_ax]\n",
    "    mask = shapely.vectorized.contains(region.geometry, x, y).reshape(\n",
    "        (DS.latitude.shape[0], DS.longitude.shape[0]))\n",
    "    if np.any(mask):\n",
    "        print(\"Processing {region}\".format(region=region.attributes[NAME_ATTR]))\n",
    "        mask_array = xr.DataArray(mask, dims=('latitude', 'longitude'),\n",
    "                                      coords={'longitude': DS.longitude, 'latitude': DS.latitude})\n",
    "        masked_data = DS.mean(dim='time').where(mask_array, drop=True)\n",
    "        im = ax.pcolormesh(masked_data.longitude, masked_data.latitude, masked_data[VARIABLE],\n",
    "                           transform=ccrs.PlateCarree(),\n",
    "                           #norm=colors.LogNorm(),\n",
    "                           vmin=val_min, vmax=val_max,\n",
    "                           cmap=COLOR_PALETTE)\n",
    "        \n",
    "        # Add scale bar\n",
    "        scalebar = AnchoredSizeBar(ax.transData,\n",
    "                       SCALEBAR_LENGTH, \n",
    "                       '{length} {unit}'.format(length=int(SCALEBAR_LENGTH/NUMBER_METERS_PER_UNIT), \n",
    "                                                unit=SCALEBAR_UNIT),\n",
    "                       SCALEBAR_LOCATION, \n",
    "                       pad=-0.6,\n",
    "                       sep=SCALEBAR_SEP,\n",
    "                       color='black',\n",
    "                       label_top=True,\n",
    "                       borderpad=-0.6,\n",
    "                       fontproperties=fm.FontProperties(size=SCALEBAR_TEXT_SIZE),\n",
    "                       frameon=False)\n",
    "        ax.add_artist(scalebar)\n",
    "\n",
    "        ax.title.set_text(region.attributes[NAME_ATTR])\n",
    "        ax.set_extent([masked_data.longitude.min(), masked_data.longitude.max(), \n",
    "                      masked_data.latitude.min(), masked_data.latitude.max()], ccrs.Geodetic())\n",
    "        index_ax = index_ax + 1\n",
    "    else:\n",
    "        print(\"{region} does not intersect AOI\".format(region=region.attributes[NAME_ATTR]))\n",
    "\n",
    "# Color bar\n",
    "colorbar = fig.colorbar(im, ax=axes.ravel().tolist(), aspect=40)\n",
    "colorbar.locator = mticker.AutoLocator()\n",
    "colorbar.formatter = mticker.ScalarFormatter(useMathText=True, useOffset=True)\n",
    "colorbar.formatter.set_scientific(True)\n",
    "colorbar.formatter.set_powerlimits((0,0))\n",
    "colorbar.ax.yaxis.set_offset_position('left')\n",
    "colorbar.ax.yaxis.labelpad = 20\n",
    "colorbar.ax.yaxis.set_tick_params(width=1)\n",
    "colorbar.set_label(LABEL)\n",
    "colorbar.update_ticks()\n",
    "\n",
    "plt.savefig(\"{folder_plots}/regions_global_NO2.{format}\".format(folder_plots=folder_plots,\n",
    "                                                                format=FORMAT_EXPORT),\n",
    "           bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monthly Boxplot\n",
    "\n",
    "We downsample the datset per month and look at the temporal variation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boxes_months = [group[VARIABLE].values for _,group in DS.mean(dim='longitude').mean(dim='latitude').resample({'time': '1MS'})]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6,2))\n",
    "ax = plt.axes()\n",
    "\n",
    "ax.yaxis.set_major_formatter(mticker.ScalarFormatter(useMathText=True, useOffset=True))\n",
    "\n",
    "boxplot = sns.boxplot(data=boxes_months,\n",
    "                      linewidth=0.8, width=0.5, saturation=0.5, showfliers=False, color=\"slategrey\")\n",
    "boxplot.set_xticklabels(DS.mean(dim='longitude').mean(dim='latitude').resample({'time': '1MS'}).mean('time').time.dt.strftime('%b %Y').values, rotation=30)\n",
    "plt.ylabel(LABEL, labelpad=10)\n",
    "\n",
    "plt.savefig(\"{folder_plots}/box_plot_NO2.{format}\".format(folder_plots=folder_plots,\n",
    "                                                        format=FORMAT_EXPORT), \n",
    "            bbox_inches='tight', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weekday boxplot\n",
    "\n",
    "Grouping the entries by weekday. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monday=0, Sunday=6\n",
    "boxes_weekday = [group[VARIABLE].values for _,group in DS.mean(dim='longitude').mean(dim='latitude').groupby('time.weekday')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "WEEKDAYS_LABEL = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']\n",
    "\n",
    "fig = plt.figure(figsize=(6,2))\n",
    "ax = plt.axes()\n",
    "ax.yaxis.set_major_formatter(mticker.ScalarFormatter(useMathText=True, useOffset=True))\n",
    "\n",
    "boxplot = sns.boxplot(data=boxes_weekday,\n",
    "                      linewidth=0.5, width=0.5, saturation=0.5, showfliers=False, color=\"slategrey\")\n",
    "boxplot.set_xticklabels(WEEKDAYS_LABEL)\n",
    "plt.ylabel(LABEL, labelpad=10)\n",
    "\n",
    "plt.savefig(\"{folder_plots}/box_plot_weekdays_NO2.{format}\".format(folder_plots=folder_plots,\n",
    "                                                                 format=FORMAT_EXPORT),\n",
    "            bbox_inches='tight', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Montly maps\n",
    "\n",
    "We resample the dataset per month;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DS_month = DS.resample({'time': '1MS'}).mean(dim='time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We keep in memory the min and max value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nat_avg_by_month_max = DS_month[VARIABLE].max(dim='longitude').max(dim='latitude')\n",
    "nat_avg_by_month_min = DS_month[VARIABLE].min(dim='longitude').min(dim='latitude')\n",
    "\n",
    "val_min = np.nanmin(nat_avg_by_month_min)\n",
    "val_max = np.nanmax(nat_avg_by_month_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FIGSIZE = 1.8\n",
    "NROWS = 5\n",
    "NCOLS = 3\n",
    "\n",
    "fig, axes = plt.subplots(\n",
    "    figsize=FIGSIZE * plt.figaspect(NROWS/NCOLS),\n",
    "    nrows=NROWS, ncols=NCOLS, subplot_kw={'projection': PROJECTION}\n",
    ")\n",
    "\n",
    "plt.subplots_adjust(wspace=0.2, hspace=0.3)\n",
    "\n",
    "for ax in axes.flat:\n",
    "    ax.set_facecolor('white')\n",
    "    ax.set(aspect=1, adjustable='datalim')\n",
    "\n",
    "index_ax = 0\n",
    "for index, time in enumerate(DS_month.time):\n",
    "    print('Processing {month}'.format(month=time.dt.strftime('%b %Y').values))\n",
    "    ax = axes.flat[index_ax]\n",
    "    data = DS_month.isel(time=index)\n",
    "    im = ax.pcolormesh(data.longitude, data.latitude, data[VARIABLE],\n",
    "                       transform=ccrs.PlateCarree(),\n",
    "                       norm=colors.LogNorm(),\n",
    "                       vmin=val_min, vmax=val_max,\n",
    "                       cmap=COLOR_PALETTE)\n",
    "\n",
    "    ax.set_title(data.time.dt.strftime('%b %Y').values)\n",
    "    ax.set_extent([data.longitude.min(), data.longitude.max(), \n",
    "                  data.latitude.min(), data.latitude.max()], ccrs.Geodetic())\n",
    "    index_ax = index_ax + 1\n",
    "\n",
    "# Color bar\n",
    "colorbar = fig.colorbar(im, ax=axes.ravel().tolist(), aspect=40)\n",
    "colorbar.locator = mticker.AutoLocator()\n",
    "colorbar.formatter = mticker.ScalarFormatter(useMathText=True, useOffset=True)\n",
    "colorbar.formatter.set_scientific(True)\n",
    "colorbar.formatter.set_powerlimits((19, 19))\n",
    "colorbar.ax.yaxis.set_offset_position('left')\n",
    "colorbar.ax.yaxis.labelpad = 20\n",
    "colorbar.ax.yaxis.set_tick_params(width=1)\n",
    "colorbar.set_label(LABEL)\n",
    "colorbar.update_ticks()\n",
    "\n",
    "plt.savefig(\"{folder_plots}/france_month_global_NO2.{format}\".format(folder_plots=folder_plots,\n",
    "                                                                     format=FORMAT_EXPORT),\n",
    "            bbox_inches='tight', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maps with point of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lyon, Arcelor Mittal Dunkerque, Croisement A4/A35 Vendenheim, Forêt des Landes, Pic du Midi de Bigorre\n",
    "lon_points = [4.8357, 2.3212992, 7.7337068, 0.3380, 0.1411]\n",
    "lat_points = [45.7640, 51.0244429, 48.6716962, 43.6034, 42.9369]\n",
    "\n",
    "LABEL_POINTS = ['U', 'I', 'R', 'F', 'M']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FIGSIZE = 4\n",
    "SCALEBAR_LENGTH = 100000\n",
    "\n",
    "# Moyenne temporelle\n",
    "data = DS.mean(dim='time')\n",
    "\n",
    "# Création du canevas\n",
    "fig = plt.figure(figsize=(FIGSIZE*data.dims['longitude'] / data.dims['latitude'], FIGSIZE))\n",
    "ax = plt.axes(projection=PROJECTION)\n",
    "ax.set_facecolor('white')\n",
    "\n",
    "# Ajout du pcolormesh\n",
    "p = data[VARIABLE].plot(\n",
    "        ax=ax,\n",
    "        transform=ccrs.PlateCarree(),\n",
    "        cmap=COLOR_PALETTE,\n",
    "        #norm=colors.LogNorm()\n",
    "    )\n",
    "\n",
    "# Add points labels\n",
    "droite = 9.5\n",
    "gauche = -6\n",
    "milieu = DS.longitude.mean().values\n",
    "\n",
    "for i,(x,y) in enumerate(zip(lon_points,lat_points)):\n",
    "    if x < milieu:\n",
    "        plt.plot([x, gauche],[y, y], '-', color='black', linewidth=0.3, transform=ccrs.PlateCarree())\n",
    "        plt.text(gauche-0.5, y-0.2, LABEL_POINTS[i], transform=ccrs.PlateCarree(), fontsize=7)\n",
    "    else:\n",
    "        plt.plot([x, droite],[y, y], '-', color='black', linewidth=0.3, transform=ccrs.PlateCarree())\n",
    "        plt.text(droite+0.4, y-0.2, LABEL_POINTS[i], transform=ccrs.PlateCarree(), fontsize=7)\n",
    "\n",
    "# Add point\n",
    "s = sns.scatterplot(lon_points, lat_points, transform=ccrs.PlateCarree(),\n",
    "               color=\".2\", s=10, linewidth=0.5)\n",
    "\n",
    "# Add scale bar\n",
    "scalebar = AnchoredSizeBar(ax.transData,\n",
    "                           SCALEBAR_LENGTH, \n",
    "                           '{length} {unit}'.format(length=int(SCALEBAR_LENGTH/NUMBER_METERS_PER_UNIT), \n",
    "                                                    unit=SCALEBAR_UNIT),\n",
    "                           SCALEBAR_LOCATION, \n",
    "                           pad=0.1,\n",
    "                           fontproperties=fm.FontProperties(size=SCALEBAR_TEXT_SIZE),\n",
    "                           sep=SCALEBAR_SEP,\n",
    "                           color='black',\n",
    "                           label_top=True,\n",
    "                           frameon=False)\n",
    "ax.add_artist(scalebar)\n",
    "\n",
    "# Re-fit axis extent \n",
    "ax.set_extent([data.longitude.min(), data.longitude.max(), data.latitude.min(), data.latitude.max()],\n",
    "              ccrs.Geodetic())\n",
    "\n",
    "# Color bar\n",
    "p.colorbar.locator = mticker.AutoLocator()\n",
    "p.colorbar.formatter = mticker.ScalarFormatter(useMathText=True, useOffset=True)\n",
    "p.colorbar.formatter.set_scientific(True)\n",
    "p.colorbar.ax.yaxis.set_offset_position('left')\n",
    "p.colorbar.ax.yaxis.labelpad = 20\n",
    "p.colorbar.ax.yaxis.set_tick_params(width=1)\n",
    "p.colorbar.update_ticks()\n",
    "\n",
    "# Export\n",
    "plt.savefig(\"{folder_plots}/dotted_france_global_NO2.{format}\".format(folder_plots=folder_plots, format=FORMAT_EXPORT),\n",
    "            bbox_inches='tight', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PI evolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,4))\n",
    "ax = plt.axes()\n",
    "\n",
    "# format the ticks\n",
    "ax.xaxis.set_major_locator(mdates.MonthLocator())\n",
    "ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %Y'))\n",
    "ax.xaxis.set_minor_locator(mdates.DayLocator())\n",
    "ax.yaxis.set_major_formatter(mticker.ScalarFormatter(useMathText=True, useOffset=True))\n",
    "\n",
    "for index,(lon, lat) in enumerate(zip(lon_points, lat_points)):\n",
    "\n",
    "    DS.resample({'time': '1MS'}).mean(dim='time')[VARIABLE].sel(latitude=lat, longitude=lon, method='nearest').plot(\n",
    "        label=LABEL_POINTS[index], ax=ax, linewidth=1)\n",
    "\n",
    "# Change le style pour impression-friendly\n",
    "for l, ms in zip(ax.lines, cycle('ox+>^*')):\n",
    "    l.set_marker(ms)\n",
    "\n",
    "plt.xlabel('')\n",
    "plt.ylabel(LABEL)\n",
    "plt.title('')\n",
    "    \n",
    "plt.legend()\n",
    "plt.xticks(rotation=30)\n",
    "\n",
    "plt.savefig(\"{folder_plots}/evol_points_interet_NO2.{format}\".format(folder_plots=folder_plots, format=FORMAT_EXPORT),\n",
    "            bbox_inches='tight', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}