storm_forecast_starting_kit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"text-align: left\">\n",
    "<img src=\"http://project.inria.fr/saclaycds/files/2017/02/logoUPSayPlusCDS_990.png\" width=\"800px\">\n",
    "</div>\n",
    "\n",
    "# [RAMP](https://www.ramp.studio/problems/storm_forecast_hackathon) on Tropical Storm Intensity Forecast (from reanalysis data)\n",
    "\n",
    "_Sophie Giffard-Roisin (CU/CNRS), Mo Yang (CNRS), Balazs Kegl (CNRS/CDS), Claire Monteleoni (CU/CNRS), Alexandre Boucaud (CNRS/CDS)_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. [Introduction](#Introduction)\n",
    "2. [The prediction task](#The-prediction-task)\n",
    "2. [Installation of libraries](#Installation-of-libraries) : To do before coming!\n",
    "2. [The data](#The-data)\n",
    "3. [The pipeline](#The-pipeline)\n",
    "4. [Evaluation](#Evaluation)\n",
    "5. [Local testing/exploration](#Testing-the-submission)\n",
    "6. [Submission](#Submitting-to-the-online-challenge:-ramp.studio)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "The goal of the RAMP is to predict the intensity of tropical and extra-tropical storms (24h forecast) using information from past storms since 1979. The intensity can be measured as the maximum sustained wind over a period of one minute at 10 meters height. This speed, calculated every 6 hours, is usually explained in knots (1kt=0.514 m/s) and is used to define the hurricane category from the [Saffir-Simpson scale](https://en.wikipedia.org/wiki/Saffir–Simpson_scale). Estimating the intensity evolution of a storm is of course crucial for the population.\n",
    "\n",
    "<img src=\"https://github.com/sophiegif/ramp_kit_storm_forecast_new/blob/master/figures_pynb/all_storms_since1979_IBTrRACKS_newcats.png?raw=true\" width=\"70%\">\n",
    "<div style=\"text-align: center\">Database: tropical/extra-tropical storm tracks since 1979. Dots = initial position, color = maximal storm strength according to the Saffir-Simpson scale.</div>\n",
    "\n",
    "Today, the forecasts (track and intensity) are provided by a numerous number of guidance models (1). Dynamical models solve the physical equations governing motions in the atmosphere. Statistical models, in contrast, are based on historical relationships between storm behavior and various other parameters. However, the lack of improvement in intensity forecasting is attributed to the complexity of tropical systems and an incomplete understanding of factors that affect their development. What is mainly still hard to predict is the rapid intensification of hurricanes: in 1992, Andrew went from tropical depression to a category 5 hurricane in 24h. \n",
    "\n",
    "Machine learning (and deep learning) methods have been only scarcely tested, and there is hope in that it can improve storm forecasts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The prediction task\n",
    "\n",
    "<ul class=\"list-unstyled list-inline text-center\">\n",
    "  <li>\n",
    "    <img src=\"https://github.com/sophiegif/ramp_kit_storm_forecast_new/blob/master/figures_pynb/storm_shema3.png?raw=true\" alt= \"image1\" width=\"350\" height=\"350\">\n",
    "    <figcaption>Goal: estimate the 24h-forecast intensity of all storms.</figcaption>\n",
    "  </li>\n",
    "  \n",
    "  <li>\n",
    "    <img src=\"https://github.com/sophiegif/ramp_kit_storm_forecast_new/blob/master/figures_pynb/hurricane_pb.png?raw=true\" alt= \"image2\" width=\"350\" height=\"350\">\n",
    "    <figcaption>Feature data: centered maps of wind, altitude, sst, slp, humidity...</figcaption>\n",
    "  </li>\n",
    "</ul>\n",
    "\n",
    "This challenge proposes to design the best algorithm to predict for a large number of storms the 24h-forecast intensity every 6 hours. The (real) database is composed of more than 3000 extra-tropical and tropical storm tracks, and it also provides the intensity and some local physical information at each timestep (2). Moreover, we also provide some 700-hPa and 1000-hPa feature maps of the neighborhood of the storm (from ERA-interm reanalysis database (3)), that can be viewed as images centered on the current storm location (see right image).\n",
    "\n",
    "The goal is to provide for each time step of each storm (total number of instants = 90 000), the predicted 24h-forecast intensity, so 4 time steps in the future. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "References\n",
    "\n",
    "1. National Hurricane Center Forecast Verification website, https://www.nhc.noaa.gov/verification/, updated 04 April 2017.\n",
    "\n",
    "2. Knapp, K. R., M. C. Kruk, D. H. Levinson, H. J. Diamond, and C. J. Neumann, 2010: The International Best Track Archive for Climate Stewardship (IBTrACS): Unifying tropical cyclone best track data. Bulletin of the American Meteorological Society, 91, 363-376 https://www.ncdc.noaa.gov/ibtracs/index.php?name=wmo-data\n",
    "\n",
    "3. Dee, D. P. et al.(2011), The ERA-Interim reanalysis: configuration and performance of the data assimilation system. Q.J.R. Meteorol. Soc., 137: 553–597. https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.828"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Installation of libraries\n",
    "\n",
    "To get this notebook running and test your models locally using the `ramp_test_submission`, we recommend that you use the Python distribution from [Anaconda](https://www.anaconda.com/download/) or [Miniconda](https://docs.anaconda.com/docs_oss/conda/install/quick#miniconda-quick-install-requirements). (uncomment the lines before running them)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# !conda install -y -c conda conda-env     # First install conda-env to ease the creation of virtual envs in conda\n",
    "# !conda env create                        # Uses the local environment.yml to create the 'storm_forecast_2' env"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**OR** if you have Python already installed but are **not using Anaconda**, you'll want to use `pip` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# !pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Installation of ramp-workflow\n",
    "\n",
    "For being able to test submissions, you also need to have the `ramp-workflow` package locally. You can install the latest version with pip from github:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# !pip install git+https://github.com/paris-saclay-cds/ramp-workflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Download data (optional)\n",
    "\n",
    "If the data has not yet been downloaded locally, uncomment the following cell and run it.\n",
    "The starting kit data is 260 MB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# !python download_data.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The data\n",
    "\n",
    "The 3000 storms have been separated in a train set, a test set and a local starting kit (train+test sets). The data from `download_data.py` (local starting kit) includes only 1/4 storms of the total database; and the train set on which your code will run on the platform has another half. They are disjoined. \n",
    "\n",
    "Let's have a look at the local train data (only the first rows are plotted)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>stormid</th>\n",
       "      <th>instant_t</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>hemisphere</th>\n",
       "      <th>Jday_predictor</th>\n",
       "      <th>initial_max_wind</th>\n",
       "      <th>max_wind_change_12h</th>\n",
       "      <th>basin</th>\n",
       "      <th>...</th>\n",
       "      <th>vo700_10_1</th>\n",
       "      <th>vo700_10_2</th>\n",
       "      <th>vo700_10_3</th>\n",
       "      <th>vo700_10_4</th>\n",
       "      <th>vo700_10_5</th>\n",
       "      <th>vo700_10_6</th>\n",
       "      <th>vo700_10_7</th>\n",
       "      <th>vo700_10_8</th>\n",
       "      <th>vo700_10_9</th>\n",
       "      <th>vo700_10_10</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1979152N11265</td>\n",
       "      <td>0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>-95.5</td>\n",
       "      <td>25.0</td>\n",
       "      <td>1</td>\n",
       "      <td>5.433485e-08</td>\n",
       "      <td>25.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.000010</td>\n",
       "      <td>-1.940000e-05</td>\n",
       "      <td>-0.000010</td>\n",
       "      <td>-0.000012</td>\n",
       "      <td>-0.000002</td>\n",
       "      <td>-5.100000e-06</td>\n",
       "      <td>-0.000002</td>\n",
       "      <td>-0.000014</td>\n",
       "      <td>0.000012</td>\n",
       "      <td>0.000036</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1979152N11265</td>\n",
       "      <td>1</td>\n",
       "      <td>11.0</td>\n",
       "      <td>-96.2</td>\n",
       "      <td>30.0</td>\n",
       "      <td>1</td>\n",
       "      <td>5.433485e-08</td>\n",
       "      <td>25.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>3</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000014</td>\n",
       "      <td>-6.500000e-06</td>\n",
       "      <td>-0.000020</td>\n",
       "      <td>-0.000015</td>\n",
       "      <td>-0.000003</td>\n",
       "      <td>-3.000000e-07</td>\n",
       "      <td>-0.000006</td>\n",
       "      <td>-0.000010</td>\n",
       "      <td>0.000027</td>\n",
       "      <td>0.000026</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1979152N11265</td>\n",
       "      <td>2</td>\n",
       "      <td>11.0</td>\n",
       "      <td>-96.9</td>\n",
       "      <td>30.0</td>\n",
       "      <td>1</td>\n",
       "      <td>5.433485e-08</td>\n",
       "      <td>25.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>3</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000006</td>\n",
       "      <td>-9.000000e-07</td>\n",
       "      <td>-0.000022</td>\n",
       "      <td>-0.000044</td>\n",
       "      <td>-0.000069</td>\n",
       "      <td>-2.370000e-05</td>\n",
       "      <td>0.000018</td>\n",
       "      <td>0.000012</td>\n",
       "      <td>0.000031</td>\n",
       "      <td>0.000022</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1979152N11265</td>\n",
       "      <td>3</td>\n",
       "      <td>11.3</td>\n",
       "      <td>-97.6</td>\n",
       "      <td>30.0</td>\n",
       "      <td>1</td>\n",
       "      <td>5.433485e-08</td>\n",
       "      <td>25.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.000002</td>\n",
       "      <td>-6.500000e-06</td>\n",
       "      <td>-0.000026</td>\n",
       "      <td>-0.000040</td>\n",
       "      <td>-0.000018</td>\n",
       "      <td>-8.300000e-05</td>\n",
       "      <td>-0.000036</td>\n",
       "      <td>0.000034</td>\n",
       "      <td>0.000035</td>\n",
       "      <td>0.000024</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1979152N11265</td>\n",
       "      <td>4</td>\n",
       "      <td>11.8</td>\n",
       "      <td>-98.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>1</td>\n",
       "      <td>5.433485e-08</td>\n",
       "      <td>25.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.000031</td>\n",
       "      <td>-4.720000e-05</td>\n",
       "      <td>-0.000025</td>\n",
       "      <td>0.000015</td>\n",
       "      <td>-0.000006</td>\n",
       "      <td>-2.390000e-05</td>\n",
       "      <td>0.000004</td>\n",
       "      <td>0.000025</td>\n",
       "      <td>0.000018</td>\n",
       "      <td>0.000016</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 859 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         stormid  instant_t  latitude  longitude  windspeed  hemisphere  \\\n",
       "0  1979152N11265          0      11.0      -95.5       25.0           1   \n",
       "1  1979152N11265          1      11.0      -96.2       30.0           1   \n",
       "2  1979152N11265          2      11.0      -96.9       30.0           1   \n",
       "3  1979152N11265          3      11.3      -97.6       30.0           1   \n",
       "4  1979152N11265          4      11.8      -98.0       30.0           1   \n",
       "\n",
       "   Jday_predictor  initial_max_wind  max_wind_change_12h  basin     ...       \\\n",
       "0    5.433485e-08              25.0                  0.0      3     ...        \n",
       "1    5.433485e-08              25.0                  5.0      3     ...        \n",
       "2    5.433485e-08              25.0                  5.0      3     ...        \n",
       "3    5.433485e-08              25.0                  0.0      3     ...        \n",
       "4    5.433485e-08              25.0                  0.0      3     ...        \n",
       "\n",
       "   vo700_10_1    vo700_10_2  vo700_10_3  vo700_10_4  vo700_10_5    vo700_10_6  \\\n",
       "0   -0.000010 -1.940000e-05   -0.000010   -0.000012   -0.000002 -5.100000e-06   \n",
       "1    0.000014 -6.500000e-06   -0.000020   -0.000015   -0.000003 -3.000000e-07   \n",
       "2    0.000006 -9.000000e-07   -0.000022   -0.000044   -0.000069 -2.370000e-05   \n",
       "3   -0.000002 -6.500000e-06   -0.000026   -0.000040   -0.000018 -8.300000e-05   \n",
       "4   -0.000031 -4.720000e-05   -0.000025    0.000015   -0.000006 -2.390000e-05   \n",
       "\n",
       "   vo700_10_7  vo700_10_8  vo700_10_9  vo700_10_10  \n",
       "0   -0.000002   -0.000014    0.000012     0.000036  \n",
       "1   -0.000006   -0.000010    0.000027     0.000026  \n",
       "2    0.000018    0.000012    0.000031     0.000022  \n",
       "3   -0.000036    0.000034    0.000035     0.000024  \n",
       "4    0.000004    0.000025    0.000018     0.000016  \n",
       "\n",
       "[5 rows x 859 columns]"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from problem import get_train_data\n",
    "import numpy as np\n",
    "\n",
    "data_train, y_train = get_train_data()\n",
    "data_train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[30. 30. 35. ... 45.  0.  0.]\n"
     ]
    }
   ],
   "source": [
    "print(y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the data is a list of time instants (one every 6h). The first storm will result in x lines beginning with its stormid and the corresponding time step, with all the associated features on the same row. Then the time steps from the second storm will be below, and so on. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of storms in the local training set: 587\n"
     ]
    }
   ],
   "source": [
    "print('Number of storms in the local training set: {}'.format( len(set(data_train['stormid'])) ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of time steps in the local training set: 15777\n"
     ]
    }
   ],
   "source": [
    "print('Total number of time steps in the local training set: {}'.format(y_train.size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. 0D features from track data\n",
    "\n",
    "A set of simple features has been extracted for each storm at each time point: \n",
    "\n",
    "- latitude, longitude: in degrees\n",
    "- windspeed: current (max) windspeed (knots) \n",
    "- hemisphere:  South=0, North=1\n",
    "- Jday predictor:  Gaussian function of (Julian day of storm init - peak day of the hurricane season), see (1)\n",
    "- initial_max_wind: initial (max) windspeed of the storm \n",
    "- max_wind_change_12h: last 12h (max) windspeed change\n",
    "- basin = based on the present location: \n",
    "       0 = NA - North Atlantic / 1 = SA - South Atlantic    / 2 = WP - West Pacific       / 3 = EP - East Pacific /\n",
    "       4 = SP - South Pacific  / 5 = NI - North Indian      / 6 = SI - South Indian       / 7 = AS - Arabian Sea /\n",
    "       8 = BB - Bay of Bengal  / 9 = EA - Eastern Australia / 10 = WA - Western Australia / 11 = CP - Central Pacific\n",
    "       12 = CS - Carribbean Sea/ 13 = GM - Gulf of Mexico   / 14 = MM - Missing\n",
    "- nature = nature of the storm  \n",
    "       0 = TS - Tropical / 1 = SS - Subtropical / 2 = ET - Extratropical / 3 = DS - Disturbance /\n",
    "       4 = MX - Mix of conflicting reports / 5 = NR - Not Reported / 6 = MM - Missing / 7 =  - Missing\n",
    "- dist2land = current distance to the land (km)\n",
    "\n",
    "\n",
    "(1) DeMaria, Mark, et al. \"Further improvements to the statistical hurricane intensity prediction scheme (SHIPS).\" Weather and Forecasting 20.4 (2005): 531-543. https://journals.ametsoc.org/doi/full/10.1175/WAF862.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. The reanalysis data\n",
    "\n",
    "At each time step, we extracted 7 grids (11x11 pixels) of meteorological parameters centered on the current storm location. Their choice is based on the forecast literature, on personal experience and on known hypothesis of storm strengthening.\n",
    "\n",
    "#### a) 25x25 degree z, u and v at 700hPa-level\n",
    "First, we provide 3 maps of 25 x 25 degrees (lat/long) at 700hPa-level pressure: the altitude `z`, the u-wind `u` (positive if wind from the West) and the v-wind `v` (positive if wind from the South). These grids are subsampled to 11x11 pixels (1 pixel ~=2 degrees).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>z_0_0</th>\n",
       "      <th>z_0_1</th>\n",
       "      <th>z_0_2</th>\n",
       "      <th>z_0_3</th>\n",
       "      <th>z_0_4</th>\n",
       "      <th>z_0_5</th>\n",
       "      <th>z_0_6</th>\n",
       "      <th>z_0_7</th>\n",
       "      <th>z_0_8</th>\n",
       "      <th>z_0_9</th>\n",
       "      <th>...</th>\n",
       "      <th>z_10_1</th>\n",
       "      <th>z_10_2</th>\n",
       "      <th>z_10_3</th>\n",
       "      <th>z_10_4</th>\n",
       "      <th>z_10_5</th>\n",
       "      <th>z_10_6</th>\n",
       "      <th>z_10_7</th>\n",
       "      <th>z_10_8</th>\n",
       "      <th>z_10_9</th>\n",
       "      <th>z_10_10</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>114352.8</td>\n",
       "      <td>114374.9</td>\n",
       "      <td>114387.9</td>\n",
       "      <td>114395.1</td>\n",
       "      <td>114428.3</td>\n",
       "      <td>114492.2</td>\n",
       "      <td>114562.4</td>\n",
       "      <td>114578.7</td>\n",
       "      <td>114575.3</td>\n",
       "      <td>114582.1</td>\n",
       "      <td>...</td>\n",
       "      <td>114340.3</td>\n",
       "      <td>114356.6</td>\n",
       "      <td>114393.2</td>\n",
       "      <td>114436.4</td>\n",
       "      <td>114441.7</td>\n",
       "      <td>114387.4</td>\n",
       "      <td>114270.6</td>\n",
       "      <td>114041.3</td>\n",
       "      <td>113682.3</td>\n",
       "      <td>113547.7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>114196.2</td>\n",
       "      <td>114197.1</td>\n",
       "      <td>114204.9</td>\n",
       "      <td>114229.1</td>\n",
       "      <td>114275.5</td>\n",
       "      <td>114302.6</td>\n",
       "      <td>114324.4</td>\n",
       "      <td>114364.1</td>\n",
       "      <td>114400.3</td>\n",
       "      <td>114420.2</td>\n",
       "      <td>...</td>\n",
       "      <td>114373.7</td>\n",
       "      <td>114401.3</td>\n",
       "      <td>114456.9</td>\n",
       "      <td>114455.0</td>\n",
       "      <td>114469.5</td>\n",
       "      <td>114479.2</td>\n",
       "      <td>114438.1</td>\n",
       "      <td>114261.5</td>\n",
       "      <td>113938.3</td>\n",
       "      <td>113724.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>114387.3</td>\n",
       "      <td>114379.0</td>\n",
       "      <td>114383.9</td>\n",
       "      <td>114408.6</td>\n",
       "      <td>114428.9</td>\n",
       "      <td>114442.9</td>\n",
       "      <td>114464.2</td>\n",
       "      <td>114461.3</td>\n",
       "      <td>114460.8</td>\n",
       "      <td>114501.9</td>\n",
       "      <td>...</td>\n",
       "      <td>114362.1</td>\n",
       "      <td>114459.4</td>\n",
       "      <td>114578.4</td>\n",
       "      <td>114585.1</td>\n",
       "      <td>114592.9</td>\n",
       "      <td>114589.5</td>\n",
       "      <td>114533.4</td>\n",
       "      <td>114384.9</td>\n",
       "      <td>114194.7</td>\n",
       "      <td>113915.6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>114082.5</td>\n",
       "      <td>114092.6</td>\n",
       "      <td>114107.6</td>\n",
       "      <td>114120.7</td>\n",
       "      <td>114146.3</td>\n",
       "      <td>114157.5</td>\n",
       "      <td>114176.3</td>\n",
       "      <td>114185.5</td>\n",
       "      <td>114224.7</td>\n",
       "      <td>114269.2</td>\n",
       "      <td>...</td>\n",
       "      <td>113964.0</td>\n",
       "      <td>114174.4</td>\n",
       "      <td>114301.6</td>\n",
       "      <td>114371.8</td>\n",
       "      <td>114423.6</td>\n",
       "      <td>114402.8</td>\n",
       "      <td>114417.3</td>\n",
       "      <td>114357.3</td>\n",
       "      <td>114222.8</td>\n",
       "      <td>114046.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>114382.4</td>\n",
       "      <td>114369.9</td>\n",
       "      <td>114343.7</td>\n",
       "      <td>114330.7</td>\n",
       "      <td>114311.3</td>\n",
       "      <td>114307.4</td>\n",
       "      <td>114314.7</td>\n",
       "      <td>114335.0</td>\n",
       "      <td>114380.5</td>\n",
       "      <td>114415.8</td>\n",
       "      <td>...</td>\n",
       "      <td>113963.0</td>\n",
       "      <td>114211.7</td>\n",
       "      <td>114433.7</td>\n",
       "      <td>114525.6</td>\n",
       "      <td>114513.1</td>\n",
       "      <td>114502.9</td>\n",
       "      <td>114455.5</td>\n",
       "      <td>114354.9</td>\n",
       "      <td>114249.9</td>\n",
       "      <td>114101.8</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 121 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      z_0_0     z_0_1     z_0_2     z_0_3     z_0_4     z_0_5     z_0_6  \\\n",
       "0  114352.8  114374.9  114387.9  114395.1  114428.3  114492.2  114562.4   \n",
       "1  114196.2  114197.1  114204.9  114229.1  114275.5  114302.6  114324.4   \n",
       "2  114387.3  114379.0  114383.9  114408.6  114428.9  114442.9  114464.2   \n",
       "3  114082.5  114092.6  114107.6  114120.7  114146.3  114157.5  114176.3   \n",
       "4  114382.4  114369.9  114343.7  114330.7  114311.3  114307.4  114314.7   \n",
       "\n",
       "      z_0_7     z_0_8     z_0_9    ...       z_10_1    z_10_2    z_10_3  \\\n",
       "0  114578.7  114575.3  114582.1    ...     114340.3  114356.6  114393.2   \n",
       "1  114364.1  114400.3  114420.2    ...     114373.7  114401.3  114456.9   \n",
       "2  114461.3  114460.8  114501.9    ...     114362.1  114459.4  114578.4   \n",
       "3  114185.5  114224.7  114269.2    ...     113964.0  114174.4  114301.6   \n",
       "4  114335.0  114380.5  114415.8    ...     113963.0  114211.7  114433.7   \n",
       "\n",
       "     z_10_4    z_10_5    z_10_6    z_10_7    z_10_8    z_10_9   z_10_10  \n",
       "0  114436.4  114441.7  114387.4  114270.6  114041.3  113682.3  113547.7  \n",
       "1  114455.0  114469.5  114479.2  114438.1  114261.5  113938.3  113724.5  \n",
       "2  114585.1  114592.9  114589.5  114533.4  114384.9  114194.7  113915.6  \n",
       "3  114371.8  114423.6  114402.8  114417.3  114357.3  114222.8  114046.2  \n",
       "4  114525.6  114513.1  114502.9  114455.5  114354.9  114249.9  114101.8  \n",
       "\n",
       "[5 rows x 121 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_train_z = data_train[[col for col in data_train.columns if col.startswith('z_')]]\n",
    "data_train_z.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "sample_id=20 # sample number plotted - you can change it to see other storms and other instants\n",
    "grid_l=11 # size of all 2D-grids (in pixels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAEICAYAAABoNzG1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu4HFWZ7/HvTwIkQCAJAUMIEIeLykURwsUDKspFUEccRrmMCAwyoCPiHMFR1FE84mW8HHFEQVAGERGRA4ooFxWDgCBGJ0KAgBEiCUmAkAQIEJKQ9/yx1ia1O71vq3f17mT/Ps+zn93dVW+tVZe36+2q6i5FBGZmZmZWj5cMdQfMzMzM1mUutszMzMxq5GLLzMzMrEYutszMzMxq5GLLzMzMrEYutszMzMxq5GJrLSPpBEm31jBdSfpvSYsl3TnY0x9gX2qZxybtHCBpbhvaCUk71N3OcCZpcl7OI4a6L9adpKmSTuph2MclfWcQ25ot6aDBml6nGor3FEnnS/qPdra5LnGxVZET9TlJSyt/5w51v9pkf+BgYFJE7N04ML8pVpfLc5JWSRqfh28o6SJJT0laIOnDDfEHSpop6VlJv5G0XWlH85v3styPhZKukrRV6fSsO0kfkTRD0tOSHpL0kYbhjXlyYz+n+1ZJt0pakreRCyWNrgy/WNLyhu1svcGev+FM0lmSLh3qfnSJiM9HRNNCrC95ezl7sPvUqe22U7MPvBHxvoj47FD1aW3nYmtNfx8Rm1T+Th3qDrXJdsDsiHim2cD8pvjicgH+E5gaEQvzKGcBO+bpvBH4d0mHAuSC7CrgP4BxwDTgRy3299Tcj52AMcDXWpyerSbgOGAscChwqqSjG8ap5skh/ZzuZsDZwETglcAk4MsN43ypIf9eKJ+N1uUjvn6fzHzksHX+ADFMRYT/8h8wGzioh2HnAVdWnv8n8GvSjmkscC3wOLA4P55UGXcqaSfzO2Ap8DNgc+AHwFPAH4DJlfEDOA14EFhI2iG9JA87Abi1Mu4rgF8Ci4D7gSN7mb+JwDV53FnAv+TX3wssA17I/ftMH8tJwF+B4yuvPQIcUnn+WeDy/Phk4HeVYRsDzwGv6GH6JwC3Al/Jy/Mh4LCG5XlS5fkHgBn58VuB/8nLdQ5wVi/zcQAwt2H5/L+8Hh8CTqu8/hwwrjLua/K6WT8/PxG4L/f3BmC7hvW5Qw99GOi28fU8X08BfwReVxl2FnAlqZB9GvgT8OrK8I/m9fR03lYO7Gde/BfwjX7myS+Ar1ae/wi4qIdxjwDurjy/GDi7n31aL28fC0l58oG8nEfk4ZsB3wXm53k+G1ivEvvVHPsQcGpD7FTgc8Bteb3v0Nv0elv/pFz5GvAY8CRwF7BrP+dxF1bn9qPAx/PrLwE+RsrBJ4ArurZNYHKel+OBh/M8fiIPOxRYDqzI29qf+7GsTsjL4Wu5H2f3Y3s/GJiZ5/dc4GYq+dowj2cBl/bV9yZxJ+f5WJ7n5WeVbfOMvJyfJG1/IytxbwOmA0tIOfeqXt7j1lhvvbT7yrzdLAHuAd7esF2fR8qNZ4CD8mvfAq7L07kNmACck5fpTOA1vWwbL76n5PV3Cel962/AJ8n7izz8X/K6ehq4F9gjv961DXW9/g+VeanuD5Y0y8883Vl5u7gGmNjQv/cBf8nz801A/dnu19W/Ie9AJ/3R+05kI+AB0pvP6/IbwaQ8bHPgH/M4o4EfAz+pxE7NG+X2OTHuzdM6CBiRE+W/K+MH8BvSUaBt87gn5WEnkIstUtEyB/jnPJ09cr926WEebs4JPhLYPSfngY3T7cdyen1Owk3y87G5zy+tjPNO8o6UVCCc1zCNGcA/9jD9E0hvaP9C2jG+H5jXlaxUii1gPHAT8P38/ABgN9IO6VWkndQ7emjnAHKxlcf/I/ApYAPg70g78Tfn4TeRi9P8/MvA+fnxO/L6fWVeD5+ke3HZV7E1kG3jWNL2NgI4HVhA3pmQdlwr8rJfn7TTeSg/fnneVibmcScD2/djXYtUvL6vIU8ezdvPjXQv6CaQdlBvAt6dl+HoHqZ9Drkgz88vJr1xL8rroun2kcd9H2mHtA0pT35D94LpJ8C3STmyJXAncEol9l7SkbWxwK9Ys9h6mFTsjMjLr7fp9bj+gTfneRmTl+Urga36sdxHk4qf00n5OhrYJw/7N+CO3P8Nc79+WFmvAVwIjAJeDTwPvLKyjVza0FZv83YCsBL4YJ63UX3M73jSB4GubfB/5/iBFFtN+94k9mIainPStnkn6QPSOFKR8b48bA/StrkP6X3l+Dz+hk2m3eN6a2w3z+cs4OOk9443kQqYl1fGfxLYj/Q+MzK/thDYMz+/iZSrx+W+nQ38ppfto1psXQL8NG8jk0nvH+/Nw95FKqD3yvOxA6s/CLwrL6eXAEeRCsGueTyBhv1Bdb7zPC7My3RD4BvAbxv6d21eftuS3isO7Wu7X5f/hrwDnfSXE28p6dNJ1191B7s3aUfwN+CYXqazO7C48nwqlU9opE/V11We/z0wvfI8qhsm8K/Ar/PjF5MgJ8gtDW1/G/h0kz5tQ/qkMrry2heAixun24/l9N2uuMq0g+6fIA8mnZbsGv+LDdO4DTihh+mfAMyqPN8oT39CZXk+m9fPI6SjQFv0MK1zgK/1MOwAVhdb+wAPNww/k1zoACcBN+XHIhUur8/PryO/ueXnL8n9266yPnsrtvq9bTSJX0wudkg7rjsa+jGf9OFgB9KO5iDy0bh+ruvPAH+mskMi7TRG5fVyJqngG1MZfkRePguB/XuY7sG57ztVXtuD1YXkW0g7rP16iL+J7gXgIXk5jwBeStpJj6oMP4a888qxp1SGHcSaxdb/qQzva3o9rn/STukBYF8qRxv6sdyPAf6nh2H3UTkqCWxFKrJHsLpgqR5ZvxM4urKNXDqAeTuBNfOit/k9rmEbFDCXgRVbTfveJPZimhdbx1aef4nVH4rOAz7bMP79wBuaTLvH9dbYLim/FtD9aNIPyUfV8/iXNJnGhZXnHwTuqzzfjXxEqYd5D1JOr5fX386VYaeQLvGAdNTxQ/3c5qYDh1fWe2/F1ndJp/y7hm2St8HJlf7tXxl+BfCx/m7/6+Kfr0VY0zsiYkzl78KuARFxJ+mTukgbDwCSNpL0bUl/k/QU8FtgTMO5+Ucrj59r8nyThn7MqTz+G+kTSKPtgH3yBcdLJC0hHU2Y0GTcicCiiHi6YbpbNxm3R5JGkT4Rfa/y8tL8f9PKa5uSdpZdw6vDXhwu6XWVi6HvqQxf0PUgIp7ND6vL6LS8fraOiHdHxOO5f/vkC/Afl/Qk6SjG+H7M2nbAxIZl+XHSzgjS6bnXSppIOrIXwC2V2K9X4haRtpH+Ltt+bxuSTpd0n6Qnc1ubNczfi9tNRKwi7egmRsQs0hGRs4DHJF2e56VHkk4l7TzfGhHPV6Z7W0Q8FxHPRsQXSEXv6yqh15J2AvdHxBrfKpW0L3AZ8M6IeKAy3T9FxBMRsTIifkEqoo/ooXsTWTNHumxHOtowv7JOvk06atMstvq42Wt9Ta/H9R8RN5FOpX0TeFTSBZIac6GZbUineJrZDri60t59pA9SL62Ms6Dy+FnWfH/p77zBmsunt+2927KNtKdttnx709++DzR+O+D0hhzfhibvrQNcbxOBOTnfujS+tzZbBgPdJzQznnQ0rbr9V9vucTuSdJyk6ZVlsSv9e6+ENM8vthkRS0mntKvz3Op6XKe42BoASR8gHTKdB/x7ZdDppNM0+0TEpqSdMaQ3oFLbVB5vm9tsNAe4uaE43CQi3t9k3HnAuOq3v/J0Hxlgv44gvblO7XohIhaTjqC8ujLeq0nXLpD/vzhM0sak02b3RMQtsfpi6F0G2JdmLiNdP7BNRGwGnE//1sMc4KGGZTk6It4CEBFLSKfMjgT+iXTaJiqxpzTEjoqI3w3C/LxI0utI110dCYyNiDGk0xPV+dumMv5LSKea5uV5uCwi9iftdIJ03WFPbZ1IuqbjwIjo6+cxoqEPnyMVAFtJOqZhuq8hrZ8TI+LXA5xu1XzWzJEuc0if9sdX1semle1rPmm5dKlOp9p2f6fX6/qPiP+KiD1JpyV3Arp9u7MHc0g50tOwwxraGxkR/cnlaHje17z1FNPT/HZbL5JE8+U7GBr71Zc5wOca+r1RRPyw6cR7Xm+N7c4Dtmn4IkXje+tA+9pfC0lHlLbroe2m21H+NviFpOsVN8/vJTNYnW999Xdetc38nr45A9+fDBsutvpJ0k6k8+jHAu8hfdtu9zx4NOmTyBJJ44BPD0KTH5E0VtI2wIdo/u29a4GdJL1H0vr5by9Jr2wcMSLmkC4I/YKkkZJeRbow/gcD7NfxpEPijcl4CfDJ3OdXkK63ujgPuxrYVdI/ShpJui7qroiYOcC2+2M06QjeMkl7kwqj/rgTeErSRyWNkrSepF0l7VUZ5zLSkZ5/zI+7nA+cKWkXAEmbSXpX67OyhtGk618eB0ZI+hRrHjHcU9IR+Vtj/0bakd4h6eWS3iRpQ9LFr8+RjoasQdK7gc8DB0fEgw3DtpW0n6QN8nb0EdKn4dvy8NeTriE8Lv99Q9LWediuwPXAByPiZ03afaekTSS9RNIhpFy7podlcQVwmqRJksaSCkMAImI+qTD+qqRN8/S2l/SGSuyHJG0taQypgO1RP6bX4/rP+biPpPVJ18R0XXjc9fX62T00ey0wQdK/Kf2symhJ+1Ta+1zeYSJpC0mH9zYPFY8Ck7sKg37MWzO9be8/B3apbIOn0fxI+2B4lHRtZX9dCLwvrw9J2ljp50hGN47Y23pr0u7v8zj/nt+DDyCd/r984LM0MJG+rXsFaXsYnbeJDwNdP+/xHeAMSXvmed4hj7MxqaDqOiPwz6QjW10eBSZJ2qCHpi8D/lnS7vk95fPA7yNi9iDP4jrDxdaafqbuv/NzdX7TuBT4z4j4c0T8hXSK6ft5QzuHdA3LQtKFq9cPQj9+SrpAczrpDey7jSPkU4KHAEeTPmksIB2t2LCHaR5Dui5iHqkA+nRE/LK/Hco7zTeRCqtGnyYdrv4b6UL8L0fE9bmfj5MKlM+RrtPZJ/e5Dv8K/B9JT5OKuiv6GJ/cxxdIb5C7ky5UXUh6o9qsMto1pJ+3eDQi/lyJvZq03C9XOo08Azis9VlZww2k62UeIC3nZax5euKnpGv5FpM+FBwREStI28QXSfO1gHSa6OM9tHM26VPqHyp5cH4eNpp07cti0qfYQ0lHWZ5QOs1yCelnOR7JpxC/C/x3PsJxOrAF8F01P3X8oTzNJaQvIPxLREztoY8X5uXxZ9K3Lq9qGH4c6fTKvbmvV5KubeqKvZH0DbP/IX1LbCU9FJ99Ta+P9b9pbm8xaZ09QfoWJaQjPrc1ayzn9sGkbXIB6Vtdb8yDv07aFm/M2/kdpJzqjx/n/09I+lNf89ZD33qc30g/BfMu0rb2BClfms7jIPgusLPSabCf9DVyREwjfQg8lzSfs0jXJjXT23rr1m5ELAfeTloGC0lfQjqupg+TzXyQVOw9SPoW92XARQAR8WPS++5lpMs6fkL65uq9pOtDbycVVrvRfT3dRDojsUDSQhrko9L/Qfr29nzS0bN+vacrXzoy4Llcy3V9u8s6iKQAdszX2Zj1i6SzSBfiHzvUfVmbSDqMdBH1dn2OPLjt3ki6ePm+drZrZu3nI1tmNqzk08RvkTQiH639NOlIb1tFxCEutMyGBxdbZjbciPSTFotJpxHvI51yNjOrhU8jmpmZmdXIR7bMzMzMauRiy2wISPqhpHf0c9zZkg7qZfh1ko4fxL5NlXTSYE1vsEn6uKTvFMYeIKmv3w0bFJIulnR2O9qqtLlU0kB+DqEa++J6l/R2SbX/dIHZcOFiy6zNlH7j7NWkn2kYaOxZki6tvhYRh0XE9/LwEySt8avt65KI+HxEdGwxOJTyjwM/2PeYfU7nGtJv471qELplNuy52DJrv1OAHzT5YVizTvJD4OSh7oTZusDFlln7HUb64VcA8i923yTpCUkLJf0g/7J5N5IOJf0Q6VH5dNGf8+tTJZ2kdOeA80n3cFyqdL+zNU4LNh79knSwpJlK91s8l4Zb5Eg6Uel+jIsl3ZB/gbpPjUfhJE2WFPlHgpuN/zdJe+bHx+Zxd87PT+r64crqdCvTPF7Sw3n5faIyzVH5dN5iSfcCezVputqHr0uaI+kpSX9UukVSdX6ukHSJpKcl3SNpSmX4ayT9KQ/7ETCyhzY2zD+IuWvltS0kPSdpyybj/7Okn1Wez5JUvTfrHOW7WeRlsUN+fLGkb0r6ee7T7yVtX4nrdb2Tbsn11t6Wl5n1j4stszZSuofYy4D7qy8DXyDd3PWVpF8WP6sxNv8i/+eBH+XTRa9uGH4f6cbbt+fhaxRsTfoznvQr0J8k3Xbnr8B+leHvIBV4R5B+/f0W0hGPOtwMHJAfv570i9hvqDy/uUlMl/1J9yc9EPiUVt+y6tOkX7feHngz6XZTvfkD6S4C40i/uv1jpVtMdXk76TYsY0i/4n4ugNJtTX4CfD/H/ph014Q1RLqp91WkOzp0OZJ0n9PHmoTcDLxO6VY6W5FuHL1fbvfvSDf4vauH+TmG9DMXY0m/mP65HNfres/uI93apz83zjazXrjYMmuvrgLo6a4XImJWRPwyIp7Ptzb6v6wuMur2FuDeiLgy39bnHNLtYbqcAnwhIu6LiJWkYm/3/h7dGqCbWT3fryMVoF3P30DvxdZnIuK5fBulP7P6xudHkm4+vCjfH/S/eutARFwaEU9ExMqI+CrpNkcvr4xya0T8It/e6fuVdvYlFUHnRMSKiLiSVLj15DK6F1v/RPf7bVb79CBpe9mdtBxuAB5RugfpG4BbImJVD+1cFRF35nX3gzwN6Hu9w+pttM+i3cx652LLrL2W5P8v3vxW0paSLpf0iNK95i4lHW1oh4lU7q+YryOr3m9xO+Dr+bTXEmAR6Ujc1q00qnx/NHW/P2LXEZwJwHqkm6/vJ2ky6R6V03uZZLVQeJZ0tAca5o90n7ve+nV6PmX6ZJ7fzei+LhrbGZlPi04EHmm4Dq+3tm4CRind7Hg7UhF0de7DdZVl8+48ftdRv64jfFNJhVZfRWi/lkuT9Q6rt9ElmFlLXGyZtVFEPEM6ZbNT5eUvAAG8KiI2BY5lzetnXpxEX000ee0ZYKPK8wmVx/NJpy0BkKTqc9IO+JSIGFP5GxURv+ujH722GxG35FOdm0TELvm1WaSC4DTgt/lmzAtIF2nf2svRm950mz9g255GzNdnfZR0NGxsPg37JD2vi8Z2ts7Lr8+28rxcQTq69U/AtXl+u75d2rVsfpBDuoqt1+XHXUcB+yq2eutvb+sd0int2RHxVMH0zazCxZZZ+/2C7qcJRwNLgSVK9+r7SC+xj5Kuo+kpdx8FJuVriLpMB46QtFG+ePq9lWE/B3aRdEQ+QnMa3Yux84EzJe0CIGkzSe/qGqh08f1ZPfRlOvB6SdtK2gw4s5f56nIzcCqrC4ipDc8H6gpS/8dKmgR8sJdxRwMrgceBEZI+BfT3eqXbc+xpSvdcPALYu4+Yy4CjgHfTwynEipuBNwKjImIu6dq5Q4HNSbccGqi+1jukbfS6gmmbWQMXW2btdwHw7spRkM8Ae5COovycdPF0T36c/z8h6U9Nht8E3AMskLQwv/Y1YDmpEPse6dodACJiIfAu4IvAE8COwG2V4VcD/wlcnk9xziB9m7LLNtXxqyLil6RTgXcBfwSu7WW+utxMKnp+28PzgfoM6XTeQ8CNpOusenIDqbh4IMcsY81Ta01FxHLSlwhOIN1z8Sh6X49ExO9JR/8m0kdRExEPkAryW/Lzp0hfILgtXz82IH2t9+wY4NsDnbaZrcn3RjQbApIuA66IiJ8MdV9K5SNFP46I1w51X2xwSfp74D0RceRQ98VsXeBiy8zMzKxGPo1oZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1GjHUHajadNPxseWWk4tin3uuvN2nn24l9oXyYJa2ENuKjYsjpfJNJmJVcWxry+qxwrhniViuFhpu2fgxY2LyhAllwSrv+rL1yreRxx8vDmX99ctjJ4x6sjj2+VmzimNXFEfCBi3EtqJ0y5gDLIoY0pwYNWp8bLbZ5KLYRYsGty/9tWLF88Wx48ZtWBz7shFzimNb2anGjjsVxz74YHEoS5Y8Ux5cvI9ZRMQzfeZERxVbW245mS99aVpR7IwZ5e1OnVoe+5vfLC4P5o4WYluxZ3HkiBFbFseuWNFKItzWQuw3CuNuaaHNwTF5wgSmXXBBWfDIkcXtPjBm7+LYc88tDmXrrctjP7rrz4tjH3zb24pj5xVHwqQWYlt58y49pXFYC20Ols02m8xxx5XtJy6/vLzdlSvLY+fPf6g49s1vfllx7GXjTyuOZebM4tDl195YHHv00cWhXH11K/vU0n3MOf0ay6cRzczMzGrkYsvMzMysRi62zMzMzGrkYsvMzMysRi62zMzMzGrkYsvMzMysRi62zMzMzGrkYsvMzMysRoNSbEm6SNJjkmZUXhsn6ZeS/pL/jx2MtszWBs4Js+6cEzacDdaRrYuBQxte+xjw64jYEfh1fm42XFyMc8Ks6mKcEzZMDUqxFRG/BRrvOnU48L38+HvAOwajLbO1gXPCrDvnhA1ndV6z9dKImA+Q/ze9qZ6kkyVNkzTtySdbuIOtWecbcE48vmRJWzto1mYDzonnnvN+wtY+Q34j6oi4ALgAYMruu8cRB5Tdln3y5HHFfVi2rDiU228vv8Rg2bIdi2N3222H4th3v7s4lN12K4+FjYsjZ848pDj2yivLYu+6a0pxm63olhO77BKMGVM2ofHji/uwbGFxKLvvXh574it+Vx78gU8Wh25b3irLW4h9rIXYuS3E3lUYN1RlTjUnNthgSlx6adl05s9vZam1Elu+n7jsjD8Vx9615zeKY3ctjoQNrr2qOPaoo44ojl25ct/i2GnTymIXLvxhv8ar88jWo5K2Asj/W3lfMVsXOCfMunNO2LBQZ7F1DXB8fnw88NMa2zJbGzgnzLpzTtiwMFg//fBD4Hbg5ZLmSnov8EXgYEl/AQ7Oz82GBeeEWXfOCRvOBuWarYg4podBBw7G9M3WNs4Js+6cEzac+RfkzczMzGrkYsvMzMysRi62zMzMzGrkYsvMzMysRi62zMzMzGrkYsvMzMysRi62zMzMzGrkYsvMzMysRi62zMzMzGrkYsvMzMysRoNyu55Bs2gRXHppUegexx5b3OzsKeOKYw8+uDiUrbfeoTj2/e8vb/dVK/9UHnz77eWxS5YUh75liy2KYz/8sQlFcVM+XN7fQbN8OcyZUxa7Q/n2tXBmcSjvfGd5LCedUx57zz3FofPKW2WDFmJb2cIWtBC7ZWHc+i20OVhWrHiG+fPvKIxe3ELLjxVH7rbb8X2P1JP9RhWHLipvlVktxO5UuB8HOOqTk4tjD7t0j+LYGTPK4k48sX/j+ciWmZmZWY1cbJmZmZnVyMWWmZmZWY1cbJmZmZnVyMWWmZmZWY1cbJmZmZnVyMWWmZmZWY1cbJmZmZnVyMWWmZmZWY1cbJmZmZnVyMWWmZmZWY1cbJmZmZnVyMWWmZmZWY1cbJmZmZnVaMRQd6Cbp5+GqVPLYidPLm52ypS3F8eOH18cyuv3X1UefO655bHnnVccunTmzOLYJcWRsLKF2FLLh6DNNay/Pmy1VVHoNb/aqLjZBQuKQ3kTN5UHjxlTHrv55sWh277nPeXtXn55cejIOXOKY7ctjiz/lP2tFtocPI9R3pPdWmh3y+LIu84pz4lbD1xWHDu1OBJOayGWCROKQxdN3qM49vpri0OZPbss7umn+zde7cWWpNnA08ALwMqImFJ3m2adzDlh1p1zwtZ17Tqy9caIWNimtszWBs4Js+6cE7bO8jVbZmZmZjVqR7EVwI2S/ijp5MaBkk6WNE3StMeff74N3TEbcv3PicWLh6B7Zm3X75wA7yds7dOO04j7RcQ8SVsCv5Q0MyJ+2zUwIi4ALgCYMnZstKE/ZkOt/zmx887OCRsO+p0T0ubOCVvr1H5kKyLm5f+PAVcDe9fdplknc06YdeecsHVdrcWWpI0lje56DBwCzKizTbNO5pww6845YcNB3acRXwpcLamrrcsi4vqa2zTrZM4Js+6cE7bOq7XYiogHgVfX2YbZ2sQ5Ydadc8KGA//0g5mZmVmNXGyZmZmZ1cjFlpmZmVmNXGyZmZmZ1cjFlpmZmVmN2nUj6v555hmYNq0s9oADipvddv9F5bE7LCuO5azzy2O//OXi0OnLyvv8cHEklC9lWNVCbKmnhqDNRo8/sxHfumOPoti77y5v9/3vL4+F8eWh73tfeewZZ5THLlhQHjtpUnHoxN13H5JYri/7ZYWRZ55Z3uageQEovY1V+e2vxo79SHEsl55YHDq7vFX+VwuxS1uIHTdhQnnssnnFseutN7E49jvfKYtb2M9bp/vIlpmZmVmNXGyZmZmZ1cjFlpmZmVmNXGyZmZmZ1cjFlpmZmVmNXGyZmZmZ1cjFlpmZmVmNXGyZmZmZ1cjFlpmZmVmNXGyZmZmZ1cjFlpmZmVmNXGyZmZmZ1cjFlpmZmVmNXGyZmZmZ1WjEUHegm1Wr4Omny2IXLChvd9q08tiZM8tjL7ywOPTOZcuKY+8tjoRFLcSW9xiWtxC7NrXZ6NFH4ZxzymLf+97ydl81qYU1PX1heewmmxSH3jh7p+LYSZPKY3/F64tjv3h0cWhLb3mTJh1Z2OaXyhsdNBsCOxTGvry41bPPLg6F6eW72n1baLaFTGzN9tsXhz47ZmJx7I9+VBzK/PlCom2KAAAOrUlEQVRlcStW9G88H9kyMzMzq5GLLTMzM7MaudgyMzMzq5GLLTMzM7Ma1V5sSTpU0v2SZkn6WN3tmXU654TZas4HGw5qLbYkrQd8EzgM2Bk4RtLOdbZp1smcE2arOR9suKj7yNbewKyIeDAilgOXA4fX3KZZJ3NOmK3mfLBhoe5ia2tgTuX53PzaiySdLGmapGmPr1pVc3fMhtyAcuKFFx5va+fM2qzPfIDuOQHPta1zZoOl7mJLTV6Lbk8iLoiIKRExZYuX+Hp9W+cNKCfWW2+LNnXLbEj0mQ/QPSdgVBu6ZTa46q5u5gLbVJ5PAubV3KZZJ3NOmK3mfLBhoe5i6w/AjpJeJmkD4GjgmprbNOtkzgmz1ZwPNizUem/EiFgp6VTgBmA94KKIuKfONs06mXPCbDXngw0Xtd+IOiJ+Afyi7nbM1hbOCbPVnA82HPiKdDMzM7Ma1X5ka0BeeAGWLCmLnTmzvN25c8tjp04tDn14wYLi2FauIC1vFVr5cY6h2thGFsatN6i9KLNqFSxdWhZ78MEtNHzGGeWxCxeWxx5wQHHopEP3Lo496KDiUObP/3Z5MLe1ELu8OHLOnMMKI58obnPwbAzsWRT52tceW9zqv064qjiW8eOLQ3fYbbfi2HF3310cO6Y4Ethrr+LQ888vb3ZECzuZD36wLO6SS/o3no9smZmZmdXIxZaZmZlZjVxsmZmZmdXIxZaZmZlZjVxsmZmZmdXIxZaZmZlZjVxsmZmZmdXIxZaZmZlZjVxsmZmZmdXIxZaZmZlZjVxsmZmZmdXIxZaZmZlZjVxsmZmZmdXIxZaZmZlZjUYMdQeqlgNzC2M3ufrq4nZXFkfCYy3EPthC7MMtxC5vIXZcC7HjW4gdCusPdQeAnXaCq64qjL31ovKGFywoj122rDz2jjuKQ3d+2wPFsfPnb1EcC9NbiF3cQuy+LcTuXhi3UQttDo7Ro8ex777HFsXe+J0W3jkvnlEeu2RJeeykScWh4044obzdyZOLQ2+cvVNx7IgWqpJjyzYLAMaMKYu75pr+jecjW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmVqPaii1JZ0l6RNL0/PeWutoyWxs4J8xWcz7YcFL3L8h/LSK+UnMbZmsT54TZas4HGxZ8GtHMzMysRnUXW6dKukvSRZLGNhtB0smSpkmatqjmzph1gAHlxOLFj7e7f2bt1Gc+QPecWLHCOWFrn5aKLUm/kjSjyd/hwHnA9qQ7ns4HvtpsGhFxQURMiYgprdzk2KwTDHZOjB3byg2SzYbWYOQDdM+J9dd3Ttjap6VrtiLioP6MJ+lC4NpW2jJbGzgnzFZzPpgltV0gL2mriJifn/4DMKOvmOeA6YXttXKIblULsUtbiG3ltOlTLcRu0kLspBZiJ7YQu7KF2NLlvF4LbTZTkhMjR6xkpzGPlTW4SQtr+oADymMnTCiP3XXX4tBF43cqb5cnWojdq4XYnVuI3a848o1vfHVR3LRpo4rbbFSSD5A2rzPOKGx0yZLCwBZjp00rj91//+LQh9/54eLYK68sDmX27PLYV7yiPPbt+7ewV124sChskxHL+jVend9G/JKk3YEAZgOn1NiW2drAOWG2mvPBho3aiq2IeE9d0zZbGzknzFZzPthw4p9+MDMzM6uRiy0zMzOzGrnYMjMzM6uRiy0zMzOzGrnYMjMzM6uRiy0zMzOzGrnYMjMzM6uRiy0zMzOzGrnYMjMzM6uRiy0zMzOzGrnYMjMzM6uRiy0zMzOzGtV2I+oSy4GHC2NbqRpXtRC7rIXYVmzQQuyWLcRObCF2TAuxS1qILV1WaqHNwbJk6QiuuaNsjY2ZcGRxu7OWFodyx63lsZtML4+dNKk8duTIzYtjX3jhxOLYU08tDmWffcpj//rXsrj77y9vc7BsOmoFh+w6ryz46lvKG25lAzv22PLYWbOKQ6dOLW/29NOfKY597Ws3Lo5tJSe47rry2Ne8poWG++YjW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmViMXW2ZmZmY1crFlZmZmVqOWii1J75J0j6RVkqY0DDtT0ixJ90t6c2vdNFs7OCfMunNOmMGIFuNnAEcA366+KGln4GhgF2Ai8CtJO0XECy22Z9bpnBNm3TknbNhrqdiKiPsAJDUOOhy4PCKeBx6SNAvYG7i91+kBqwr7srIwjhbabDW2FRsNUWwry3lJC7FLW4htp8HOifnz4eyzy/oyfXpZHMCKFY+VBzO3OHLHHfcojj3wwOJQxo4tj/3KV8pj/+no8neQSy4tPzFx5ZVlcYsXDzxmsHOClSth4cKBdwRgm23K4gD23784dN6yccWxrWxfXzv+mfJgyt8DvvKVlxXH7nTrRcWx7Lprceiflu1cFPfsqpH9Gq+ua7a2BuZUns/Nr61B0smSpkmatrbsUM0KFOXEypWPt6VzZkOgKCceL6n4zIZYn0e2JP0KmNBk0Cci4qc9hTV5LZqNGBEXABcAbCs1Hcesk7QzJzbeeIpzwjpeO3Niyi67OCdsrdNnsRURBxVMdy5QPV47CZhXMB2zjuOcMOvOOWHWu7pOI14DHC1pQ0kvA3YE7qypLbO1gXPCrDvnhA0brf70wz9Imgu8Fvi5pBsAIuIe4ArgXuB64AP+hokNB84Js+6cE2atfxvxauDqHoZ9DvhcK9M3W9s4J8y6c06Y+RfkzczMzGrlYsvMzMysRi62zMzMzGrkYsvMzMysRi62zMzMzGrkYsvMzMysRi62zMzMzGrkYsvMzMysRoronHt6Snoc+FtB6Hhg4SB3p4T70V2n9APK+rJdRGxRR2f6yzkxaNyPNTknhob70V2n9ANqzImOKrZKSZoWEVPcD/ejJ53Ul3bolPl1PzqzH9BZfWmHTplf96Mz+wH19sWnEc3MzMxq5GLLzMzMrEbrSrF1wVB3IHM/uuuUfkBn9aUdOmV+3Y/uOqUf0Fl9aYdOmV/3o7tO6QfU2Jd14potMzMzs061rhzZMjMzM+tILrbMzMzMarTWFluS3iXpHkmrJE1pGHampFmS7pf05jb36yxJj0ianv/e0ub2D83zPUvSx9rZdkM/Zku6Oy+DaW1s9yJJj0maUXltnKRfSvpL/j+2Xf1pJ+dEj+07J5wTzonu7Tsn2pwTa22xBcwAjgB+W31R0s7A0cAuwKHAtySt1+a+fS0ids9/v2hXo3k+vwkcBuwMHJOXx1B5Y14G7fwNlYtJ673qY8CvI2JH4Nf5+brIOdHAOQE4J5wTFc4JYAhyYq0ttiLivoi4v8mgw4HLI+L5iHgImAXs3d7eDZm9gVkR8WBELAcuJy2PYSMifgssanj5cOB7+fH3gHe0tVNt4pxoyjnhnHBOdOecGIKcWGuLrV5sDcypPJ+bX2unUyXdlQ9VtvPwfCfMe5cAbpT0R0knD1Efurw0IuYD5P9bDnF/2q0TtgvnhHOik3TCduGcGEY5MWIwJzbYJP0KmNBk0Cci4qc9hTV5bVB/36K3fgHnAZ/NbX4W+Cpw4mC231vXmrw2VL/tsV9EzJO0JfBLSTPzpwlrgXNi4F1r8ppzYh3inBh415q85pyoWUcXWxFxUEHYXGCbyvNJwLzB6VHS335JuhC4djDb7kPt895fETEv/39M0tWkQ9dDlUSPStoqIuZL2gp4bIj60TLnxIA5J5pzTjgnujgnklpzYl08jXgNcLSkDSW9DNgRuLNdjeeV1OUfSBdotssfgB0lvUzSBqQLQK9pY/sASNpY0uiux8AhtHc5NLoGOD4/Ph7o6dPuuso54Zxo5JxwTjgnuqs3JyJirfwjbaBzgeeBR4EbKsM+AfwVuB84rM39+j5wN3BXXnlbtbn9twAP5Pn/xBCtm78D/pz/7mlnP4AfAvOBFXn7eC+wOenbJX/J/8cNxXJpw7w7J5q375xwTjgnurfvnGhzTvh2PWZmZmY1WhdPI5qZmZl1DBdbZmZmZjVysWVmZmZWIxdbZmZmZjVysWVmZmZWIxdbZmZmZjVysbUOkNTRdwIwayfng9manBdDy8VWB5A0WdJMSd/LNya9UtJGedinJP1B0gxJF0hSfn2qpM9Luhn4kKS/l/R7Sf8j6VeSXprHOytP90ZJsyUdIelLku6WdL2k9Rv6MlHS9MrfC5K2a/tCsWGrk/KhMu0p+fF4SbPbtzTMkk7KC0mvlHRnQ9/uauPiWOu42OocLwcuiIhXAU8B/5pfPzci9oqIXYFRwNsqMWMi4g0R8VXgVmDfiHgNcDnw75XxtgfeChwOXAr8JiJ2A57Lr78oIuZFxO4RsTtwIfD/IuJvgz2zZn3oiHww6zAdkRcRcR+wgaS/yy8dBVwxiPO5znGx1TnmRMRt+fGlwP758RvzJ5G7gTcBu1RiflR5PAm4IY/3kYbxrouIFaTbQ6wHXJ9fvxuY3KwzkvYDTqJ9d6I3q+qofDDrEJ2UF1cAR+bHRzW0Yw1cbHWOxvsmhaSRwLeAd+ZPGBcCIyvjPFN5/A3Sp5vdgFMaxnseICJWASti9T2aVgFrnMdXuknqd4GjImJp+SyZFeuYfABWsvq9cmST4Wbt0kl58SPgSEk7pbD4S+E8DQsutjrHtpJemx8fQzrc25UICyVtAryzl/jNgEfy4+N7Ga9X+dz8FcBHI+KB0umYtagj8iGbDeyZH/fWplndOiYvIuKvwAvAf+CjWn1ysdU57gOOzxcZjgPOi4glpE8pdwM/Af7QS/xZwI8l3QIsbKEf/wvYC/hM5SL5iS1Mz6xEp+QDwFeA90v6HTC+xWmZtaKT8gJSkXUsvl6rT1p9pNCGiqTJwLX54kazYc35YLYm58XazUe2zMzMzGrkI1tmZmZmNfKRLTMzM7MaudgyMzMzq5GLLTMzM7MaudgyMzMzq5GLLTMzM7Ma/X85wQ067J4z0AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "params_25x25=['z','u','v']\n",
    "plt.figure(figsize=(10,4))\n",
    "for p,param in enumerate(params_25x25):\n",
    "    image=np.zeros([grid_l,grid_l])\n",
    "    for i in range(grid_l):\n",
    "         for j in range(grid_l):\n",
    "            image[i,j]=data_train[param+'_'+str(i)+'_'+str(j)][sample_id]\n",
    "    plt.subplot(1,3,p+1)\n",
    "    plt.imshow(np.array(image),extent=[-12,12,-12,12],\n",
    "               interpolation='nearest', origin='lower', cmap='seismic')\n",
    "    plt.xlabel('param '+param)\n",
    "t=plt.suptitle('Example of 700-hPa level maps 25x25 degrees, centered in the storm location.'\n",
    "         +'\\n (altitude, u-wind and v-wind)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### b) 11x11 degree sst, slp, humidity at 1000hPa,  and vorticity at 700hPa\n",
    "We provide some more localized maps of 11 x 11 degrees (lat/long) at the surface: the sea surface temperature `sst`, surface level pressure `slp`, the relative humidity `hum` at 1000hPa (near surface). These grids are sampled to 11x11 pixels (1 pixel = 1 degree). We also provide the vorticity at 700hPa `vo700`. \n",
    "\n",
    "NB: `sst` is only defined on the sea, so land has NaNs values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAloAAADRCAYAAAAUsxzYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xm8FNWZ//HPA4ioICCoiGxxj5rRKC5xGR23qDGaMa6Jjsa4TeIvOokxi1mIS8zmxEyYaNRE3Pdx10RNRIOJC+4bLlEURVQQEBQE5Pn9cU5D0Xbfe6q76va93O/79bqv29311Dmnq5+uPn2quo65OyIiIiJSvB6tboCIiIjI8kodLREREZGSqKMlIiIiUhJ1tERERERKoo6WiIiISEnU0RIREREpiTpa0jQzO9LMJpRQrpnZRWY208weKrr8VjGzcWZ2RqvbIZ2fmY0yMzezXnWWP2NmOxdU185m9noRZXVmZe2vEuqda2brdHS90nrqaHVyZjbZzObFN2nlb2yr29VBdgB2B4a5+9bVC81sLTO72cymxg+jUVXLDzKzv5vZB2Y2Pk/FZna6mT1lZovMbEyeemX5Fd+Pu7W6HRXuvom7j29k3Zi76xXcpE5bb0cys/FmdnT2MXfv6+4vt6pN0jrqaHUNn49v0srfCa1uUAcZCUx29/frLF8M/An4Yp3l7wLnAD9roO6XgFOA2xqotyXqjXpI5xBHaLXPbYJyXLoivem7MDM718yuy9z/uZn9Je7QB5rZrWb2Tjz0dquZDcvEjjezM+KIz1wzu8XMBpnZ5Wb2npk9nB2pid9Cv2FmL5vZdDP7Zb0PDTPbyMzuMrN3zex5MzuojecwNI4OvWtmL5nZMfHxrwIXAp+J7ftJ9bru/pa7/w54uFbZ7n63u18DTK1R78Hxuawa7+9lZtPMbPW47sXufgcwJ2+9Ner6tJk9amZzzOxqoE/V8n3M7HEzmxVfj3/JLNvCzB6L615rZldXDjtWDvWY2XfMbBpwUUJ5Q83s+pgXr5jZNzLLtjazifH1f8vM/rvO86nUe4qZvW1mb5rZF8xsbzN7Ib6W368q9x+xPW+a2Vgz651ZXje3zGw9M7vXzGbHZVenbPO47jFm9lzcds+a2RYJ22CMmV1jZpfE9Z4xs9Fx2aXACOCWmJOnxMe3jdt5lpk9YZlDefF9dqaZ3Q98AKxjZv3N7A9xW7wR34c9Y3xPM/tVfK4vA59r5zkuGWFrq+011rsv3nwiPpeDM8u+lXldv5J5fMXYttdifpxnZivVKb/m61av3vhavRRz52YzG5opy83s62b2IvBi5rGvmdmL8bmebmbrxjx7L26H3tXtqtPW7Szs72bH/9tllq1m4fSFqRb2ozfGx+vuX83sTGBHYKxljkBYZiQv5sAlcf1XzewHmZw/0swmxG09M+boXinPRTopd9dfJ/4DJgO71Vm2MvACcCThjT2dcJgNYBBhxGVloB9wLXBjZt3xhFGbdYH+wLOxrN2AXsAlwEWZeAfuAVYjfNi8ABwdlx0JTIi3VwGmAF+J5WwR27VJnedwL/A7Qudjc+AdYNfqctvZRr1i+0bVWX40ML7G45cD4+K2mgrsUyPmMmBMI/XGmN7Aq8B/ASsABwALgTPi8i2At4FtgJ7AEfE1XzGz7olx3f2BBZl1dwYWAT+P8Su1U14P4BHgR7HsdYCXgc/G8v4BHB5v9wW2rfOcKvX+KLbrmPi6XRFzbRNgPrBOjN8S2DZur1HAc8BJibl1JXBqbHsfYIfE982BwBvAVoAB6xFGSNvbBmNi2/eO2+8s4IF670dgbWBGjO9BONQ9A1g98z57LW6TXnF73Qj8nvBeWQN4CDguxh8PTAKGx+1xT9w+vdrbP7TX9hrrOrBejdf1tNjOvQmdw4Fx+TnAzbFd/YBbgLPqlF33datR7y6EfcQWhDz9LXBfVfxdsd6VMo/dDKwat+2HwF/i61nZnx1Rp21HsnR/tRowEzg8vj6HxvuD4vLbgKuBgXGb7JRj/3p0ve1N2L/eFNcdRcj5r2bat5DwvuoJ/Cdh/2Qpua+/zvfX8gbor50XKOxI5wKzMn/HZJZvTThE9ipwaBvlbA7MzNwfD5yauX82cEfm/ueBxzP3Hdgzc/9rwF/i7eyO62Dgb1V1/x74cY02DQc+AvplHjsLGFddbjvbqNGO1gDCh+BTwO/rrNtsR+tfq3eSwN9Z2lk6Fzi9ap3ngZ3ium9UrTuBZTtaC4A+meVtlbcN8FrVsu8RO9TAfcBPgMHtbO+dgXlAz3i/X9wO22RiHgG+UGf9k4AbEnPrEuB84heIHO+bPwMn1ni8vW0wBrg7s2xjYF7V+zHb0foOcGmNuo/IvM9Oyyxbk9ApWCnz2KHAPfH2X4HjM8v2IF9Hq27ba6xbq6M1L1sXodO+LaGz+j6wbmbZZ4BX6pRd93WrUe8fgF9k7vcldDRGZeJ3qVHG9lX59p3M/bOBc+q07UiW7q8OBx6qWv6PGLMW4TSBgQn5Vmv/WrOjReg8fQhsnFl2HHEfFet+KbNs5bjukDzvAf11nj8dOuwavuDuAzJ/F1QWuPtDhG/kBlxTedzMVjaz38dh6fcIH6IDKocoorcyt+fVuN+3qh1TMrdfBYbycSOBbeJhlFlmNgv4MjCkRuxQ4F13zx6ee5UwSlA6d59F+Ca6KWHHXIahwBse95jRq5nbI4FvVW2v4XG9WutmXwOAd9x9fmJ5I4GhVcu+T/jwB/gqsAEwKR5C2aeN5zXD3T+Kt+fF/zXzx8w2iIdWpsVc/CkwuKq8erl1CiG3H4qHwo5qo01Zw4F/1ni8vW0AMC1z+wOgj9U/N2gkcGBVeTsQPqRrPbeRhJGRNzPxvyeMbEF43tXbIo88ba9lhrsvqiqjL7A64QP/kUy7/xQfryXP6zaUzPN097mEUcHsfqA67yH//qvduqPKPmg4Yf80s3qlxP1rPYNZOlpdXWfFktfR3T+IN1Oej3RC6mh1cWb2dcJw+1TCzq3iW8CGhFGGVQmjIxB2fo0anrk9ghrnPhF2iPdWdQz7uvt/1oidCqxmZv2qyn2jiTYmM7PNgaMIhzn+p6Rq3gTWNrPsdh+RuT0FOLNqe63s7lfWWTf7GkD4ppvVVnlTCCMQ2WX93H1vAHd/0d0PJXzo/xy4zsxWaXoLhFG2ScD6MRe/z8fzsGZuufs0dz/G3YcSvvX/ztJ+sTaFcFi81uN1t0GCWtv70qryVnH3n9VZZwphNGNwJn5Vd98kLn+Tj2+LzmA6ofOySabd/d295od/ztdtKqEDCkDMuUEsux+o3u5FWabuqLIPmkLYPw2osV57+9e22judMGKXrbfD9nvS8dTR6sLMbAPgDOAwwhD4KbHzAOFwzjxglpmtBvy4gCq/HU8CHU44b6jWicm3AhuY2eFmtkL828rMPlkd6O5TCIfRzjKzPhZO2v4q4dypJGbWh9DRBFgx3q8s6xnv9wJ6xDpWyKx3GeFD/yuEDs3XMuuuEGN6AL3iuj0zy+vWW+UfhPNevmFmvcxsf8Lh3ooLgOPNbBsLVjGzz8XO5z8Ih1ZPiOvuV7VuLW2V9xDwnoWT51eK22dTM9sqPqfDzGx1d19MOERNrL9Z/YD3gLlmthHhnJNqNXPLzA60pT/imEn4APsoLhtvVZfeyLgQONnMtozbYT0zG0k72yDBW4TzgCouAz5vZp+t5JuFHwsMq7Wyu78J3AmcbWarmlkPCydx7xRDriHkyjAzGwh8N7Fdjah+LnXFnLgA+LWZrQFgZmub2Wdrxbf1utWo9wrgK2a2uZmtSBjxfNDdJ+d8Po24nbC/+lJ8jx1MOOR6a3yt7iB0EgfGfUKlQ9Xe/rXuto0jwdcAZ5pZv5iX3yTkkiyH1NHqGiq/cqr83RAPB1wG/Nzdn3D3FwmdhkvjzuocwsnR04EHCMP8zbqJcC7E44STRP9QHRAPA+4BHEL4tjiNpSdr13Io4WTQqcANhHO57srRpnmEc9ggjJrMyyw7PN4/l/BjgXmEDwsI54K97u7nuvuHhM7qGWa2flx+QYw/lHBS77xYXkq9S7j7AsJJ7EcSPnAOBv4vs3wi4aTXsXH5SzE2u+5XCR2fwwgd2Q/rbYx2yvuIcO7d5sArhNy4kHDyMMCewDNmNhf4DXBI1WHJRp0MfInwC84LqN1Br5dbWwEPxjbdTDjv6pW4bDhwf60K3f1a4EzCh/gcwgnoqyVsg/acBfzAwuGzk+OXhf0I7713CKMg36btfet/EA4dPUt4ja5j6aHGCwjneD0BPEomV0owBrg4Ppe6vwzO+A4hnx6Ih8vuJozq1NLW67ZMve7+F+CHwPWEEb11CfuP0rn7DGAfwgjVDMJRgX3cfXoMOZww+jSJcL7aSfHx9vavvwEOsPCrwVqj5f+PcM7by4TzLq8A/pjSZgu/9jwv6QlKp2DLnv4hUpuZOeHQz0utbkt3ZmYPAue5+0WtbktRGsmtOFpyrbt/pryWiYg0TyNaIp2Yme1kZkPiYY0jgH+hmNHJLs3dX1cnS0S6Al1lV6Rz25BwPkdfwq/oDojnjoiISBegQ4ciIiIiJdGhQxEREZGSqKPVCZjZlWb2hZLKPsPCXGPT2o/unuLP/G+xMNfZta1ujwTxMgmvd0A9oyzMQ1f6qRRmNiL+crjmhS0tzFd4WUpsK5nZOItzbnZwvePN7OiOrrdoZvZ9M7swIe48M/thg3Xsa2ZXNbKuFEsdrRazcO2ozQg/by+67OGEny1v7O61rsxeHb9kJ9/NHEC4Mvggdz+w1Y3pbMxsLQsT/U6NHZJRVctXNLM/WpjMd5qZfbNq+a5mNsnMPjCze+J1g5LWbadd48xsQeyMvGthIvONmn2+ZXL31+IFfNu9Pll1bFfpZMRORPZyNPPMbLGZDY7LG86XduqtdJgr9U42szKvQ5ak1hcGd/+pu7f7Wrr78e5+er1y2ln3ZmBTy0wqL62hjlbrHQdc7gWfLBe/nY8kTKnxdpFlL0/iaMFI4IWqqUeKKr/lPzgpoA2LCb90/GKd5WOA9Qnb8d8IF87dM9Y9mHAtqB8SJvCdyLLX0aq7bqJfxKuTDyNc52hcjnWlBLET0bfyR7iO3vjMtanG0Hi+pBgQ6z0U+FHOfCpUJ3j/Xwkc2+I2iHeCCRe78x/hgnXZme3XA+4FZhMuhnd1fHwUVZPLkpm4lHBRyvuBXxMmmZ5AuIjmYsKFNce10449CRMUL4zxT8TH+xMuHvkmYYqIM1g6mXC2zlnxuWwXH59C+OA7IlPHOOA84C7CRSTvBUY2uN32JlzwcU5s18mZNk2oil0yiW1sw7mEK0K/H9uffd5fTah7POHClQ/F1+kmwsUws6/TVwkTVt8XH9+WcBX8WYSLUe6cKe/IuO3mEC6i+eWCc6EyCfVRwHOEi2T+Oe+2p84k2nH775G5fzpwVbx9LPD3zLJVCHm5UcK6OwOvE0Zl3ybk4Feq8umMzP3PAXPj7a0JV9afFdcbC/Su87wq2/OI+JpNZ9kJ16vr2ZlwsdvK/cmEi5Q+GXPqD4QR0jvia3o3cWLi6tcO+ER8jecQ3hdjgcuqYwkXX/0ImE/I07HA/wJnVz2XW4CT6jzP3xDel+8RLg67Y2bZGMKvWy+JbXkGGJ1Z/mnCBVTnEDo+V2W3SRs5Y4Rfy2b3A83ky/gYf39sy53ESdCrt2187GGW7hvqPv+qNm9LuNByz8xj/w48GW9XLgg9Nf6dA6xYlbPfiWVcy7L74bmE+RXHVF7nuN4OLN0/TAGOzOZeZjtUl/MBYSS+Us6WhAvnrhDvb0+dib/113F/GtFqIQtzen0CeD7z8OmEncdAwrf03+YochvCB/YawO7AXsBUD98sj2xrRXf/E2Hqi6tj/GZx0cWEKWTWI+xs9wCyQ97bED5gBhGubnwV4arQ6xGuZD7WzLLzoX05PsfBhKuAJ0+3U+UPwHHu3o8wKfRfc6z7JcIHVz9gV5Z93h+72n0d/0HouAwlbJ/qqz/vBHwS+KyZrU242vkZhG/pJwPXm9nqMQf+B9grPpftCNsFisuFM+M5gN8nXGl+deBvhG+7TbEwTcxQQuex4gmgMnffJtll7v4+4YN3k4R1IUxG3p8w4e5Xgf+N61W3oy8htx6LD30E/Bchzz5DeJ2/Vr1elR0Il9PYlTAS8rFpo9rwRcJ7bgPClefvIGzvwYQjB9+os94VhA/9wYTX+4haQe5+KuE1OyHm6QmE9+ahZtYDlowG7Ur91/VhwhXxV4v1XmvLTh21L+H9O4BwNfexsdzehCvrXxrXvZb6o5vVdiR0Oq+PZTWcL5n4LxGmzVqDcIX9k6srtWD7uF4lJ9p7/pU6HyB0mHepqvOKePtUQmdsc8JpH1sDP8jEDol1jCTsJ7L74b7uvswcsWY2gpAvvyW8Nzdn6T4gux1qlTMeyF7V/zBCp3VhvP8cMMrMVq1+ntJx1NFqrcpkpXMyj1UmGx3q7vPdfUKO8qa6+2/dfZG715wSJg8zW5Pw5j7J3d/3cAjy1yw7PcYr7n6Rh/NIriZMi3Kau3/o7ncSRouyk8ne5u73eZj25lTgM/FcsrwWAhub2aruPtPdH82x7k3ufr+7L/bGp5i51N2fjjvAHwIHVZ20PCZus3mEnd/t7n57rPMuwiGRykTGiwnnUqzk7m+6+zOZ51hULhwHnOXuz3k4RPpTYPPU81/aUOlEz848NpvQia0sn82yKsvbWxfCNjjN3Re6++2Eb/LZaV9ONrNZhKlh+rJ0uqFH3P2B+PwnA78ndH7b8hN3n+fuTxA+7DdrJz7rt+7+lru/QegQPejuj8U8v4HwJWUZ8QN2K+CH8f1yH2FEKom7V0ZUd40PHUI4RPdWnfjL3H1G3CZnE0ZmsttyQszRjwidqsrz3xZYATgnvg7XETotKY4ArnP3ynRVzeRLxUXu/kLM62sIHZOs6YSR3AuB73qY4ifl+WddSTj0iIV5QvdmaQf2y4ScfNvd3wF+wrLTcy0mTCX2YeJ++MvA3e5+Zdy+M9z98XbXCi4m7F8qp0EcSnjtKiqfLbUmxpYOoo5Wa1Um7s3uRE4hDLc/ZGbPmNlROcqbUljLgpGEHeybFuYlm0X4wFojE5Pdqc8DqNrRz2PpznWZNsad77uEb7h5fZGw83vVzO41szxXCS9iO2XLeJWwnQbXWT4SOLCyDeN23AFYK3bUDgaOJ2zn2zIndBeZCyOB32TqfzeWvXaOMmupfIBmvzGvytId/NyqZdnl7a0L4RzD7LlzH7BsPv3K3Qe4+xB339fd/wlhwnUzuzWebP0eoWOZfX1qyf4yt7qe9lTnfFvvgYqhwMyYAxWv5qgTMh+08f+l9QLN7Ftm9pyFX9fOIowUZrdJ9fPvE88xGgq84e7Z80jbbaeZrQQcGNtY0Uy+1Gtn9bYd7O4D3f2T7r5kpDnh+WddAexvYd7Y/YFH3b3ynIey7PN/lWX3Ye/k/AI3nDBq14ibCF841yGMqM6OHfCKymfLrI+tKR1GHa0WygyLb5B5bJq7H+PuQwmjEL8zs/UIQ9kAK2eKqP4lYbMn1FevP4UwgfHg+GE2wN1XdfdNaqybasnoVTzcsxrhPId8DXV/2N33I3T6biR8s4WwnZZsIzOr9WvLIn54kB2FG0EYeZmeeSxbxxTCCNiAzN8q7v4zAHf/s7vvTphYeBJx4uuCc2EK4VBrtg0rufvfcz/zbCXuMwnnQGVHfzYjnOND/L9kWTxUui7wTMK6zTiXsC3Xd/dVCYfxrMGylskpPr6tG/UmMDBuk4oRbcTXytvLgP3MbDPCoeoba61oZjsSzhs6iHC+2ADCSFHKNnkTWNvMsrFttbNif0KHfvySJ9BEviTUV1fe5+/uzxI6UHux7GFDCPur7EjwCJbdh1W/Tu3tb6YQnmN7PlZO7NBdQxgVO5yPd7Q/CUx29/cSypeSqKPVereTOaRhZgdamDAXwknLDnwUh6jfAA4zs55xdCPlzZnHW4Tj+T0APEz1cidwtpmtamY9zGxdM2vvEExb9jazHeJ5H6cTDrHUHGGKP9Xeucbjvc3sy2bW38O5CO8RzsmBeL6HmW0ez78Y00Rb23KYmW1sZisDpxEOj9T7yf5lwOfN7LPxtesTf6o9zMzWtHC9m1UIndq5ledScC6cB3zPzDaJZfc3syWXsrBw6YAx9VaO23LFeHfFqnNbLgF+YGYD42jcMSz99d8NhMOiX4zr/IhwUvGkhHWb0Y+QF3Njuf/ZRFmPE/J2tdhxP6mA9hFHSCYCP4k5vQPh/K563gLWqSrjdcJhvEuB69s4VNWPcC7hO0AvM/sRHx85qucfcd1vWJhzc3/CeUntOQK4pGokDJrLl0Y18vyvIJxb96+E89IqroztXz2eF/cjwnu8nreAQWbWv87yy4HdzOyguH0HmVn14dC2yrmEcMh83xrt2Ilw/pe0kDparXc+8OXMt8WtgAfNbC7hhNQT3f2VuOwYwq+bZhBO8sw9GmHh+jI71llc2ZnMMLPKOU//QTjh9FnCh/11hJGXRl0B/JjwTXdLwjexStueMbMvx9vDCJ2Op+qUczgwOR4WOp54+MTdXyB0fO4GXiT8+jIXM9sxbv+2XEr4cJgG9KH+yc7EjuR+hFGVdwjfYL9NeP/1IPyqbiphm+zE0pO2C8sFd7+B8DP7q+I2e5rwbb1iOOGXXPXMY+lhn0nxfsWPCSOzrxJ+QfdLDz+uIHYKv0j48cFMwkn6h6Ss26STCSMRcwgjhHkvEZB1KaEDP5nwxaOZsqp9ibBN3iVsi0vaiP0NcICZzTSz7I8vLgY+RRuHDQm/Mr0DeIGwreeTeAjd3RcQRqeOJLyGBxMuwVCXhR+A7ELt59NMvjSqked/JeFXhH/1pZemgPCjlomEHwE9Rfg1Zt2Lt8ZO4pXAy/HQ/dCq5a8RToP4FiEPHqfG+YH1ynH3+wnnhT0az0fMOpRwuoe0kOY67ATM7ArgGnevOey/vDCzcYSfxf8gIfYwYBN3/17pDcvJzMYTfprd7pWdu4LYqb3W3fOc5yadhJn9K2EkY5S7L251e6TjmdlfgSuy+yQz+zxwuLsfVH9N6QitvpiaAO7+pVa3obNx9+54hfqWiIef1MnqgsxsBeBE4EJ1sronM9sK2IIwar6Eu99Cjl+xSnl06FBEpAuycJ2vWYRD+ee0uDnSAmZ2MeE0iZPcfU578dIaOnQoIiIiUhKNaImIiIiURB0tERERkZKooyUiIiJSEnW0REREREqijpaIiIhISdTREhERESmJOloiIiIiJVFHS0RERKQk6miJiIiIlEQdLREREZGSqKMlIiIiUhJ1tERERERKoo6WiIiISEl6taLSwYMH+6hRo1pRtXSQRx55ZLq7r15mHYMHD/ZRI0akBZulF/zhh+mx8+enxy5YUE65ALNmJYcuWrgwOTY9EvK0eF5i3CzgA/ccL14+Zqs5DEuKHTGid3K5eV7qtRe9mhzr06cnx+Z57QA+yBG7OEfsijli83zzT93EbwOzS8whgB49BnuvXqOSYtdeO73cnj3TY99/Pz122rSP0oOZmyMW8ryKZv2SY/v0SW/BvNQdDOkfDe6TcZ/eUB61pKM1atQoJk6c2IqqpYOYWfqnR4NGjRjBxAkT0oJ75Uj1yZPTYydNan25ALfemhw6fcqU5NhpOZqQp8VPJ8adn6PMxgwDbk+KPPXUtA4Z5Hupfzr92OTYBRdckByb57UDeDxHbJ5O2To5YnN8lvJ6YtyJOcpsVK9eoxg8OO0z7Sc/SS+3b9/02DwfqWedNTM9mL/niIU8r2KvXrsmx260UXoLHnssvSPZq1dab3bRotHpDaiiQ4ciIiIiJVFHS0RERKQk6miJiIiIlEQdLREREZGSqKMlIiIiUhJ1tERERERKoo6WiIiISEkK62iZWU8ze8zM0i/oI1JFeSTNUg5JEZRHUpQiR7ROBJ4rsDzpnpRH0izlkBRBeSSFKOTK8GY2DPgccCbwzSLKlO4ndx716JE8L8OCRenfKXoPGZIcm3uqnFTTcl7Xe5VVkkPzXH07zw4ifYIaGJAYl2MGEqCRfdFiIG3ukuOOuy65HQceeEByLEPSX5Fbr/fk2P33yTEPEDAiz2XIc8wN83aO98gl6S1InoInz1XsK/LmUf/+sM8+DVTUjhwTPnDRRW/nKPnh5MgDD/xcjnLhml+9lh48/dHk0J/+aYvk2MceyzHN2MJnEyNzzOtTpagRrXOAU8g3BZZINeWRNEs5JEVQHklhmu5omdk+wNvu/kg7ccea2UQzm/jOO+80W60sZ5RH0qxGcgjyzPkm3UEjeTR/vvZFUl8RI1rbA/ua2WTgKmAXM7usOsjdz3f30e4+evXVVy+gWlnOKI+kWblzCAZ2dBul88udR336aF8k9TXd0XL377n7MHcfBRwC/NXdD2u6ZdKtKI+kWcohKYLySIqm62iJiIiIlKSQXx1WuPt4YHyRZUr3ozySZimHpAjKIymCRrRERERESqKOloiIiEhJ1NESERERKYk6WiIiIiIlKfRkeJGOtGgRvDsr7bvC3Lnp5Y4YlmMqkjzT9fTK8XZbtCg9Nqe+EyYkxw7729+SY9/N0YbUSVnyTOvTmAXAq4mx6Rc3vWbT05Jj//7j3ybH7vvb9Fi+/e30WMg338vD6VO4LPrBD5Jjc7xNGZoY19k+5M45Jz32scdmJMearZEcu3jRXumN6GnpsQDT/y059OaT/poce+qp6dP1wA05YlPf/zmnRcvQiJaIiIhISdTREhERESmJOloiIiIiJVFHS0RERKQk6miJiIiIlEQdLREREZGSqKMlIiIiUhJ1tERERERKoo6WiIiISEnU0RIREREpSdOzE5jZcOASYAiwGDjf3X/TbLmNmGzpUwU8lKPcg9zzNybRhBxtHpWj3GFnn50evNtu6bEDBiSFbfmpT22ZXmhjeeQO8xPncnnzzfS2TJqU/v1j+vT0aS/69MkTu3FyLMCg/fZPjt3myBeSY/uecUZy7NaXXpocOzgxbqXkEhvLoXXXXZVf/GKPpPL3H5VjCpCnn04O3e6ZZ9LLzTF/y9+/8Iv0coHtt784R3TaNgM488xTk2O32DRXNKxVAAAWgUlEQVS9BakzBn1ww+j0Qmksj2bPTm/PtBwzuay//qDk2McfTy+XcePSY2+6KUfBwJ57JoeenOP1hvdzxKZvN3gqMe6jHGUuq4hpoBYB33L3R82sH/CImd3l7s8WULZ0H8ojaZZySIqgPJJCNX3o0N3fdPdH4+05wHPA2s2WK92L8kiapRySIiiPpGiFnqNlZqOATwMPFlmudC/KI2mWckiKoDySIhTW0TKzvsD1wEnu/l6N5cea2UQzm/jOO+8UVa0sZ/Lk0YwZyiP5uDw5NHu2ckhqy5NHixcrj6S+QjpaZrYCISEvd/f/qxXj7ue7+2h3H7366qsXUa0sZ/Lm0aBByiNZVt4c6t9fOSQflzePevRQHkl9TXe0zMyAPwDPuft/N98k6Y6UR9Is5ZAUQXkkRStiRGt74HBgFzN7PP7tXUC50r0oj6RZyiEpgvJICtX05R3cfQKQfjEokRqUR9Is5ZAUQXkkRdOV4UVERERKoo6WiIiISEmKuDJ8qabmmKJmVo5y88SOydGGvBt06xyxQ/IU3Ldveuym6fMgLFiU1jf3FXqn19+ghQvTp7N45JH0cm+8MT32rrsS5wACIM+FpfNMNwGwRXLkb36zQXLsN046KTm219y5ybEjbrghKW7F5BIbM316+mwkR45P38Zz5sxLjj377PTplr45eXJy7Brb5zv69WGO2N7HHJMc+/IhWyXHHnlkehuGJO4QP2p85pRkvXvDqFFpsZttll5unpnUFi1Kj/0jRyXHHtXztvSCIdeLOG7cFcmx5523Y3Lstdemx6ZO4wY5PkSqaERLREREpCTqaImIiIiURB0tERERkZKooyUiIiJSEnW0REREREqijpaIiIhISdTREhERESmJOloiIiIiJVFHS0RERKQk6miJiIiIlKQlU/BMngxHJc4A8Mc//zm53KH77Zccu3L6dfeZnByZv+e6To7YXrvvnh48enRy6H0T0ludOsvKe+8lF9mw+fNh0qS02AceSC/3rrtm5GjFuByx43PErpIjFmCf5MixYw9Ljh09On3ame123TU5ts+ECUlxNnNmcpmNmD37Q2655aXE6NNylHxwcuQ3Hzgovdgc+601Znt6uUDviX9Njp260S7Jseuu/XZy7O67r5EcmzrlTe/yZwOjVy8YPDgtdp/0tyob93ohOfateelTax11477JsWNuuSU5FmDvHLHbnXxycuysQ9L3Ra+/nt6Ge+4ZmRjZ+IRgGtESERERKUkhHS0z29PMnjezl8zsu0WUKd2P8kiapRySIiiPpEhNd7TMrCfwv8BewMbAoWaWPh29CMojaZ5ySIqgPJKiFTGitTXwkru/7O4LgKuA9JOlRALlkTRLOSRFUB5JoYroaK0NTMncfz0+JpKH8kiapRySIiiPpFBFdLSsxmMf+7mLmR1rZhPNbOL8+e8UUK0sZ3Ln0Zw5yiNZRu4cgnc7oFnSxeTOowULtC+S+oroaL0ODM/cHwZMrQ5y9/PdfbS7j+7TZ/UCqpXlTO486tdPeSTLyJ1DsFqHNU66jNx51Lu39kVSXxEdrYeB9c3sE2bWGzgEuLmAcqV7UR5Js5RDUgTlkRSq6QuWuvsiMzsB+DPQE/ijuz/TdMukW1EeSbOUQ1IE5ZEUrZArw7v77cDtRZQl3ZfySJqlHJIiKI+kSLoyvIiIiEhJWjLX4YwZMG5cWuxZZ+2RXO6a112XHLvBI4+kx06fnhyba5IlgC23TI/dccf02E03TQ7dPH36NObNS4tbeeX0Mhv14Ydh3swUqXHBWzli388ROzBH7Fo5YgE+kRw5YEB6qXlSn5VWSo8dMiQtLnVyzYa9BZyTGJvntT4vOXLf+elzyQ0Yld6CS/u/mh4MPP98+vyF0yfnKXlQcmTfvumlJk6XWX4KAT17prc9dU5EAG69NTl04Le+lRx7Z44mJL5Tl8h1ZdfUCSuB6U+nF5trv0XqZJi1foyaRiNaIiIiIiVRR0tERESkJOpoiYiIiJREHS0RERGRkqijJSIiIlISdbRERERESqKOloiIiEhJ1NESERERKYk6WiIiIiIlUUdLREREpCQtmYIH3sf9gaTI731v2+RSd9jhc8mxi4akxw7ZPDk0z4wCuU2blh47ZFJ6bI7Zelh10kNJcSssyDNdSWMWLUrfJvmm4RiWI/aAHLEzSmoD7Lhj+hQ8hxySXu7o0TkacUPi/EwAc+akxX30UY4GNGJF0qcvWiVHuackR+aYOSzf9C28mSeYl14amRw7Kcf+5fOf75kcm2d6qNQpeD78ML3MZvRK/DRNjQN47YBvJseO2Dz9g2qPMWOSY9/929+SYwH6jh2bHHv+daslx/7sZ+lteOWVJ9KDSZ2qakGOMpelES0RERGRkqijJSIiIlISdbREREREStJUR8vMfmlmk8zsSTO7wcxyHGEXCZRH0izlkBRBeSRlaHZE6y5gU3f/F+AF4HvNN0m6IeWRNEs5JEVQHknhmupoufud7r4o3n2AvD+VEkF5JM1TDkkRlEdShiLP0ToKuKPeQjM71swmmtlEmFVgtbKcSc6jefPe6cBmSReSY19U/mVIpMtKzqMPP9S+SOpr94oeZnY3MKTGolPd/aYYcyqwCLi8Xjnufj5wfoj/pDfUWumyysijNdccrTzqRsrZFw1XDnUzZeTRoEHaF0l97Xa03H23tpab2RHAPsCu7q5kk5qUR9Is5ZAUQXkkHa2pK8Ob2Z7Ad4Cd3P2DYpok3Y3ySJqlHJIiKI+kDM1OwTOWMIfFXWYG8IC7H9/+au8Bf0mq4KKL0s9FHDcuPTbPVA9rrpke++lPp8dCvik1+vYtJ/bpp9NjR4/eOilufs8805U0lkdm0KdPWgV5XsNPfWrV5NiFCzdOjh04ML0Nm22WHguw117psdtskx675vsvpwc//3x67OTJ6bFpGtwXDSR9GqX0KZT+8pdBybFXXZUcypw5t6UHc3+OWLjwwvTpzg47LL3cPFOn5DFuXFpczjGphvJopZVgk03SKnj88fTGXHhheuxJJ+2SHLvL0a8nx652QJ5pxuCUV7+eHPurX6WX657WZwhuzRGban7DazbV0XL39ZpZXwSUR9I85ZAUQXkkZdCV4UVERERKoo6WiIiISEnU0RIREREpiTpaIiIiIiVRR0tERESkJOpoiYiIiJREHS0RERGRkqijJSIiIlISdbREREREStLsFDwNWgikTgFwdXKp7hslx86cmTh3CzBzZvrcKZMmjUyOBejTJ32qjvXXz1NurmYk23zztLi33iqn/qzevWFY4qxLebZHnmmU8kyhtF6Oa05vuml6LMA6oxanB0+fnh6bZ76QHLFvJ8YtSq+9IUOH9ua449Les4sWpb+399wzvQ0LF34xPZiv5Yg9KUcsnJQj/Lrr0mNT9xkA6/R6LTn2hBNGJMVdnf4R0rCPPoLZs9NiN0r/mGLb9FmRkqcAAmDT9AS9+p41chQMvzxxZo7oPNPqPJEj9qUcsWslxuXYx1bRiJaIiIhISdTREhERESmJOloiIiIiJVFHS0RERKQk6miJiIiIlKSQjpaZnWxmbmY5foMlsizlkRRBeSTNUg5JkZruaJnZcGB3IP13uSJVlEdSBOWRNEs5JEUrYkTr18ApgBdQlnRfyiMpgvJImqUckkI11dEys32BN9w9z5XERJahPJIiKI+kWcohKUO7V4Y3s7uBITUWnQp8H9gjpSIzOxY4Ntzrm9xAWT6UkUcDB6ZdGVqWH0XkUTaH+vdXDnU3ZeyLVl1VeST1mXtjo6Nm9inC9fM/iA8NA6YCW7v7tLbXHeCwY2JN6dPfQM8csXkG81bNEZv33Ml1csRunCM2/Y2/1lrp23i33dLibrttNDNmTLT24prJow03HO3nnjsxqT1z5yaFAfmm1RlSa3ddQOzKc1MnqYmmtbmplpVnY9x6a3rs2LHJoc/OmZMUdxDwtHtpedS//2jfbru0HPrnP5PCAHjxxdQpxqBPn8R5pIAdU3eb5JsGCGCffdJjN9wwfYqTV19Nn3tqxLSHkmPvm791Utyxx45m0qRy90XrrDPazzgjLY/mz08KA6BXjknyBgxIj+2bY6zj6KPTYwFeeeW2HNFTc8SOzxH7Xo7Y1Dmi/oD71HbzqJaG5zp096eAJZMgmdlkYLS755hITbo75ZEUQXkkzVIOSVl0HS0RERGRkjQ8olXN3UcVVZZ0X8ojKYLySJqlHJKiaERLREREpCTqaImIiIiURB0tERERkZKooyUiIiJSEnW0REREREqijpaIiIhISdTREhERESlJw1PwNFWp2Rzg+Q6vOMyP06qr/Laq7lbVu6G79yuzgm6YR92tXig5j7phDrWybu2LitfdXstW1t1wHhV2wdKcnnf30R1dqZlNbEW9ray7lfV2QDXdKo+6W72VukuuolvlUCvr1r6oeN3ttWxl3c3kkQ4dioiIiJREHS0RERGRkrSqo3V+N6u3lXUvz/Uuz89N9XZM3dqmqrer1KF6W1t3w/W25GR4ERERke5Ahw5FREREStIhHS0zG2Nmb5jZ4/Fv7zpxe5rZ82b2kpl9t4B6f2lmk8zsSTO7wcwG1ImbbGZPxbY1/MuC9tpvZiua2dVx+YNmNqrRuqrKHW5m95jZc2b2jJmdWCNmZzObnXkNflRQ3W1uOwv+Jz7nJ81siwbraUkOxTKX+zxqZQ7FspfrPOoOORTLXe73RbEs5RHaFyXnkbuX/geMAU5uJ6Yn8E9gHaA38ASwcZP17gH0ird/Dvy8TtxkYHCTdbXbfuBrwHnx9iHA1QVt37WALeLtfsALNereGbi1hNe2zW0H7A3cARiwLfBgV8qh7pJHrcyh7pBH3SGHWp1HHZVDyqNy82h53Bd1pkOHWwMvufvL7r4AuArYr5kC3f1Od18U7z4ADGuyjW1Jaf9+wMXx9nXArmZmzVbs7m+6+6Px9hzgOWDtZsstyH7AJR48AAwws7VKqqvwHILukUedPIegi+dRd8gh6PR51JE5BMqjhnTyHIIG8qgjO1onxGG2P5rZwBrL1wamZO6/TrEb9yhCL7QWB+40s0fM7NgGy09p/5KY+GaZDQxqsL6a4tDtp4EHayz+jJk9YWZ3mNkmBVXZ3rYr8nVtdQ5BN8ijFuQQdK88Wu5zCJb7fREoj5aJ0b6ovsKuDG9mdwNDaiw6FTgXOJ3wBE4HziYkyTJF1Fi33Z9EtlWvu98UY04FFgGX1ylme3efamZrAHeZ2SR3v6+9uqubUuOx6vY39ByTG2DWF7geOMnd36ta/Cgw0t3nxvMJbgTWL6Da9rZd8nNuVQ61V3d3yqMW5RAsB3mkHMoU3sX3RaA8QvuiQvIICuxouftuKXFmdgFwa41FrwPDM/eHAVObrdfMjgD2AXb1eIC1RhlT4/+3zewGwpBp3qRMaX8l5nUz6wX0B97NWU9NZrYCISkvd/f/q16eTVR3v93Mfmdmg929qTmjErZd8uvaqhxKqbs75FGrciiW1+XzSDkULA/7oliO8kj7oqbzqFJo6X/AWpnb/wVcVSOmF/Ay8AmWnni3SZP17gk8C6zeRswqQL/M7b8DezZQV7vtB77OsicOXlPQ9jXgEuCcNmKGsPS6aVsDr1XuN1Fvu9sO+BzLnjj4UFfKoe6SR63Koe6SR90hh1qZRx2ZQ8qjcvNoedwXNdWoHI2/FHgKeBK4uZKkwFDg9kzc3oRfGPyTMEzabL0vEY6lPh7/zquul/CLiifi3zPN1Fur/cBpwL7xdh/g2tiuh4B1Ctq+OxCGLp/MPNe9geOB42PMCfH5PUE4iXK7Auqtue2q6jXgf+M2eQoY3ZVyqLvkUatyqLvkUXfIoVbmUUfmkPJI+6K8eaQrw4uIiIiUpDNd3kFERERkuaKOloiIiEhJ1NESERERKYk6WiIiIiIlUUdLREREpCTqaHVhZjbKzL7U6nZI12BmR5rZ2Fa3QzovMxtvZqNb3Q5ZfpnZ38zs8fg31cxujI+bmf2Pmb0UpzbaIrPOEWb2Yvw7onWtb0xhV4bvzsysly+d6LMjjQK+BFzRgrqlQC3MIenClDeSV6tzxt13zLTleuCmeHcvwjQ66wPbEKY52sbMVgN+DIwmXF/rETO72d1ndmjDm6ARLZaMDE0ys4tjT/o6M1s5LvuRmT1sZk+b2flmYWby+M3vp2Z2L3CimX3ezB40s8fM7G4zWzPGjYnl3mlmk81sfzP7hZk9ZWZ/ilMNVLfnG2b2bGzLVfGxnTLfAh4zs37Az4Ad42P/1WEbTD6mK+RQ1fJxZnZe/Hb5gpntU/Imkho6W95EB5rZQzEvdoxlLTMaama3mtnO8fZcM/u5hUl47zazrWMbXzazfUvcfN1SZ8oZM/ukmT1U1bYn4+1dY/lPWZh4e8WqdfsBuxDmKQTYD7jEgweAAWa2FvBZ4C53fzd2ru4iXCG/y1BHa6kNgfPd/V+A94CvxcfHuvtW7r4psBJhjqmKAe6+k7ufDUwAtnX3TwNXAadk4tYlXLZ/P+Ay4B53/xQwLz5e7bvAp2Nbjo+PnQx83d03B3aM634X+Ju7b+7uv27y+UvzOnsOVRsF7BTXP8/M+uR9wlKIzpQ3AL3cfWvgJMJIQntWAca7+5bAHOAMYHfg3wlXEZfidYqccffngN5mtk586GDgmrgvGQccHNftBfxn1XP4d+AvvnTewrUJV72veD0+Vu/xLkMdraWmuPv98fZlhGkAAP4t9vyfIvS+N8msc3Xm9jDgzzHu21Vxd7j7QsLl+nsCf4qPP0X4sKv2JHC5mR1GmKEd4H7gv83sG4Q3jA4XdD6dPYeqXePui939RcKcZhslPEcpXmfKG4DKJL6PtBGTtaCq3HszdaasL/l1ppy5Bjgo3j441rMh8Iq7vxAfvxj416r1DgWuzNy3GmV7G493GepoLVX9wnnslf8OOCD2yi8gzO1U8X7m9m8J3yY+BRxXFfchgLsvBhb60nmPFlP7PLnPEeZS2pJwPLqXu/8MOJrwLeUBM9OHYufTqXMopb31npiUqjPlzZJ1gI8yMYtY9vMiW0d1udk6dR5wOTpTzlwNHGRmG4TV/EVqd46WMLNBhMmgb8s8/DowPHN/GDC1jce7DHW0lhphZp+Jtw8lDK1Wkm+6mfUFDmhj/f7AG/F2w7+KMLMewHB3v4cwnDsA6Gtm67r7U+7+c2AiYfRhDtCv0bqkcJ06h2qEHmhmPcxsXcJkqs83Wqc0pVPkTTsmA5vHfBlO+JCU1uk0OePu/yR0yn/I0lGzScAoM1sv3j8cuDez2oHAre4+P/PYzcB/WLAtMNvd3wT+DOxhZgPNbCCwR3ysy9C3jaWeA44ws98DLwLnuvsHZnYBYch0MvBwG+uPAa41szcIs4l/osF29AQuM7P+hG8Fv3b3WWZ2upn9GyGhnwXuIHzDWGRmTwDjdJ5Wy3X2HKqOe56w81uTMDP9/OoA6RCdJW/acj/wSmzP08CjJdQh6TpbzlwN/LJSjrvPN7OvxDp6xbacl4k/hPBjrqzbgb2Bl4APgK/Est41s9Mzz+c0d3+3yfZ2KFs6Kth9mdkoQu960xY3RbqorpZDZjaO0N7rWt2W7qyr5Y20nnKm69GhQxEREZGSaERLREREpCQa0RIREREpiTpaIiIiIiVRR0tERESkJOpoiYiIiJREHS0RERGRkqijJSIiIlKS/w9WT5bSQm0xgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x216 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "params_11x11=['sst','slp','hum','vo700']\n",
    "plt.figure(figsize=(10,3))\n",
    "for p,param in enumerate(params_11x11):\n",
    "    image=np.zeros([grid_l,grid_l])\n",
    "    for i in range(grid_l):\n",
    "         for j in range(grid_l):\n",
    "            image[i,j]=data_train[param+'_'+str(i)+'_'+str(j)][sample_id]\n",
    "    plt.subplot(1,4,p+1)\n",
    "    plt.imshow(np.array(image),extent=[-5,5,-5,5],\n",
    "               interpolation='nearest', origin='lower', cmap='seismic')\n",
    "    plt.xlabel('param '+param)\n",
    "t=plt.suptitle('Example of 11x11 degrees maps, centered in the storm location.'\n",
    "         +'\\n (surf. temp., surf. pressure, 1000hPa humidity and 700hPa vorticity)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://github.com/sophiegif/ramp_kit_storm_forecast_new/blob/master/figures_pynb/pipeline.png?raw=true\" width=\"70%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For submitting at the [RAMP site](http://ramp.studio), you will have to write two classes, saved in a specific file:   \n",
    "\n",
    "* a class `FeatureExtractor` in a `feature_extractor.py` file.\n",
    "* a class `Regressor` in a `regressor.py` file.\n",
    "\n",
    "You can look at the simple examples provided in /submissions:\n",
    "- starting_kit : using only the track data\n",
    "- starting_kit_pressure_map : using both track data and image data\n",
    "\n",
    "### Using data from previous time steps\n",
    "Of course, you can use the data from previous time steps, e.g., for the prediction of the intensity of storm S at t=3 you can use data from S at t=\\[0:2\\]. However, it is completely forbidden (and we check it!) to use future data like S at t=4,.. This is illustrated in the figure below, where the estimation of the 24h-forecast of the time instant 2 (red line) can use blue but not red features.\n",
    "\n",
    "- `illegal_lookahead`: this simple submission illustrates the error you will have if you are illegally looking ahead time of the same storm.\n",
    "- `legal_lookbefore` : this simple submission illustrates how to use information from previous time steps of the same storm.\n",
    "\n",
    "<img src=\"https://github.com/sophiegif/ramp_kit_storm_forecast_new/blob/master/figures_pynb/illegal_lookahead.png?raw=true\" width=\"70%\">\n",
    "<div style=\"text-align: center\">Data from previous steps are allowed, but data from future steps are forbidden.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluation\n",
    "The framework is evaluated with a cross-validation approach. The metric used is the RMSE (root mean square error) in knots across all storm time instants. We also made visible three other metrics: `mae` is the mean absolute error, in knots. `mae_hurr` is the MAE using only time instants corresponding to hurricanes (windspeed>64 knots), while `rel_mae_hurr` is the relative RMSE on hurricanes. These metrics are interesting because the current forecasting practice is to exclude all other stages of development (e.g., extratropical, tropical wave...), see [this page](https://www.nhc.noaa.gov/verification/verify5.shtml?)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing the submission\n",
    "You can test locally our pipeline using `ramp_test_submission` command line (`-h` will give you all infos). For that, open a terminal in your `storm_forecast/` folder and type on a terminal `ramp_test_submission --submission starting_kit`. You can then copy a submission example in `submissions/<YOUR_SUBMISSION_NAME>/`and modify its codes as you want. Finally, test it on your computer with `ramp_test_submission --submission <YOUR_SUBMISSION_NAME>`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get to see the train and test scores, and no errors, then you can submit your model to the ramp.studio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some warnings when building the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-danger\">\n",
    "\n",
    " <ul>\n",
    "  <li>If you want to use the features from previous time steps in your learning (for example using LSTMs), you will have to use the 'stormid' and the 'instant_t' columns. Moreover, you will have to handle separetly the first time steps, which are not provided with past data.</li>\n",
    "  <li>The intensity value to predict is the max windspeed. However, this value was measured empirically with a precision of ~5knots. </li>\n",
    "</ul> \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Submitting to the online challenge: ramp.studio\n",
    "\n",
    "Once you have found a good model, you can submit them to [ramp.studio](http://www.ramp.studio) to enter the online challenge. First, if it is your first time using the RAMP platform, [sign up](http://www.ramp.studio/sign_up), otherwise [log in](http://www.ramp.studio/login). Then sign up to the event [storm_forecast_CI2018](http://www.ramp.studio/events/storm_forecast_CI2018). Sign up for the event. Both signups are controled by RAMP administrators, so there **can be a delay between asking for signup and being able to submit**.\n",
    "\n",
    "Once your signup request is accepted, you can go to your [sandbox](http://www.ramp.studio/events/storm_forecast_CI2018/sandbox) and copy-paste (or upload) [`feature_extractor.py`](/edit/submissions/starting_kit/feature_extractor.py) and [`classifier.py`](/edit/submissions/starting_kit/classifier.py). Save it, rename it, then submit it. The submission is trained and tested on our backend in the similar way as `ramp_test_submission` does it locally. While your submission is waiting in the queue and being trained, you can find it in the \"New submissions (pending training)\" table in [my submissions](http://www.ramp.studio/events/storm_forecast_CI2018/my_submissions). Once it is trained your submission shows up on the [public leaderboard](http://www.ramp.studio/events/storm_forecast_CI2018/leaderboard). \n",
    "If there is an error (despite having tested your submission locally with `ramp_test_submission`), it will show up in the \"Failed submissions\" table in [my submissions](http://www.ramp.studio/events/storm_forecast_CI2018/my_submissions). You can click on the error to see part of the trace.\n",
    "\n",
    "After submission, do not forget to give credit to the previous submissions you reused or integrated into your submission.\n",
    "\n",
    "The data set we use on the backend is usually different from what you find in the starting kit, so the score may be different.\n",
    "\n",
    "The official score in this RAMP (the first score column after \"historical contributivity\" on the [leaderboard](http://www.ramp.studio/events/storm_forecast_CI2018/leaderboard)) is the RMSE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More information\n",
    "\n",
    "You can find more information in the [README](https://github.com/paris-saclay-cds/ramp-workflow/blob/master/README.md) of the [ramp-workflow library](https://github.com/paris-saclay-cds/ramp-workflow)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contact\n",
    "\n",
    "Don't hesitate to [contact us](mailto:admin@ramp.studio?subject=Storm forecast CI2018 ramp)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.6.4"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}