Chapter9_Fire_and_Teleconnections.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-25T19:53:07.854160Z",
     "start_time": "2018-11-25T19:53:05.289059Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data already available\n"
     ]
    }
   ],
   "source": [
    "# All imports go here. Run me first!\n",
    "import datetime\n",
    "from pathlib import Path  # Checks for files and so on\n",
    "import numpy as np  # Numpy for arrays and so on\n",
    "import pandas as pd\n",
    "import sys\n",
    "import matplotlib.pyplot as plt  # Matplotlib for plotting\n",
    "# Ensure the plots are shown in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "import gdal\n",
    "import osr\n",
    "import numpy as np\n",
    "\n",
    "from geog0111.geog_data import procure_dataset\n",
    "\n",
    "if not Path(\"data/mod14_data\").exists():\n",
    "    _ = procure_dataset(\"mod14_data\", destination_folder=\"data/mod14_data\")\n",
    "else:\n",
    "    print(\"Data already available\")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Group project: Fire and teleconnections\n",
    "\n",
    "There is much public and scientific interest in monitoring and predicting the activity of wildfires and such topics are [often in the media](http://www.bbc.co.uk/news/science-environment-15691060), or [here](https://en.wikipedia.org/wiki/Camp_Fire_(2018)) for a more recent event.\n",
    "\n",
    "Part of this interest stems from the role fire plays in issues such as land cover change, deforestation and forest degradation and [Carbon emissions](https://atmosphere.copernicus.eu/global-fire-emissions) from the land surface to the atmosphere, but also of concern are human health impacts, effects on [soil, erosion, etc](https://www.sciencedirect.com/science/article/pii/S001282521100002X). The impacts of fire should not however be considered as wholy negative, as it plays a [significant role in natural ecosystem processes](http://www.pacificbio.org/initiatives/fire/fire_ecology.html).\n",
    "\n",
    "For many regions of the Earth, there are large inter-annual variations in the timing, frequency and severity of wildfires. Whilst anthropogenic activity accounts for a [large proportion](http://www.pnas.org/content/early/2017/02/21/1617394114)  of fires started, this is not in itself [a new phenomenon](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3263421/), and fire has been and is used by humans to manage their environment.\n",
    "\n",
    "Fires spread where: (i) there is an ignition source (lightning or man, mostly); (ii) sufficient combustible fuel to maintain the fire. The latter is strongly dependent on fuel loads and moisture content, as well as meteorological conditions. Generally then, when conditions are drier (and there is sufficient fuel and appropriate weather conditions), we would expect fire spread to increase. If the number of ignitions remained approximately constant, this would mean more fires. [Many models of fire activity](http://www.nasa.gov/images/content/492949main_Figure-2-Wildfires.jpg) predict increases in fire frequency in the coming decades, although there may well be [different behaviours in different parts of the world](http://science.sciencemag.org/content/334/6057/787.full).\n",
    "\n",
    "Here's Richard Feynman talking about fire and the carbon cycle.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkzODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQEDEQH/xAAbAAACAwEBAQAAAAAAAAAAAAACAwEEBQAGB//EADsQAAICAQIDBQUGBQQCAwAAAAABAgMRBCEFEjETIjJBUQYUM2FxI0KBkZLBFVJTVHIkQ4KhRGI0orH/xAAYAQEBAQEBAAAAAAAAAAAAAAAAAQIDBP/EABsRAQEBAQEBAQEAAAAAAAAAAAABEQISITFB/9oADAMBAAIRAxEAPwDwC6mpCiDS7kfyMpGxp7XypTjhPpIlBy09fJlQj+RXdMM+CP5F5bporzjhkZDVTW5b1x/IsV6epv4cPyBqXmNjsmBE40xXwq/0oDT01Ssf2UH+BE3lljRx3bARZp6lb8OC/ACdFae9cfyLWpWLEzrod1MGqsaK2vhx/Is6eipxa7KH6UdpVmbi/MbHuX8pTQy0tX9KH6R1GmoccOmvb/1QzBNe0io73Wj+hX+hBLS0f0K/0IagkQ0n3XT/ANCv9CJ91o/oV/oQ4kBPulH9Cv8AQiVpNP8A0K/0occgaUtJp/7er9CJ900/9vV+hDkcDSvdNP8A29X6EctHp/7er9CHnA0r3PT/ANvV+hHe56b+3p/Qh5wNI90039vV+hHe6ab+3q/QhxICPdNN/b1foRPumm/t6v0IecUI9z039vV+hE+5ab+3q/QhxxE0n3LTf29X6Ed7npv7er9CLBxTSPc9N/bU/oR3uWl/t6v0IsHA0MVyrEdl6Inf1ZxINRv6s4nB2Aagjf1YRwNDv6s7f1YRANRv6sjf1YRwTQb+rO39WFg7AXQHYCOBocEknACQHggqaH8TvxCOC6DDBaG4BwQKaPOe0kcaip+sWenccnm/adYuo/xf7BYwzY4Y/eNLKt7uJjrqbOgq7J9pF4a6krae9XJxl+DIkWdRy2RTxuiqzLJlTwWJ7VlWEuV7dRnM5PfzKBwMqlKO6OewUEuXIEWz5mslnHNBFayIUL5YUQhlS5dQhl0cWqSFPPOpehabU10A6O+B0Y7C4LA1ASgkCSgyI5HHIKIlEEookkg4qJRKIJRFEQcSBBJBIRxJxwEnHHAcSQcBJxxTv4jTRLlk9wi4cjJnxypPuoOrjNM/JpkXGocV6dXVbjEln6llF0cQFgjARBxxwEHHHFHHHHAQzjiGBJxxwEHHYOwBxxxwVxxJwAYPNe1Px6P8X+x6c8z7VfHo/wAX+xGowV1R6GuCxttk89HxI9JCOyFbA1vgTOOGWLItMXYtjLJUd5Da4/aJCo+ItUR+8AFqxNi8tB3vNjwCoNoqIyFFpMFxwyYoB3a7YG12S22K+yHV2+WALcJbDY9CvGz5DYSCGkoglICTkcciiSUcSgOQQKCCOJRxyCpOOOIjkSQiSjiSDgJRJxwECr740xzIKyxRRgcS1cpS5c7EDtXxKVkWq+6jHtzJ5byzpXOTxElVzl5BvFZrDDrzkse6y9TlppRM6uCqnJPOcM0KOI3V4zJtFCNbRE5cqLpj0+i4jXe+VySl8zQ5TwcNRySym0eg4TxlWyVNzxLyfqVm8tlxwQ0NAkis4A4k4Igg44DmcccUQcccBxBJBBxJBwHHHHIDmea9q/jaf/F/selPNe1fxtP/AIv9g3GDDxo9Qo4ijy8PGj1vL3V9BWy5w5olZrYvY2KVndbIwr+Y6u1xWPUUw4IgJ4yMgKYUCgpxyzlAYtyeiAGNTbHwpWwtNsNRl6gWFFYDihMJNeIdHoEESiCUBJyOORRISBJQBHI5HIIIkgkqpIOOMo5BEHAcTggVPU1V+KaKHATmorLKNnE4/ciU9TrJ3RxnCGLiNbxHvuKexmTm75OMd8gXScp4W7bNHRaTso5e8mZtakIo0vKt1llqNJaVWAlE53p2nKp2e3Q50lvkO5DnrXlnzraEShlbo051la2vBqdF5ZV1DW8egumXLLOcM0JxKOohyScl5nSVzsej0fGsUxjJZkvMsPi2ekUeY0tmMFztNjcc7Gu+K2fyxFvis/RGTzs7nybTGr/FbPkD/FbfSJmHbjDGn/Frf5Yk/wAXn6IzOY7JMMan8Vn/ACxJ/is/5UZWDssuJjV/isvOCJXFs/dRkp+pIxcbUOJ1vxLBYjrKJLxnm+bc52Pos5GGPVKytrKksFW/iFNHXL+hj1NpdWLtnl4aJiY3IcSomurRi+1Tzdp2v5X+wpFTik5SdSk28J4GNRRh40eviu6voeQh419T2UF3UStBcdinqF6F+XQrWruMjDPY9R7gl9RyewASW4UDsZkNUSjohy8JEUHJfZkHQjmIyOY9QNM9sFhxygjscwUQY5QaAI5HHIAjkcciiUEDlLrsBLUVR++gGokpz4jVFbblefFs55IYYGnOyNUHKTwkY13FbZ2Pk7sfIr3ai23xzbEeYwxYevv/AJ2THXX58bK2DvoXFaMOKXR67lqvibl1gYyLNXQuC7qNbOzaPdRSksvIx7oU5fMuI5CrZbsYUtRJ9ulkzWh6aGdUsm5CK6mVpfjIvzvcFiKyznW4s4IZnrWTi+8mi/p7671s9zjXaJOyN7NM51rGW9jDSu3kXNJoOdtSbXOhTshLozQrWwwUdTHMGaFmcFK/ws3GOop6Xq0WVIqQ7ksj4y3OuuNWOqOwAm2FzHSMpCzkEkqJ2Ix6ELYlSAnOEcnk7qQ9gJxuSpIWmznuwDk1jJOnWczaIhVnqP5eVfQCW8R2EzeSZTeNheH5kVyKfEetf4l1IpcR8Vf4kVUr+JH6nsq/CjxtfxI/U9nX0X0JVqZCZ+FjmKn0ZGGZNbjal3QLFhsdQu6AdUN8jHEOCwicAAohcuxyW4xLYIr092TLa6CJRxLYdECcBI45AScC5KKyyldqeZlwW7NRXW92U79fPGKsIqWPneQGXFMlfbPxyYtt+p2VghIojd9TtjuZA43zkDm9zsojmSBbyVUyOT9QcnSw4hD6WnMtLAiivEEx2eVbgC5tC8JvOCbHnoQtluETK1JFDUP7dMZY9xFnekseRitxocObstaXkjXhTFRWcGXwGD7ezPXBqW0Smn3mvocq6QFkKZLDSYqFMKp80Hj5CXpOzk5dpJs6EZc2+TnXSNNTWOpU1znOvEZYC3wKty4tIjTPhpLbOsh/uM4bqTYrn1EHiKLNV12FzQY1nCm5R2kVdTLEDStrjY+bzMviK5a0kzfLNUnJBVTK7Z0ZNM6451owmMxkowsyW68+ZqM4PDJ5sMk7BpEkeQMs+RKe25Udvkk7IDjv1Al9BlFbl3mgIRy8GjVXyxWwAxjt0FTe7RaefQRbDMsoKrtAsY1gBmVQUeI+KH4l4o8R8Vf4gVK/iR+p7Kvwx+h42v4kfqexr6RCmMVPoNYuRhhQ1EcSG0eAHUobSu4ihq6HMk5AcgyEtiUVHYJRxJBKInJQi5PyJK98+Z48kWKTO+U/oIkNk0V5J5NIiQIYMuhVCDzErPmTgAHFZBcsPBOXkh7sKGSTB6MlvcjqET1QL8SS8ych6Rc1+X5EVdhmKSl5EWd4J+onm3KOUt8ESjnzJxgjmQQixblfpPJcshzCKq4zUnLyMdNRq8Bg5Oyfka8o91mbwxdhTj+Y0oy5kcLXaRVsrS3fQXCtNh6xSmsQeCpRO2u1uWXF7GHTFqytxjkDlTRErZSi4oivnxieMmW5AurfIcIsdFZQXKkQ8kyhhZMTi3VfU3b3iJjcRr560/NG+L9c+ox8ENDuRrd9Bc0emVwsTW+8aNa7pQ08OexeiNFPBUqeXBHMwubYhmmXZIayClgnm3wUS9kRnJPVHKO5UPoh5mpVHYp015waENkAi2a5msdCs+hduimmypIyES3FtD5rbYUyVotoo8R6w/EvsocR8UPowKlfxI/U9lX0R42v4kfqeyh0QKNi2MYtmWVTVNbE6afdx6CtT4jtPzPp0AuJ5GJAQWwaCJOOOCpJQKCQQF1irrbZW8UfqN1a5oqIp91GoFOHKL5kMc+YW4rJoC1kHyDeEBJZQVAO4WMI7bBFC+grfIeGCwAYD6hPOQfkEc3sWdJX3ef1KnLKUlFbt+RvaTRTccSjsS1qRV3ZNWnndZyxX4s2a9FBLvRTG11QhLuRSZm9NTlUfBkqm+0TaXozClNRslBreLPXXt16eUvRHinZ2ltk/Vidat5N7XK2G6LTu2E/JiI4xknR8QjRqn2m1cupazGzVHkSj6D1NxiV6LY3x54PMfUfy5Rw6dolSUkdiPqhFtcn0k4leGjkrOZzlJfUw6Sr+F8iHDDEKlx8PMvxOqepgn2mJ/Ulah62OzuTF5W6wQ9jnWir/CULY9p3fUu6l9wRp6+eTl5I1HOqk6FGiXMt0ZDg2bPGLo1VdnF96RlUvn2PRz+OPQ9PHl29SyKhHDCcu9g6xyqZZ8jk8LcnJGMs0gk8k4B6IlMo5toOlc0voCxug7/OVGlp1sWOy2zkTQlFIt57pBWnns8FZrYt2Y/Irz6kUiSFsdIUzKlNGfxHrX+JpNGdxLxQ+jKKdfxI/U9jDwr6Hjq/iR+p7Ct91EpTPIBh+QtkZUNS++x+lX2aK+offH6SXcwBaRxCCQRxJBIRyJRCCKqvqJ8tkUA8SWGRZ3tQ8+RE3hbGoBlFJCebfA3mb6gNFANZBewTluQ1kNBe6A5cBvZEdUQA3kTy4Y2a5VkXnIEZR1dTssUY9WQ0aXBXHtZSaW3TYlIt6HhsK5KU3lr5F66yNaWMti7tSox7r3fogYVN4nPqceunfnkcZSa3DVvZgSlyoUk5Rcmc9dJAa/WuUJVeTMb3aUY5S2NOjTKWoUp9F6su6rs7KZxWM4N81nqPMyexSt8Radbjzee5Gl0lmq1irgss6uLR4Jb9lKl9Vua8ehn08Nt0F/aWcri/QtqxHLp05P5VJEdhjzYMLFkcpZRyrrC+z36snswzmzLpIXJYAYUpZFzfKskxKr6h+RlW8Rt085V1Ncr6lvW3Yg5IwpZcm2duOXLquttlbNym22zqJct0WDgbTFZO7jatRkmTjcGKOTeTTAsbk53OOxuVBIlrYB5DRQi6TjBl7h0eXTKXm2ZuslvhGrokvdK0BowjlIfGOI/URU1sPz3QhVsdiu1hFi3yES6mapT6C30GyFyIpTM7iiw6vozTZncV/wBn6MooV/Ej9T18PCjyEPHH6nqI3d2JKVbFyYmeqwtkVp3Sl5mWUX7TGaeyMI7leTbZBRf95gT71AoEgaEb4SfUanlGVkt6S7OYsC4jsg8+wLYCbm/eNujJIn8TIEubOxuI6zoJTbY7fzBayigMC5PDCzLPQl9NyNA6oiWyClt0O69QEvdCmuUfb3V0E5yUA3kv8OjiqxooyNDh/e07h5ma1EQlKzUxb3UTT959dilNKiOV/wBFSFjskl908/UduW1p8XyeeiLMoRjDBX08owgkmHZb3TDorXTjF4C0qVsnOWcehSk3ZZmReorkq2o5WTUSsfXxUdVJR8LLHCM1znOOMll8OgpOU5ZH06VQWa9jrb8cf6syl22kafoZaeUadacU4szJ9yyUTlrtgHNxZK1mBc9yI1Z8hjOrEdYmH72mir2OAVX3iWNeltXr1Bttc1iIMYLAXKiYaoa2PLRl+ZltI3dfpZ26flWz8jIeh1al8Cb/AAO/PLj1VWXUOnqybeet8tsHF+jWDqUnlnSRhYi9g8AQwkEnhlR2Dk9yckYRUGjpZB3Cb7jAz75Zkbuijy6av5nn59Weh0mHpa/oBcisjorEGJqXQdKS7PqEBb3opCZLYNgPJmqW+gDGMWwpckUeMrHYfRmgUON/+P8A4v8AYDMj4kegT2R5+PiRvroiVUTYsKQtsyyI4HJHMAZDlgU7AXLIDXMtaHfMvQz8l/h/gYF7yIJOKhF6fPCS6Ep7hW+ATGOcPJuIKxPGwCTXUY3gGW/Q0BeBcs5D5WiNiKEF/Il5OIoXuInsPmJZQtlrh1ihrIJ+ZWl12NrhfCcqN1zxJ9BYsoNWueTiVq6VXlsu24q1s6pvp0ZVvlzPETh1HXmh985dsodp9R28uTJUdTyN08WprBydNaENPzPZZwWoQ5I79Bug02Ksye7QvUvtKJ117P1NSJaVNqyOEwYWOjEZdH5lDT3zjZ2NixNf9lza3Cl5Geujnla5kyjfS02/JlhJpjFHmWDnK62fGZCGXuWo1po7UR7OQMb0tjeuVgZVrJHIkFzczDhU5vEVll+qR0eDS0fD3jtLU0vJFnScPjBqU1ll6eFE7ccOfXTA4nJV2V+W5pUULGWljBmcVj2uspr9WaepsVembT+R6OeXmvTz3tJT7xOLpr3XmjDjpra+sWeqteY/Uqyr5tuUuErBi2nuHnJrvh7s+6Ilwe3qTGtZxOSzPQ21+RXsrlV4lgYol0BszygxmTbL7NkRSl1N7RpOiuSMCXU3tBL/AEdZFX4PY6UMw6sXXvJPJYe1TAXLdpgSDz3UAyUAwGG0A/mRQso8e66f/F/sXngo8e6ab/F/sQZUfEjd8kYMfEjadiURVTOQmU8Cp3C+0yZZOcyObIvJJARIKJAk0NB8NmeaPD/hsou+RDORBUT1WGVa5dV/KWSs49nZL5moJcOZ5yc3g5WJ7HSjzGkc2sCeV5DeIIW7IhRANbgSe+QXYFMbWBM+uxDfM8kwhKc1GO7ZcRd4PpveNYubdR3PRXz7HEl4URwjhy01HM/HLqFxWPLpmaYYnEZqzV9p8iipcrHWz8ys3uefp34PducGjw+nLUupiTnjBu8GkpVpnGx21s4xp5Y9NjzOm12XOtvEoPozd1F06qZOKbSMJ103z54YjL1QMWlGN2JfeR0ozi8oRXGymfqi5XPLwzlXWJqs5sKXUsITOjHeXVETvwiSKHWYljHUptbljvW7ofCjlaaO/HFce+gafSSlHnm1CC9SzRxPQxv7Hnxj72NiLNPO6lrJ57VaW6Frj2Utvkejnhwvb3kMOOVvkVcmeV4PxuWjn7vqG3DPXzR6tWRsgpReYtbM6yONtYmorcuLU/LqM1qctFIsKvm4i5eiJ1Ff+maNIr6TTwu0sJyzktQ0kIrZB6CONNFehbwiop9ikQ61gttIXKG5BRuohJdDK1eljNYcTeaRm63MZPC2QJXmNZpVp8Yb39Sq33S7xCx2P6FGXgMtylOLNnQvOlivQxeZo1+Grmo5vmZbaFbwhympQaER3G1wxllQKzklkSyksdSTFUEugDGMWyKHG5S9oFvp/wDF/sXWVfaPxab/ABf7EVirqXJ2d0pLqWX4RQD6hIDzCRlDUwhaCRAxEpgIIgIv8P8AAzOyXtBYoppvBpGj5EeYKksdQiogVcu9FjQL/hNrqiwLUUC7Etia22tyZJG0Km+aJWksMbZLlYpyyVQOWSMEuOAeYAo7ySXVnouA8M5Zdtau95JmFw5xlxCqMt02e1oi6pv+VmozauRMH2h1MZpUweX54Nqd0a6nKT2R5DV2O6+c85WdidVmTVaxip7Ib1GQpUvI83VenmM+cm2bXAO0bceV49RENFGUvCb/AAnSx08G8JPyJJrV+LdtUZaSVeN2jw8ZS0mpnRZlNPKPedWYPtRoKp0K9Yhans0bvLPv6Tp7Iy5e8W7YRrjzZR5vT3uqcVY39TRlqnZHGW0cOuXo56WZavCwmV8uc0s9WKDqeLI56ZHMOq9VRw+uFUNt8ByphHyG1SzVFr0BsXMeriPJ3SGv5UB2KTbwm2Wox2Ikjs468VxXg961crILKkaHANfKt+5alvmXhyehdak9zL4nweFlkNRS+SyL8vMLWlSs2SZGpj9m0N00fsYt+LG5Fz8ish0seWtIs+QutYQzyIAe4M9kFkHqApxTKXEY408uXqaWCvqa+aLKPJ+6dpW/Uyro8jcHs0emdXZyaMTjVXLNWLzM1qMpmxwl/wClln1MY1+Fb6ea9GZdGjWu8WI9CvDyHrdEQqE+dP6kgYxJ4CySqhgMN9AGZVy8SKntJ49P/i/2LcF9oir7TfFo/wAX+wGGupYfQrLqWWKpbCQPmSjKDQaFoNEBIJAIlFBkpgolBDq75wfU0NPqVYseZlIZVJxsjgDZRzWURHogyxFdpJbCHOWR81jIuR0gTPdCp7BzznYB7mgrORV0sIbPEY5KU5uTAt8Ls5eJVN77n0SCXJk+ecGod/EK0l0eT6HKD7Fr5GmOnmLeKWau++uK5a4PH1K7WUVKfs9fqqvNvYuxOPVb5hXLg6PUY0dCHmcLXoi1p13kbEJ7pGXpY573ki4rVHzN8RjqtSKwed9reZ1QxJpGzVe5oy/aamV3D3OOcxOlcv6806VhNdR0JNYEaW/mWGWFHmexy6j0c06M8jYdehFdGwzl5Wc9ar0HDdap0qqTScehfTyeXqeHk1tDrFlRsf0Z246ce+WowJBrdENbnolefAxRFu+EMSA6yKgorYRavt/qWEKsX2qYBxWCWcjmUA+pyBU+dtLyDSKiAZrYPAF01BGdVmayCitup57jKxXuekdcrJc0jB4/XzRWPIK82zV4R8Oz6mZNYeDR4T8O38Dm6tOI2tdRNZYUkkVC3smR1Jby5ehBlUMBhMFmVFTvbEpe0nxaPo/2Lun+PF+hT9pPjU/R/sBiIsMrrqPfQVQEkHGQSCQKCQQaCQKJQBEoElBBIKOzT9AUWNJV2lm/RAakHmKYQK2QRpCrV1K/K1ktWbLJXck0aiFy6FezONhkoyz1K19nkmbCLLOZCPMmXUs8N0r1esjXjK8wr0XshpJRjK+cevTKPVfdYnSURopjCPRDJywjeOdrw3GIe58Y7XHdmOU0jR49poXyg7GvkZUtLZU+7vE49x05pnNljltHJVqUubvJjr5dlp5NnDy7619Kk9LFrfJFkdthPArXZoeWW3L0L8458jrzHLqkUTcH9S9HlurcZbplJwwx9E+VnTGHmuM8DnoJPUafvV+cfNCdNJSSa3PYylG7utJr0Z5jjPCp6F9vpN631j6HPrlvno2tPAxxMrTcSTSUtjSrtjYs8xw65d5ZRJbjVHunQUevUXq7eWvC6snOlbfDLpy0zzuk9maSw0VeFU8uhrT81ktzrcVserh5+0MXjchT9diO0jzcvMsvyydXIeRdm7TCbAl0CDXQC6XcwurOUsIV2nfbLEOqr5K0gnJR8xSslPbOCUor5stBOXMttgOTL3DS2AlbGPzZlXSrTi8dTzuvq7RTTW6PQdrLGehnXRUpSeOpB4bUR5bpL0Zc4Wsxt/Aji9HZX83lI7hUsdr+Bl1n41YNDMZERWSxFbBC+jO6nQ7zbJfUzWkNgMIFvcypukWb0UfaX4tP0f7Ghofj/gZ/tL8Wj6P9ijERYkV0WJEqlnI5nGQSCQKCQQSCQKDQHBIjAaiETFZaL2kXLJ/Mr1V+ZZo8bKLaZ2TkSVAT3RWklEs2bRyVpy5kbiFTmuV4M25Zlk0LVyxKNjNCvjc9P7J6fEp3yX0yYdGklc15RPZcK0/u+liujZqQtbC2gKt+G5+h3P3RGpsUKm5PZG3NmcWl4Y+ZWTzgi+7t71J9Ey5qFU64utYkc6sVeRfyoC3TxuwpLYccYyOno7QVxog1BYTLilsYGs112nw6MPHVMCnj9i2tq/FFkSt99QOblKlXE9NYvE4v5oa5wuj3JZXmaYXtN05n5jbIRti4ySaYjTyzD5IfFtFweU43wR6V+8aaLcPvJeRnUXuPRn0FRjZHlkspnhuP6aOj1ualiMn0OXXLrz0t6bU7d42dDw16pdpYsLyRjcA0E9Xapz+Gv+z29EVCKSWEiThu9FcOnHsnV0nXs0XeqM3X02ws9602O0j1j5SGaHidd7ccONi61y6o6SONpmpoU4PDw35nj9fTxLhurWpqlO6tf9HtpSUhFlcZZTSa+ZUYXD/aCrVYVi7Ofp5GpC6M84M7iOhhXh1RSbfTBe0lPZ1pS64KgpdBS3ZZuq5q3jZiaq5RWZdTWphkY+hFtirWFvJkSnjZDIQi0m+oCE7bf/WIyNcYLfcc8JCpfUlC7ZZ2RXdeSw2snbGVeU9pa+RR28zK4d45r1Rve1cfsUzz/D5ctz+hlufjYq6IsQw0VIYkizRtEqhguWUkc+rBlPE3ghPmZhUsFhNgN7GVWdBvc/kjP9pfi0fSX7Ghw/4kvoUPab4mn+j/AGKMMsyRWXVGhOvZEqqpAcotMhR3Mo4klRDUQIihsIfImCQ1YRRyrSJxvsRKQyuO2QhkVsHV8UEmr4gFtPY7JCJKiWVrEOszy7FeTfmagqW2Jtr0D4XpI6zXRhLogLksmx7PaZxrd+O8+jOkStmrh1cEsJMe4qKwvIntdgHLJ0c7UTtVabk8JGHr9XZqdTGvdVos8Wta7OuPWXUrwpXaKUt0iWkKhHm3LMeiRKqUVsDkw0PCK+tuVNWfNlmEW/IxuO2OF8IZWCKSr1J94u6SuFqcmk0YkZ4ZoaK+fLyQXVhV7V1VcqjCPfZZ01a09Kgnl+bEd2lczac2RC7PVmmWrp7uWWH0NGmcZvCMKuRe0s5RlkI1LO5Xt1ZkcV4bHW6fDW63ND3hY7wqeoXkEitwzU1V1qiC5XDZxawzYhblZPP6vTqUu3qfLYvMdode8qFu0vULr0lXLJb4Kev0FOpz92a6Tjszq580U0NTGDJb4loeqjqaV6bSRC9oNP8A7ldsH6OBruSwVraKb9p1pgY+s43VNZqhKcl0ysGXZ7Q8QUvDGH4GzrNDXXfHkikmPjw/TX0qE68/iBk6T2mvXx4xl9NjX0nF9Lq+5lwm+mTK1fsrKVmdHbjPlLyMbiWi1fB7Urd/SUSNY9m4SjLdjq5bGJwjjtOqqVV75bV0fqbVTrk9i6zhmcgyWUM7uAZZQFadbBy4jpSFvDIMP2nXNo8nlqHizHqes9oI50UjyVLxaiNxsabwIt1NFCh7MuaZbMiulFc8vmClgJ+KQLe5mq55BZOQGzKr3Dl436Gd7TfE0/0l+xpcM8EzN9pviaf6S/YDDXVGrjZGVHqjWXRCqGMU3uE6lkmHUJ9TKBVUTuyiMRIQvsUT2I1Eoor9niQ6JH3iQCJrf2oKCr+MBZCQJKZRzexV1DxF4LFi5o46FSzZGoipmUml6vB7Ph1Kp0dcMeR5jh1HvOthB9Fuz2SSWMeR15Y6Ktj6FedvZRcpdEW28MwuOah9i6q/FI0xCYXPUWO2XrsWO0SRUoqdNMYegeG31MtHK2UnhLI2FfnMXXKFfTqdO6UvoZaMnfyru52PNcWvd+r/AMTbsy4NIzP4ZzTcpybbIM6mDssS9TehXVoac9ZMXDQQj02K2slOuWJt48jS6my/n3yArXF9Sr2mehe0mmU481ssfIqNbhyjZTzzw0XVZFJJGdXZXVDljJHT1sYfMJVy67kjlsrRvcngztTrpWPbZIVDUyi9hqY1NVqpRikn1Ag+1hv1RnTvdjyx9WpUVkarX0PEXTJQseYevobsbVKKa3TPCdvKOpeXmEv+j0PBddmXY2tcv3WwY2d2FBbhRxkYkEUtXV2klL0Op7pdnBNFSSUHjIQ6E3F5QrV6erWx5b61JfMhWLpkeiNSvJcV9mHVm7Qt+vIP4FrJXr3a/MNRXth9WejsWTK4lw2UrIavSLlvg/L7yEVoRWVyy3KGt970j59O3ZD+R+RehJuMZtcsn1j6BuSwVln6PVy1LlGaxJeQ9tJ4YF2k+07ajuz/AP0nUZlUrMYfmBU4rV2uhnjyPC9JHs9TqOSme/keNm+abfqzFb5X9Lst2alDXqZdHhRoaeL6kUdnifzFZGy6sTJbmarmwWcyCK1OFr7CTMv2m+Jp/o/2NXhuVpvqzK9p/i0fR/sBhrqjW8kZK6o1vJCq6HUN+ICHUN+IygkSQiSgkSiESghf3gkR5kgEiYfFQJMPioosoIFEhATly9Spf3k8Fq1ZKdslE3Be9no4vsm/JHo657HneDSxRZL1eDZU+VYOnLn0s22xjByk9jFlWrLpWz39PkdrLZXahRi2oREWOx9C2pIK6yMdkV3c28IlaayTGxoVazLBh0dWnjLDbb2iskxxKWPI0dPXXXDpuBSlU61Hm6sHAq/XRlqJLd8vngh2OXQrI+rwjp6SFvj3Boknak2aNlcIx7pUZn8OpW4FmkaXcZdmKeQus6yNlW8hE7SxxGzuqJmtvJAzmyycgLoTkjQsncwJwB5H0XOMl3mseZVTJyUen0fE7MKNsunRl96y6Me6mzx9Oo5XhvY2uH8bdHLC1Kdf03QRcs4hqseBoR73c3umb2mt0usjmvD+XmFZoa5dERlgLWWLyY+PFJpdDRlw+GegP8Prz0KqnHic/wCXI1a+2XStFuOjrj91BqqK8kEUo+8WvOeVFuFKjHdtjUkLutrrWZSSLolAXR56pRXVore/Qk9s4HV2ZaeSDxvErp1xsqn3X5GMeo9rdKko3RPLmK6cr+j8Jq09DH0r7pp6dPqRTZ+YhsZOeHgR5gxILOIzsZVtcPX+jj82ZPtP8XT/AEf7GzoY40sMmL7UfGo+j/YDDXVG1GicktjGh419T1/Y80Fvgqs+OhtW7wc9NLPVF15isZbFtoYyQqHjqjlTvhPIfLZOaUenmaFNUat8b+pcFSvQynjOyLC4fDG7ZY7Q5zMKqS4dDOVJiZ6CXlLJf5iVIyYynpLY+WRD5qppuLNvIuyMZR3SZdMZvvUPQj3yGQtVw+NmXU+RmXqNLdp95LK9TUpjQlqIyKlzU/MpdqSrNzTONjhE3Gi6Po9jUhbz1pow+E3qOpcZ9JrBo6eTo1bqs8EnsdJWbDlXPD+YcFyrvDrbYwTKV17b9Ah8r4xRSnOVktgMuT8y1TCEI5kS1R11ciyyb9Ty19SvqNS28R2RmcQ1j5Ywg9yauGxkk8yY9ylZX3OhhuxzTcrHnyRMdRbBYU3gaYu03yha1GW5rU3trvHmVN83Nnc2tFqaL4KM5qFi9fMspY1VholxQv3aaWYTTIUp17TTNMMriUX2+xR8z0c667vEkZ9/D98wI1KzcnB3UzreGmKyRRZOyBk7IB5IyQcB2QoWuD+QsgDT02onHE6puLR6XhnHdlXqfwkjx2kuVVve8L2ZpOt4zF5T6FR72E42RzFppkTR4ijiWq0csQsePRm1pOPwtSjZ3Jf9BGy0LnLlQENXGf3kx65ZL1IjPtutn3YLHzFPQSsfNOeWajgs9AWtwM33DHRjYUygupbexUv1EY+YFPjsY28Okpb4R4Q9XxfXL3ScV1PKGa6crOj8zXqfdMfS+Zp6eM8dSKZPoJYzyFslEMglg+aIr0Wm209f0ML2n+NR9H+xvUrFMF8jB9p/i0fR/sBiQ8cfqeuVqUUsnkYeJfU9Oo5wzUKZN7NioVztns8R9R1cObr0H7Vx+QtZx0YRgvmQ5iOfmYWTGtD5mFzCsk5Io+YnIvODuYgZzAufkCpCdTnk5l5GVPyRKKnHEllMRprlY+u5ZLKMjWcJTfPRt8jIkpQk1JYaPWlXVaCrUrdYl6o1KmPOV2OElJPdHptFqKtfCuWyuh1RmXcEklmqWfkzOxdprPOEkdJWbHtrKU05PoZmpcM7NCuGcTt1dcqbXulsyZ6eT6yN6zgFNRfqBbqlFbsDUV9nHu5bM6cJSl1Ma1ht+snPaCZTkpt5aZs8P0HaJy9C/wDw+vBnVeW5JtZUXgjDPVx0qguVRWPoV9Rw2q3y5X8iaPOYJRe1PC7qsuHfiUoRfNh7MuizTrtTQu5bLHo3k2dJxGV9SlNQeOq6Hn7FyvBFNzqnnqvQ1OkvL1+aJYasSyTyx+60zK0L0erWOecLfRlmei1FTzCeUjpK54fdplYt0ZGs0Mq23FNot+830yxLO3qH/EMrvwyXDWDJNEZNW+3TTT+zw/kZk8J7GWg5JTBbRMOaySjCOWwqWDk0aODam7d4gi5D2c2zK1kGDzGjw7WqP2dr7vk/Q0IcCpi+9OTG/wAJ00V0bGmA5K7l6iJaSSfdNCqiutYisDGoxRqVMZ1MtRTJcrZu6HiWcRtjyv1yUUovyOcUkNZx6PtYcuW9ilZxTTro8mTXdbBOPNzRfkxVsY8jkliXoQX7+JOS7myMfWa2T6Ff/UTW6S+gMqZtE1pS1N0pp5ZRNK7S2creMmdJNPcjUWtFjzNerwGLp3tt1Nal/Z7kHSxgU3uMmKMkczo7zivmQ3sHRvfBfMK9FFYjFeiPP+03xaPo/wBje5tzA9pX9tR9H+wGLHxI9LGaWEeaj4kejhFSwWLV2PdRXvuzsHOXLAqx78iVFipbBNnLCQuUtzKiT3C5kKyDKQDnIjmyJ58sNdDKu7TleBnjg/mUtQ+VpljTWc0CDPjN6bUvz3NaqznimjN4jXhqSI4dqeWWJPKA2DiE8rKIeQJbFWVV2LvxTDcsA86yXURTRXUsQSWR6rUlsJyFGxw36o3qYrajRWSbezM+3T9nJuxSikehrsjOOYsiyMLYuMkmgYytHq4xhyw6FyGqy+oXudChyqCRUu0FkXmuW3oyK0I3Rl5hPBir3it9Ht8htevlV8WLSINJxRVv0NVyeY7jK9XVZ0khycX0YHndbwy2l5hmcf8AsouMk8NNM9gIv0lV6xOCY9FeXxJb7m9wfV6pVZt79flnqIt4U4v7OeY+jLELPd6+Tlwom50mNGNun1eYrl5192XUGeirx4Dy1+onPUu2Pdflg2OF8dxirV7p9JnSdOd5WrOHwaeI4Zl36C6Lb5co9VyqyGa5Jr5FG+ix53NI8vTpp3W8vRefyPQ6PTafTxSjFN/zMqcqqtljqMru3JjTUepUUdC+Ut1uillSQyqXZ9Aq7G2MtnswbEn5iGlNpjnXnzMgFXh7ATeHhrI/woWpKZAKQMubnx5ByylsDGTxuFRjAuT72Ghk28ZQGcpZGmOwsC8Rk8oOUsbArBDEOKx0Kd+iruzlYfyLknsQBgz08tNdh9PIuwn3EP1EVbW8rdFWKxBJ9QGc2wtvc59ASVUtjdJ/8iH1EjdJl6iPyEG83uYPtE83U/Rmy+ph8febafowMlPDRpV8WcEvss/iZhJVa0uJO1LZR/EKnV4l938WZUOiGx6gbL1cWvFH8ya5xnnvxX4mWh1XRmRendCPSUX+JEZRsTfaRWPmZ9vUOrwS+oFzup/Fh+YcbYJYc4v8TPkcMFy5wsXij+YOn5am32kX+JV8iPImKvamyFtbjzR/MpVURU0+2X5g+pEeqLg1Y6lUwxzxl+IuHElOfK4pfPmKNnhEVfFJiNeepi1lOK/EQ7+Z/EivxKsvAKLitimyON7I/mMlfXH70X+Jl1BT6DEWaLuyuny2RcZfdyPesS6yj+oyYfFJn0Lg1ff4prwv/kNjq4S+9Bf8kYS6obWTBq++QzhuH6g37vbHvOt/8kYj+IWK+iGCxbodNJpwujB/KR0NJyvbXRS+qES6EFwaFcOTrrIP8UFZqI1478Z59GjOBn5EwaHvkceFfqRXu10OjhF/iVGV7uowNtoptllThDPzOjwyqX/lVor+hYrNI0tBz6HprYzh/LlFqzjMWmuy/wDsY/kBPxMus4O3VKEpPl5pP5h1amMo55oR+TZRu6/iLiXTGpDXb/dx9S0tTXKOVOP6jDQ6vwIauNmvUxj/ALkPzGLXRivFF/8AIxUExqNpa+E8puK/5ALUVp92cP1GRHqSiK2lfW/92H6iXbV/Vh+pGPA6RBrdrV/Wh+pEOyrHxYfqRkLoS+jKNJ21P/ch+pHO2r+pD9RkkyIrRU68/Fj+o6Vtf9WH6jMBfUI0XOt/7kPzFyjT/Vj+ZRXU6QVbcaf6i/MBqtffj+ZW8gZAWnyfzr8zqbo02c6cHj5lVCvJkVrPif8A7Q/MzuI6jt5VvK2XkVvI6z7oR//Z\n",
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://www.youtube.com/embed/ITpDrdtGAmo\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.YouTubeVideo at 0x7f2e08e7a3c8>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import YouTubeVideo\n",
    "YouTubeVideo('ITpDrdtGAmo')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "[![](http://www.nasa.gov/images/content/492949main_Figure-2-Wildfires.jpg)](http://www.nasa.gov/images/content/492949main_Figure-2-Wildfires.jpg)  \n",
    "\n",
    "\n",
    "Satellite data has been able to provide us with increasingly useful tools for monitoring wildfire activity, particularly since 2000 with the MODIS instruments on the NASA Terra and Aqua (2002) satellites. A suite of [‘fire’ products](http://modis-fire.umd.edu/index.html) have been generated from these data that have been used in a large number of [publications](http://modis-fire.umd.edu/Publications.html) and [practical/management projects](https://earthdata.nasa.gov/data/near-real-time-data/firms).\n",
    "\n",
    "There is growing evidence of ‘teleconnection’ links between fire occurence and large scale climate patterns, such as [ENSO](https://www.ncdc.noaa.gov/teleconnections/enso/enso-tech.php).\n",
    "\n",
    "[![](http://www.esrl.noaa.gov/psd/enso/mei/ts.gif)](http://www.esrl.noaa.gov/psd/enso/mei/)\n",
    "\n",
    "The proposed mechanisms are essentially that such climatic patterns are linked to local water status and temperature and thus affect the ability of fires to spread. For some regions of the Earth, empirical models built from such considerations have quite reasonable predictive skill, meaning that fire season severity might be predicted [some months ahead of time](http://www.sciencemag.org/content/334/6057/787.full)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In this Session..\n",
    "\n",
    "In this session, you will be working in groups (of 3 or 4) to build a computer code in Python to explore links between fire activity and Sea Surface Temperature anomalies.\n",
    "\n",
    "This is a team exercise, but does not form part of your formal assessment for this course. You should be able to complete the exercise in a 3-4 hour session, if you work effectively as a team. Staff will be on hand to provide pointers.\n",
    "\n",
    "You should be able to complete the exercise using coding skills and python modules that you have previously experience of, though we will also provide some pointers to get you started.\n",
    "\n",
    "In a nutshell, the **goal** of this exercise is\n",
    "\n",
    "**Using monthly fire count data from MODIS Terra, develop and test a predictive model for the number of fires per unit area per year driven by Sea Surface Temperature anomaly data.**\n",
    "\n",
    "The datasets should be created at 5 degree resolution on a latitude/longitude grid, as climate patterns will probably show some sort of response at broader spatial scales.\n",
    "\n",
    "You should concentrate on building the model that predicts *peak fire count* in a particular year at a particular location, i.e. derive your model for annual peak fire count.\n",
    "\n",
    "\n",
    "### Datasets\n",
    "\n",
    "We suggest that the datasets you use of this analysis, following Chen at al. (2011), are:\n",
    "\n",
    "- MODIS Terra fire counts (2001-2011) (MOD14CMH). The particular dataset you will want from the file is ‘SUBDATASET_2 [360x720] CloudCorrFirePix (16-bit integer)’.\n",
    "- Climate index data from NOAA (e.g. see [this list](https://www.esrl.noaa.gov/psd/data/climateindices/list/))\n",
    "\n",
    "If you ever wish to take this study further, you can find various other useful datasets such as these.\n",
    "\n",
    "#### Fire Data\n",
    "\n",
    "The MOD14CMH [CMG data](http://modis-fire.umd.edu/files/MODIS_C6_Fire_User_Guide_A.pdf) are available from the [UMD ftp server](ftp://fire:burnt@fuoco.geog.umd.edu/modis/C5/cmg/monthly/hdf) but have also been packaged for you and can be imported using the following code (this has already been done in the first imports cell above):\n",
    "\n",
    "```python\n",
    "    from geog0111.geog_data import procure_dataset\n",
    "    \n",
    "    _ = procure_dataset(\"mod14_data\",\n",
    "                    destination_folder=\"data/mod14_data\")\n",
    "\n",
    "```\n",
    "\n",
    "The data are in HDF format, and you ought to be able to read them nto numpy arrays an operate with them. Note that there is data for MODIS TERRA and AQUA sensors, and if you want to use them together, you need to figure out the overlap period (AQUA only started providing data halfway through 2002).\n",
    "\n",
    "#### The teleconnections data\n",
    "\n",
    "Teleconnections data are available from a large number of places on the internet. You can find some sources of inspiration [here](https://www.esrl.noaa.gov/psd/data/timeseries/monthly/). The data can be processed in two different ways: either as it is, or as anomalies (where you define a baseline temporal period, calculate some average value, and look at the residual between the actual index and the historical average). It's up to you what index you may want to use, and whether you want to use anomalies or directly the index value.\n",
    "\n",
    "#### The predictive model\n",
    "\n",
    "The model is very simple: we assume that the there is a linear relationship between the teleconnection at some given lag, and the recorded number of thermal anomalies. Bear in mind that the aim is to **predict** fire counts some months in advance using a teleconnection. As pseudo-code, for a pixel location `i,j`, you'd have something like this:\n",
    "\n",
    "```\n",
    "    i, j # this is the pixel value\n",
    "    # Read in the peak fire month\n",
    "    peak_fire_month = get_peak_fire_month(i, j)\n",
    "    # Read in peak fire counts for all years for the pixel of interest\n",
    "    y = get_all_fire_counts_for_all_years(i, j)\n",
    "    # Loop over some lags\n",
    "    for lag in 0, ..., 12\n",
    "    do\n",
    "        # Get the lagged teleconnection\n",
    "        x = get_teleconnection_for_all_years(peak_fire_month - lag)\n",
    "        # Perform linear regression and store the results\n",
    "        m[lag], c[lag], r2[lag] = linear_regression(x, y)\n",
    "    done\n",
    "    best_lag = argmax(r2) # Select best lag\n",
    "    store_model(i, j, best_lag, m[best_lag], c[best_lag], r2[best_lag])\n",
    "    \n",
    "```\n",
    "\n",
    "## Splitting the tasks\n",
    "\n",
    "You may want to assign tasks to individual members of the group. A reasonable split might be\n",
    "\n",
    "* One person is responsible for the **satellite data**. This includes creating a 5 degree global resolution monthly dataset, and from it, derive for each grid cell, a peak fire month dataset, as well as a dataset with the fire counts at each peak fire month for all available years (more hints below)\n",
    "* Another person might be in charge for getting hold of the **teleconnections dataset**, and process it into a suitable array\n",
    "* Finally, some other person could be in charge of combining both fire counts and teleconnections datasets together and  developing and testing a **linear model** to predict fire counts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The satellite fire counts data\n",
    "\n",
    "* You should probably start using the TERRA data. Files are named `MOD14CMH.YYYYMM.005.01.hdf`. By now, you should be familiar with the naming convention, and figure out that `YYYYMM` is the year and month the dataset refers to.\n",
    "* Remember that this is a dataset of **fire counts**. You may want to check for weird data (e.g. can fire counts be less than zero?)\n",
    "* The satellite data need to be aggregated to a coarser resolution of 5 degrees. This means that you have to **sum** the fire counts for every 10x10 original pixels, discarding missing values and so on. There's two ways of doing this: a loop and an array-based one.\n",
    "    * **Loop based aggregation**\n",
    "     ```python\n",
    "    for i in range(int(nx/size)):\n",
    "        for j in range(int(ny/size)):\n",
    "            aggregated[i,j] = aggregate_function(original_data[i*size:((i+1)*size), \n",
    "                                                               j*size:((j+1)*size)])\n",
    "\n",
    "     ```\n",
    "     * **Array based aggregation**: If the window size fits neatly in the data, you can just reshape the array and aggregate over as\n",
    "     ```python\n",
    "     aggregated = aggregate_function(original_data.reshape(\n",
    "                                     (nx/nsize, nsize, ny/nsize, nsize)),\n",
    "                                     axis=(1,3))\n",
    "    ```\n",
    "* You ought to discard 2000 as there are only two months of data available for that year.\n",
    "* A reasonable data model for this is a numpy array of size `n_months*n_years, nn, mm`. You may also want to store the months and years as a 2D array (e.g. `n_months*n_years,2`)\n",
    "* Once you have this, you can loop over your fire array, selecting all the fire counts for each year (e.g. 12 numbers) for each pixel, and finding the location of the maximum (using e.g. `np.argmax`). You'll end up with an array of size `n_years, nn, mm`.\n",
    "* So now you need to decide which month is, on average, the peak fire month. How could you do this? The mean is problematic, as you may end up with something like e.g. 6.5. Are there other statistical metrics that might results useful (e.g. see [this](http://blog.catchpoint.com/2017/05/18/using-mean-performance-analysis/))?\n",
    "\n",
    "\n",
    "At this point, you should end up with one main array of size `nn, mm` (e.g. 36, 72), where every grid cell is populated by the peak fire number defined from the data, as well as the `n_years,nn,nn` array with the fire counts at peak fire month. Make sure it is clear what you mean by month number!!! Note that we also have data available for the AQUA platform, and you may want to use it too.\n",
    "\n",
    "If you plot them, they should look like this:\n",
    "\n",
    "<img src=\"./images/fire_counts.png\" width=400px />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The teleconnections\n",
    "\n",
    "* You can start with one teleconnection, but you may want to explore others.\n",
    "* A data structure for the telecon data that might be useful and convenient is to stack two consecutive years together. It then makes it easy to loop over different lags (e.g. if your peak fire month for a pixel is February, then examining the 12 previous months can be done by looking at positions 13 (Feb), 12 (Jan), 11 (Dec previous year) and so on. Here's how it looks like graphically\n",
    "\n",
    "![](images/telecon_sorting.png)\n",
    "\n",
    "With this in mind, you should aim to have an array with your teleconnection (or a dictionary of teleconnections!) with size `n_years*24`.\n",
    "\n",
    "You can use the [ESRL NOAA webpage](https://www.esrl.noaa.gov/psd/data/climateindices/) to plot time series of your index (... indices) of choice, and make sure you have processed the data correctly.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Developing the model\n",
    "\n",
    "The model is a simple linear model that relates the teleconnection value at some lag $l$ with respect to the peak fire month ($tc_{l}$) with the number of fire counts for a given cell, $N_{counts,\\,i,j}$\n",
    "\n",
    "$$\n",
    "m\\cdot tc_{l} + c = N_{counts,\\,i,j}\n",
    "$$\n",
    "\n",
    "* You have to split the data into a testing and training set: select a number of years to fit the model, and another one to test the model.\n",
    "* The training will produce estimates of the slope $m$ and intercept $c$ for every pixel and time lag $l$.\n",
    "* There are different ways to select the best lag, but the simplest one could be in terms of the coefficient of determination $r^2$: just choose the highest!\n",
    "* You should store the per grid cell model parameters, as well as probably the $r^2$ (why?), and the optimal lag.\n",
    "\n",
    "You can see an example of how this works on a particular grid cell in the following plot\n",
    "\n",
    "\n",
    "<img src=\"./images/fire_example.png\" width=400px />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model testing\n",
    "\n",
    "The previous plot shows that for a particular grid cell, the model is very succesful in the training phase (open orange circles), but the predictive phase is a bit of a mixed bag, with half the predictions being an important underestimate. However, this is a single cell validation. You ought to think of how to best validate the model. Scatter plots similar to the one above are useful, but maybe you have other ideas.\n",
    "\n",
    "A useful test of the robustness of the model in the light of limited data is to use bootstrapping: a fraction of the dataset (e.g. a few years) are left out for testing as you've done above, but the procedure is repeated by replacing the training and testing datasets, to give you a distribution of the predictive performance of the model. There is a practical and very readable introduction to bootstrapping [here](https://machinelearningmastery.com/a-gentle-introduction-to-the-bootstrap-method/)\n",
    "\n",
    "You may want to pack all the model development, teleconnection gathering and analysis in one top-level function that allows you to quickly repeat the fitting and testing procedure using different sets of years."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Refinements\n",
    "\n",
    "* You may also want to consider using the AQUA fire counts data, which is also available in the data folder. \n",
    "* Different teleconnections might explain different variations in weather in different regions. Fire is enhanced with high temperatures, low rainfall and high winds. You may want to check other connections and see how the models.\n",
    "* Several teleconnections might be used together. For example, [Chen et al (2011)](http://doi.org/10.1126/science.1209472) use two indices together in a multi-linear model. This can be written as\n",
    "\n",
    "$$\n",
    "N_{counts\\;i,j} = \\beta_0 + \\beta_1\\cdot tc_{1,L1} + + \\beta_2\\cdot tc_{2,L2} + \\dots + + \\beta_N\\cdot tc_{N,LN}.\n",
    "$$\n",
    "\n",
    "* Other metrics of fire activity might be more representative, such as the sum of all fire counts around the peak fire month.\n"
   ]
  },
  {
   "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.6.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}