Dust_Bowl_GDP-Pandas.ipynb

{
 "metadata": {
  "gist_id": "a5c969ca1c23ab22ec34",
  "name": "",
  "signature": "sha256:9b46ba4a09196d1c3bdc1fa333194050deb521c84a3b95335dfbf7f39d7095c2"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Exploring Climate Data: Past and Future\n",
      "##Roland Viger, Rich Signell, USGS\n",
      "First presented at the 2012 Unidata Workshop: Navigating Earth System Science Data, 9-13 July.\n",
      "\n",
      "What if you were watching [Ken Burns's \"The Dust Bowl\"](http://www.pbs.org/kenburns/dustbowl/), saw the striking image below, and wondered: \"How much precipitation there really was back in the dustbowl years?\"  How easy is it to access and manipulate climate data in a scientific analysis?  Here we'll show some powerful tools that make it easy.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.core.display import Image\n",
      "Image('http://www-tc.pbs.org/kenburns/dustbowl/media/photos/s2571-lg.jpg')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "jpeg": "/9j/4AAQSkZJRgABAQEAlgCWAAD/2wBDAAYEBAQFBAYFBQYJBgUGCQsIBgYICwwKCgsKCgwQDAwM\nDAwMEAwODxAPDgwTExQUExMcGxsbHCAgICAgICAgICD/2wBDAQcHBw0MDRgQEBgaFREVGiAgICAg\nICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICD/wAARCAFfAlgDAREA\nAhEBAxEB/8QAHAAAAgMBAQEBAAAAAAAAAAAAAwQBAgUABgcI/8QARhAAAQMDAgQEBAMHAgUDBAEF\nAQIDEQAEIRIxBRNBUQYiYXEUMoGRQqGxBxUjUsHR8GLhFjNygvEkNJIXQ6LCUwglRGPS/8QAGAEB\nAQEBAQAAAAAAAAAAAAAAAAECAwT/xAAeEQEBAQEBAQADAQEAAAAAAAAAEQESAiEDMUETMv/aAAwD\nAQACEQMRAD8A9pbKkrz0NeR6EFQK/qKAt3CFFM/SooJiZqoHB/uaCQYxQVVqxA8x2qiySdCSfpRF\nwoqATGO9FcRBgZ9RRELVsOnQVRUOK0nr2NECcJiTmqgCoIqogJI+bb9aCxmggJXtQdEAg1RwUQTn\neg4KEkR9aIIgSfT+lFX61BOI7npVFUEnHSoIkxtFUV9O1B2YzQCUc+nSqjjJ9qCqgRtQBcJHoetA\nIr3np8tBGqSCc9xVBUroLzB71AZC4zQEDih7UFdXrQWTBO+N80VYuREUC7qvOdOR3oio/wDFBYqP\nSirpUZ32qAnMOJ3oCBwxk/SgIgq2qAyTiijiIqArWD70DzenFRT7S0p646CoovPGOvpQXS8PqKBh\nFx2oLm4n19KAKnQfrvVAlctWfzqIEplBVFVQHGEdKIUXa1QMWwmoKu2uKBF1sye9VC6kLxjNAZFo\noiqHm+GEplGR1qKsnhzgMgYqAvKJ9qBG5ZWlVUBUlYFEQG3DH5UFuWsY6+tFVVq2j60FdJoF7eSF\n9wM1jWnDVOBsfm2oJuHFCNW85oBh7eegqiefIj7mgjXBwc96CSoffpQcFiKIIVQZ/wDNBU9T1mqI\n1DBoiivf6VREkiO9EQkZlXSqJ06k4qDtBg4360HJbkYqoqpozn70FeTv2qi3KI/rQFbb6VAQtKNB\nKmTBFBTl5FBBb7fWgFoqi3Jke1BQtd+2KCug9veqiqmzvQBdRnBx60Cjgg7VRXHsaAraulAXrPeo\nCIP1oLz5qCJJFBJKh/egot6E4oKSTmg6c+nU0EqUIBoq0iO1RHJUZFAUe/0ooiVwOxqAgdFAcOgH\n0oo6XyQM4qBlq5WBk79KAyLozvHvQMfEDeZqKs3cTFA624ojf3mipU6QPSoAqfJM9aqLIWSN89aC\nVO4oBpXkdTQSdNAM6BQcIiipNm2sTtQVFi1G0+lBblJTECiDJKUjFFVcfhGKIVSsqVq+wqijjevN\nBPwqYTUFyhvaigqbGaIAppsjtQVLSRHp1oMtjHMg4NZaD1KDgIqiLlfagTUo/eqiUOKTg0BS6iUw\nAD3oDBMtKWVgFO471AAL2z6VQbmDvnqetBcacebfv09KIqFpjO/TrVA1LVqB6URUKEf1oDBBMHac\nE0BENyBRReWPaoIUiKqLBHTaKC6W+m01RcWhogjdqftUB/gfTPSiqrtDO1AL4Y7kbVRyrTPv2oig\ntO+3Wgv8PiI3oKm0k+lBU2Sv9qAblme2+9UKvWhoEF2vm2x2FVFFMAbiRuKAbaUpVmc9Kot+cVBb\nXJxig7mGd4FBBWkUFVXE4TVRUbT2oq2vGPrUEhZg96C3bP8A5oODkCetBxdG9AVK0kTNRUFXagsl\nRidx3oDJKiYmagMjVHtQOMoViTQF0rEAVFcHV7H6n9KqG7bG/XvUaOpfSMAioBO3IPWgHzcmqji+\nodYoOL+pVAXmRQVU/wD+KCEuSqgIQUme9FNW6SpWO0moo6mVlOpOwNQW+F/gyfm7UC/I7VQupCul\nVAiKC6FJTQQp8RUAHHxVC67igHzJVvREOPBJ09OlFJoGkEpVkdKw0C/o1wiQOg6GqgbyF4JMk0C6\njiqB/igb0HcyiJLvrRXc3f70RdLk7nf/AAVRZK1R2jegkn+00RUnGc0EjagK0Z+brGaBpkeWO1FM\nxtG9QSlqczQdyoG9UXTgiiDI7Ex60DTQSFx9qBxLaSKirfByJ9KCF2QFANVsjTOaqKptwcbUBTYT\n9KChs9Bkj60EaWgNqIG6y2TAqgD1qgpOKDGfaAWaqM+4SSNPaqFFpOe9UUnP9qIHqM/rQck5/rQR\nqOT3oqEqJIB3oi6lGRP2oJ1CfSirymd9qCed2ON4qCjjk9qCmsdaAiH0p8sSfWf6UBUkkTtUUVAV\n83SgZTtMCetQMMKGxig0GVoAqKuXE9qAC3UJMneqifj9MA0ELvJymPrRakXAiSaIj4mZzRVC8rbr\nRBWnRufpRRw+Z+maglCtRE/SgabaA9KKZZQpSwd4OBQabduEK1bH9ayo7akoBxjrQSpQUnT9zQCV\nboDe3mqhG4aSlPbvRGc+QPeqFlP6frRAXH8SKACnzJoKc3VkdKChcnP5UFFv9aAzDOpRSE4x9qw2\nG+3ofLa/w9oiqgNwqUAdBRSpA0461UC0yMUFYiTEelBxEYoiFDUD+YoKpcOQDiqCJegFJyO1EV5l\nBYOkJMnrmqJCxBz7UQZl0TgYFRTrLnr/AJ3qKdaWNAoq4cRE9aIkqQentQUPQbf1qi+rbtRDKFT0\noNC1OtMflUXDSQY3qKqpXr7CqiBBOdu1BVStIkfeqgarwjptvQQp/WKAZaKtqInlhPWigPFGUg4q\noQctELzP1qhN3h8n0oEH7NQxE1UKrt4zHpVC62jEAUQIpgGgC4doFBKSNzQE1qVGo7YH9KK4xG/p\nQWJHp+tBXm/WglZAG0VBWcwM0HKBOfxCgO1qI9YqKO1Or77UDCV+vrUFuYU0FxeLAwdqC5vSrH+9\nFDNznfaiKC4ExOO9Udz/AFgGgMm59elQQXpA/vQQh8zk0DiXIAnrQHbVJkH6VFMIVEUDjKiY6VFa\ntghQUD0op5bc5H1qAC1aRg5oKpuDsaoo5eRigQubqTPSqjOdf1E9qISed8tVCq3u1AFbpBMjNAIX\nOTQVNzQUU/60V6izQ2nyEZPXrnrXJsjdW6A4NB1IGNUf0q4FnGp8p7dBVQspqImgqEmPWqKqQpPz\nHO/egGpGf6xRFFTp9Pw4oBRv2qiiiRjr1qojGmO36UROOXJPpHWg443NQHtzPXai4bQ50+/So0YT\ncCPagEq9OqAfaqysm8Uo/N/negZS56+3aimUK1YJwagM0ohETiqhlt9SIKaimkX/AEV7fegt8QIx\nt1FB3N1daIlCSo1QUWyT0/tRXJsxPt0oiTbqHSKAa7VzqaITftXQcZqgZYeAynFBGnMR9KBa4YnM\nbVRnrtCVY2qgLlnpG0miM19khW0RVQmto4xNANUaz2oOk/Q0EyJMYoKqcG1BCVydqCdQ6b0VcEnN\nQEmKAwUBJG1A0yQodQT3qKKdIUQKgqoUAVf4KARURMDeqLKUYmgGXCdv70Hc/vQT8THXFBJuBtqk\ndTQS3cA+tQOs3Mjeim0PpQiTUVdnilsVw4euKI9BaNoWnUkgo6GstH7V8N9cDFQOfEIIOk/SgC8U\n4P1qhN1Sh6VQm/cIQfMZmiALeQtJKTNUZrzpSfeqhNb29VCrrw9qBfmzn8qAK3lEQcCgCq480A0A\nHb1tCvMR7Deg+oI4W7qJKDtiuDqUubFwLJ39f9qtCjtoSr+lVAXrFSTlPt1pQuLcg6d/T3qgS7c9\nsDaiBLaI+vpVAVtkbUASkj+lVC7kAdzRFAFT207VUXbQIM7dKioV+X60BWSIjrRRx/5qALtx5YH3\nqoAlxZVP5VUHQ4Uz985oGUXPm3xUUZF8J3zQMjiO09hFBccQ+1Bx4gIogib1cSPlO1UMsXiwPWoN\nC3ugRmitFpxJTUUXnJ+tBLj7YTNVAC8VbYoIkfWqiFKB3oAOIbIxg0AVJATAz60Circ8wHv1qiXm\nG9ODJoM161Qrcbb0Qk9YjoN6oVVw8ZIx6UCr1hp9KoWWnSfpQKmZHQURIkbUVeZ/1UFgYPsOtQXQ\nokmPtQEJjB+1Bdp0oV6dRQNpcChqnFRU68UERMTQcRQCWB396AKsYmaAFxcJZEn7CgTHEmtWBNWI\nlfEE0VZviHr9fSoGm7rsqEjeoLq4tvBx+IUCS+JBK8n1FUPWvix1tGlKiAehqRaZb8ZvhxI0+WkK\n13fH5ASAynA83es8rVz49t1BJVO3SrynS3/GVm60Tp0r7GnK0ujitvdKASoA+9IVD158M5omAd6I\no7cgjXO9VCDt2gGR9qoCtc529KAKn0iRP96BZ59OlRO0UGe/e6USncjNUZDl0pSvU0H6a0IIzFeZ\nsJVuwo7bdqKRc4Y3rJTsKKp+7G9szVHHgrWntQKv8JZA82+9BmXPD0hStGR61pGZcWxT0FVCamlz\n29aqF3GF6qqKi3UkZFBBYWoUUNbatu2AaCqWVCVDNBxKtOaIqlGqgvGMbUFQPOJ+sVRIMb9DUHNB\naiQDVQQKVHYUHcz86CyVnegOh8/0mqGmbvzfpUDrd0NxQP29/AFRR/3i0owTmiDNvNmCFTVFviRq\nigvzUmgE68idIMKoFn79pkHWaoyLrjy9fkOO9VAE8fXMEyO8UDTfE2tKVrXE9KAL3F7cTBBPakGO\n/wAYklYX/wBs1UV/fAUOmaAT/EAU75opBd1qiPaaIgLGvEUFtE5/WiraT0oi8HtmooqG/LMRQWKM\n/wB6CNECdqCyVqTgfaiiJe6GoLBUb0C15dhCCZzVGYOILAyTn61UDXxAgzPr6VAnfX5dgTVCjTys\nwKCVPq6nbcUEt3GUkH3FA18YSN8VBCrqEyMigWfulOHue9UBS73OKAzdx5tWr6VAYXKiMGKKhNyd\n+vSgum4d07/0oCovlIgpVB9KA44i85BUSr3oHxxRJaCTg1AF+6Czign44ABB+agVdugT81AB69SE\nx96oy7h5SzP2HaqhF53fsaD9KJu1cyNRzkivM7DC8khR+tEX+Nbk/pUClzx+waKk6hqAkVYPP8Q8\nUrehKFctPWDGPetcpSDvHXVIEOlQ6elWJS37xunTIURpqwNN36VRrH/UT0qRTzds1ceVHmA9qirH\nhBnIz9qUgL/DVoFKQv8ABntSgS7ClFUcNUqY/wA/vVqRyuDPH8EilIGrhTqdxHrVpAFWmneqinJM\nUFTbyZoCsWxHTNBdyzMyMCiQBTA/tQV5JG1VElCwB65oJTIyfagOla/7xQc5eFsR3qhNXFHe9VDV\nrxhffbpQaCOLACdx096gqvxAiYmFVYE7jxDp/Fn/ADrSDFvONOvYJwCT96qEjeudaolV+4PxUEK4\nsuN9tjQLu8S26+lAt8WtStU+1ARNwep9ZoL/ABEg5mg7mGJJmMUF/iYVgz+tARF3Ch1qAnxZoDIu\nhj7UU4LlvSOgqCOejMGgEbuMhW2BQVFzEmZG1BxuEK96ChudB+bagzbu4dccz8p+1ULKSqN6IDyF\nEYJoF1sOg/zUVVThQOtELrcUZ7/5vRUoUncnOxogmogGFT1orvjDGnpQVK9Qxg7mgEQ5iBIoGbdK\n91delQMqSnvRVkD1Iigo44Co+hxRAy6OmTVBEOqj323oCh4kmTgb1ART5HXH+elBT4lSDvp+39qB\ndy5yTJ9dqoCq4HegBqkye9BCmNRPl9xQfooXdiZSVbxkDEjpXmdQ7q8sWEE650fg2qwecvfFbJCt\nCehzWuUrzFxxVS17nO/pW2QUvFQImTP/AIoDBwCE/eoOFyoSQfrQT8YdlGge4fxd23XqSZHrUit6\n38SK1jVASd6zFrYPEuHKa1LMGJ0zWY0zLnjXDENKCUy5+E9RVhWS5xbUfTuKsZaXDuIpSpKnBqb6\n+1TVettjavtBxmFJrLQPEmbJu1Ut2EwMd6I+c8SuTz9KFSs711YXsrgKBBzVDyA2aimW9KE+53qC\n/kX7UAV2oUJ77Uoj4Ax8tWor8ARuKtIn4CelEUdtiBQZd8UtpJUc9q0jz796dRj61plzd9G5yKCx\n4opOJ23FAo/xAq/F7VQqb1ZxqoORc9DQMpcnFFX5aHDlUE4NEKPskKICpAoFuUSrJgVUWTrbE1FS\nhweWcURYukRnbaip5p6bUHcxR+bbvQX1bn5T0FBUvK6D2oG2nzp2zUFjcrHX6UHfFqiJigAq7UVw\nKC5de09TNBRVzB86vtQDF4oGOk9aCXLrec0UJVwTQV+JVtGRRBmjMbzvRQnmEjzDIKsg7UCbjCkk\nnRpT3oFErjNATnR81AJcaoH5VQTWYhJJ6CoGUMPGDER3oDJZUnqD7VFcVOZmqjmiVnzHFQGcZSUS\nMEUUryyM7zVQVEwDGOx6/aoJhO3XvVHawJ2oKqWCvB+9BVwJKelABSBPcUHeke0UBAjyGKivpauM\nOSlSTpIEY6f71zjRO6vrh/K1/Nv7VUZ7xQSdvT0qgCUJ+YCR/eguFaIgeaqihUR5vX/J+1ANb2ZG\n1RVeadIk770QZtzyz02mijB8JnPtUFheuq2JIoofxqkz1mgum4UYJj1oNJPEQ2ynVvsakWmrDxK4\nwjlpXGdppCu4h4jdfSUlee04pErzN1eEOa/5jIrTK1nxEpcGd+9IrUHEdZhOKkWm2uIysdQKkGki\n4QUakmorRsg27oGxXIE1B6FrhLIZHM+tZqqr4IwRqbVntVoA7YMsBS3VgBPeqjCv+K8KQZABSPWt\nYm68fx/iNs4kclPl6Gt4w8w68VLzmtI4r3igGOYsYFUVWBq/WgArUDP5CguzOPyohsLIyfvRUhc7\n/NQDUoz29aCgX+W9EGQpJEelFQ4ARsAPQQKALjYnBn8qDhriMes+lBCV7iYCh/goKG4M4zQHt2XV\nmTjtQaCWkI61ANzQTQAUpCaCS+BHcbGgGq58ucd6BJ65SRAqgYeOkHcUE/Ed/tQSbgbJPloLtvJI\nFQPMPgJGZ96KvzUkEE71BHlUNEDSP0oMS8t1svKhPk/CaoVUokmqOQl1zCMDrQadu0hAJOT3qBnm\nwsDV/tQCW8ZNBHNKjn3oLhSImKC6HfmyB3BqC6NE74PWirHlo23NABTuTGJqoVduIn9aAYfESD7U\nE/ET1xQcHCcHANBcEDegOy8NK+xFRXsUOLKVfyj86wobhUIM9KqATmTQQtcdQT0igGpSYydWYzVA\n+ZGRvQBU4NOaCEvGPQ9aAodOdWOke1AZtQKZ3P8AeoCogJkYNRQvlVOyjkVRfmxmfc1ANVwTiY7C\nqB/Er2J+v9KqJ5yifbpQXU3rBB2NRWc7rad9zg1UE+NKds0DjHEDA81IHWL9SOvl/lqD0Fl4natm\n21FGstnUms8tUy94/vFavLpbO1OSsxXje65hUFnPSa1ylIX3iS/vPKp0kE/LViMl+4cjzHM1UKrT\nKdUmeg/WqgGhAUSfpQTz8du9BPOAQQOm4qhUznEqmghZQPmgR0oIZXqVO1EXW7E0UMPATHfFVFw/\n/wCaghJRnuetUcu4x/agp8RqNBfWneook61Z260EctOrNBYrt09PNQd8SmfUZoDKuZTvUC6nwKoA\nu5b2JoF13QEEH6UAFXkpjp3oAKcVuM+9VF0ugSftRQ3H1+wNBDLqy5FBoswB/eoDpcRq71FEIXO1\nBLerX5sUDTiWnEK7moMi44c5J0/L2iqKpHLGjaN6C3O0neguHgkyN5kUCz9ynUepoKJe1kAZoCqu\nVIBx5ulBwuCBvvQED+B1oJ+OPyk7UFFvCDnzRvQLKclNVAlKV03qKK2vGaAgdmTtQcoz+KgM2cKA\nMSKD3SQoIzNYUq4pwmAmY3+9AMqUG5GFbk0AVrknSJ/1dRVAtZO3vIxVHIbkAmJPTtUFVNiROffF\nBUoIwd/SgugJ32oDsuDTH2oCOveT1nYVAsp8GSMRRVS+nB/EKBU3EqMfMrp6VUVLythvNAVDpmFV\nQ2gkAEdPrUFXFa2zOO3/AJoM+5Bxyx1/WqLWwcSc+WPtQOKfSkZOagoboxAV9KCjl5iDVAfiZO8A\n9qIK28vRqMZ2FBYuJPXO1UV1nRHaiFnOZkjb9KBR54pOMKFUGaS5pBkSdpoIlQXp+uqgC4ka5H4t\nhQDCzE6hjYVUcp6YTE+tAZsBCPMcnpQV1DUUjM9aCyXE7KyaCmscsj7UVRtzA/SiCpJ2AoowW78t\nQXK8Se21AqteqT2oIQVA43oJU/8A70A1OgzFAo86RnpVC6tSjP2oKJcKZzRBQvUJoJUvyDH0qKDK\nFqIVVDTKUA+lQOIUNO09h/eiu1aFAtmT1TH9aBlDxVBP3qIvq32mirc7TmgkLSTHT0oALtgrzAUC\nDwKXD3oBqCo7dxQQdAx1oCIUhswN6CHF6uue1AOVdOvaglUgwqR60EFyOtAJbwIPfpVHDV70BW9M\nbVBy9MA9t6C4VGxoJaWiQIxRBVHEjBor6CjRpSlU6pORFYaLP6EFcHbr1ogBkjO0CqKKZE7Y6qoK\nraQkjMzuKC20nT7mgE7AWRtpNBTtsYoOAClE9T9qAmgpGqcn+tAJYXq/6sCKCgt9SpkjpQK3WpvK\ns9ooE1XCCrt+taQxw9RcSRH1qKcVq+n2qC7VyUpiN96Cq3upxNUUCwe496C+puY9INQBUU5TuelU\nV5RCdXSKDuUpcHr/AHoFykpUAsET0qoOtehOImoFuaoqzvVRYOlQ0/8AyoOwRv5f60AXWkmOsGqK\nh7PaKDlvah6igASVYoKLVCh7VUWaW3p1Hc0F9SCN6ASFaekzQUU7Bg0BEuT5evWKAiGwT3xmgPOm\nAMwcd6iuU6ZmcUAHbgyYqoBzx03oLfECQU79aiiakRvqPWgEtZPl/KgUeQtRgY01Qs4pST2IwaIX\n1+b9QaBtpatXpQFecAAHX0qKSWv+KkdqqHmXBqAiQainIMUEISdx96AqVkYNBYu+YDaMHtQEU7MR\nkRUBNaZnaiuU4AgzB96BF6deRHaiFndpTJPWqoKg7lUx6URUpdme9FWStYORntQGS8NO1BPOxBOB\nkUQu8pJ3xO4oqsdRigsh4HANBfmdelBTnA7UFeco7mguh0Aj9aA/PT1OKg+jspUXUtmSQSBWGkXV\no0JUrB9MiaIWASN8g/4KDllMRsaoXcHWfrQAddUk+9AsHyoqMxqqgiXf5u2KA6dODH0oC6pGnfTs\nBUBYbV79O9BVQ0ketQL3FmHdRAnV64qjKPBn+dB+WrUaVtZllGkpgdqihLn6dfpQBKVLyP1qgSl4\nMdKIltatMflQEIBb1j5qKq0Bq821Bc4iFROaC7WtSfmgbqFArf6AAQY60Czrh0SlUk1WSxXAmciq\no7KpEmoCJVFEQpwEH16daBZa/tVAw75jsKAKnQnNVEtAuLlR8m8UB1NNRjAoI8iRn/PagXLi9RoB\nFcmTQMNqKUYoHmcIJXknb3oJmUkq67VFDz0FAEgldULXKFIGOtEDaU4IoDNrgxRV1rGnVO3SgVJW\ns5/FtRCdyvufegElSZ38xqh5qEAH8NRQ3VJ7z2igSV82/wBaDVtGzAVkig0nQUpEbGoFQ507d6Cj\nj6p3oF1XSwqIE0DDVwU5J9jRRviVTIP0oi7a5KZOTUVLgBBmcUCvMTMRntVF41kJPWgMiyStHl3N\nQQ7w1UevWgzr5RZKdKfLt1qgKXdYPpmgK20SNaj7CiGoRgRgUUF/SXNYgKOD2NAF+dO+TuBQDEAD\nt0oBur05oCN3YMCTHY0BColM9aD6uV8tpRRkI3z1I/OubRV9+G9h/WiF1ujVKctxMev0qihcCkav\nmKfyoAO3KchO+9UKOqK5AI70EM2z2qYxtQE5IQ4czqoDfKBBmg4XCQoCZV+VAVtZUqR3qA4dbJAO\nw79KCVXCBGmJH0oKodUFajj0/wDNBynEKSqSYPQUAAhr/q9/8NBCkN6TPy9BVGdcAavKPY70RcNN\nwANx83f1oo48gGk9NulAKSVZoBuk6c/SiBJcIn9KBW6dUswn71Qu4hSUAKIkZx60QJSVGJM96KNz\nUtjHm7iqOS+SmTEURACiNQyelAOFqT1kZ+29Aq6VzNBLLLz2EpxVD44UttrWfKrqKgAoubQaoMza\nuPAqVgDc1BZ+2bA+aAKqM1LSyomcfnQMtoUpwA4FBpIISPbE0FJQfmOIqKC46kq8p9qqAAnUR9aA\nbijIiSKDkHyQKDjogQMDeKihKydO1USWdCMGaDKvlApJAAUN6qEkKcJEY7UGiHVQATtioKKcoL2b\naVvwrpmO9FehbCIGnAAqDnlHTHTqaBEhSuwAxIoAuNLjy0CSkuNzq771QRD6lAA/L0oG7cDUO9Qa\nSGNScR3JoqqkpT5enegDyUlwriPWoLhOkyE79aqDtuqCcioovPwQd6DO4g2FtnqelUZrTS2iUn5j\nBMUB9gD+f96CCsp69N+1APWSdUz3oIUr70QF1TgT6DpRQCpRwD9KAwSjrQFSSfKk0H1D4gZwNs9P\nc1zUncPjaIkbjFUIqUVLhOTVQPnKSCkZnEUAlPdZ9KoEl/UvH4sUGilYjTq6VFDU9qVFEc695OgI\nP1zRQbQqU6oK/Dj0ohxxfK1ZoruZIJGO9AJD38SO+xoCu3A5cDE4oKIfAB1H0k+tEDNyMjH0oIL3\n8Oem1BRCkqI1bjrQULy9ZM7HP/igIleSonHrVA3LiTiO80FJ1nPSgq6nU0SgZoFPNq83z/3oOUEk\nAfrtQVcQlLYgebvRFGEpWTrg+tARxGkYiiis25UmT5UmiJXZ5GnKZzQWcs2lp8wGNqqjtJDST/So\niVK1EJ9aon4donUd+lBQtEKxj2oBFrqcUAE26dc7dhQSpBCu0daCA4JMCgE4ptIiaAGsD2FVHL0F\nwbpT96CriwISDPegGVmdI+uaArURH60VTZZ61AG5c09d6ozHFgr0+tEVKVJnA9KotzZRPrn3qCvO\n1e5MVRos2/IIcCsxUU4i6kTP2oiirwExQU+IIGPl6iiobdB6+9AO6MtECgWtmlKA3gUG4y02lKfw\n+tQMNu5V5thRVy02v37VEUCEBMfaioc0DGaAAUArPX61QSNWd6CjoJBCun5VArc2yYCgc9qoTccc\nCSDkE5FVAOapRjvRVpITncVBRJKl0RD8hBnftRSaXO33oDBZoHbRQAlVB9FU7/DMwMbdz371hSXN\nJaODzOkYoFFLI8v+e9VHJlxwgJlR+9AC51oURseoqgTCgXdts0Dbtx5pGPWoAfEEHVORVFV3O8UB\nbe7QY/CdiaBwrkAqz7VFCUqFRMT0oirqxCcznH1qgb7kIEdN6BdF2o9T2NBKn0lzv3nvQSi4UkAR\nQEDpSiYoA6+ZnvQTqSickoO9ALnBKtiex96IeQwlH8RW/QH9KKvpkkpPt9KATzLaj/E6bf8AmgAt\nsJgURXkiIO2+MUVQWyUrlOKA4YxK6BptlwjaBvQSWV6cjHagA6uBneiBhVUdzE/93SgIh3Pr3oLF\n6JA/80AgrUfNigtpC+vtQLvhRx060A0paAOk52qhR/mE7UQBTitWTmgFzSSaC2ImZP8ALQB1/wAX\nPWgKl3CpnGaDkvzQBvFhQBG9FKtW7inZOOtENhoBRMz2ooLrQGce1EJhQDgVgCao00vGNNQXShSZ\nVtNBReZ6UEoUkoGaDtWkkiipaCnjpO0b0DSW+XjpUBNWrrtsKKuHQn3oGmXiUTsfWoiqlqkafqOl\nBK0SjUTntRSykEmaoIgqAxUFySaBa4hJI71RnLEKUNs7UFRAVmg55JUBpoOQ25JPWghaFJSSrPeg\nW+H1HWMGg7lgdJV2oLNLIEnAnag+oO5bGAoJJ0r6gb1zUq+UhKFAzI1SfeKoWbYDkzgD0oCsWaEu\n75CfpPvQDvLFtcAYOZ9BVCvwiGElW/agTuFHcb1UDnVHmyftRFHnEpnrRQrZ3zT1O1BqfGHTHpj2\noFzdEuAxNBIdlIV16UBl6SCFGT65FAoC2JSPmn3oBur8+0J3FEVDmlWetUXcuPKQs+Yf1qK60Upf\nsPxVUMLXPkV+VQUB0rwnA6VQz8SpwecelRRkOrRHT1/WgIka8HpsRQDdUk75EUCjijO1EU5+SZ9K\nA9urWYUr1FFaIKlDy1BV138utAm4AZzVCK1wop7GqiFHV/kUBEq0n9aCSs/WghLuKIuh86aKop1A\nEdVUAWzjeqBqQkyDvtFAndpSzAGfWiFA4qZn2qiVOq9zUFC4CDOKCNau/uaokaoE/L3oJVBwnrtU\nBBABHTFByyVdcg/5FBRxJ2Oe00CS0FKo+1UOJJUoGoDaXl4OO9BY2hIBmaK5No4cmgG6w4lQB+tA\nyzqjaDUEvu9BvQVD6jCJyKA4SBCp+lAylY60VVTpGaIgPz/Sijam1ZoCpQnl469ag7kqmTQIvMuK\ndP5VQu5aEpP81AoW1pcOsbUFkH6TQHSdKTG9AuteogKzNAVpCIP8woBvMFQ8ok0ABaq679qD6SWn\nGtQ2SeucVhQ9PMbBmTskenegunRkY9tukUFQpAWDjFBdbiYHllUxqG1Am+oKGd07f59aox3yRrJ6\n4xVQmt1CcfmOtAs46lWxyaoqm4gyfuKAq71RQCmSTjP9KDg86len5iPwjpFAywVFIJGB/NtUQ1rJ\ngnGregXdUElWdsiikXn1apEx96qKfEeTIzRTFo3zDKzt0oNPlhsANp361B3LyFKG31oBOFC1noNs\n0FkvNpVE46HvRBS5rII6bR6UURV0UpAOSe1AF16UxGes0Qu0sKBSZmqKSAdMe9AVlwock4nYVA+m\n4UBEH1NFcp/UN99x3oBKaVAVB3ohC4lLpnE71VU5gI8u4oKcxUz9qIut0RjJqiyVnTn/ALaIE4vy\n7xQQ2hSwFH7UBCPWIoKkkiO3WgTukagOlAm605q8oMVQus6MSZ3NQCU+IiqCoc1D+lAxolIqAkdh\nBoDHlgZGSM0C6Wypzy/KO9A7yNYztQCcs0iSraghKUgiBtvQNIUVdvpRTzISBEDzVB3k1UFHEIUS\nOlAJbaOn5UVn3SCII+XoPWqiLZpxSgTEUD5QNIE461AQKRpj86BZ5cAGgomSBmimcJkHegMlWPXp\nQFS4T12qAwjfvRCziCJMfWigyiMiaCnIb1TGKoDcWsDUjbtQADCtUkGKAsBOIzQEEDNQCCwVzQfQ\nXroKVkQ4r5ynAz+XSaw0Ct1sNHQkE6oJOFYqhVx4K80ZojtKyczq+1BRS1AjPXY9qoESopKdlD9f\n0oM95tSpIHy7miEzbFbh3TG9ULOMOIVHrFBBa1dPrVBk2iSgasHNQQlCWyoj69zVBua2tyDgTRBS\n4EpMGQevtUUsSpSQYg7fUVUAUnzeqe9Bdu1VOjbEE0GjbNts4nbaopoOIKBHSgEXCcRg9qAVwv8A\nhakjHWiBWtmpzzflQMBHLmM+nrQVS0VQpeAKKi5caqoVRoEmaAVxcgwEjPaglKbn+XPSg1mEL0AL\n37VAdJbAzGfrRV+ak4HtQI3rSFwnvnVVC4t0pTH59aBMghZMR6VUFTttjoKCo8x0b0RVbAI/SglK\nggaZ0j1oJB1ZBoIlImaKG8U6QDRAVFCSe9BjXukOSB7xQBQNSpiqGgEggCoH20NaAVKjtQVKUggJ\n2oozVu44s+UwOtA2lAiIECoJViJ+1AF13GmqjrdKdW0miiFUHG3agsh1ffbvUFlOztnvQclwn/ag\njUVR60AnkahgzQShGjFBdK/JETQKOrWkwKoEPPE7mgda0BHm+YVBRbgn0oDNL1J70URCyjFAU3A2\n6VBR58aN6oCFpKR0oOSqcA4qA4IwN6CqVJkqNBCm2lKmgWfBE+m1ELtmd96qvoj9m22UlDnzJCsY\nE7EE+lc2me+UckqnzTtVRzCmnR5iIoGBnB807E0CtyNVz82oHze/+RQGYbRpn/8ALemiFsNK1dxE\nY3oBGyRoKgnJzIoE3uGktagmVdKozzauj5UwraqgqrR1KSFjcCKBXkq0OeXciAd6CvKUFJB2Pyj+\ntAxyT2+XGBmg4W7nMSkYT29aB791NKBcPzTlXtQS5btpCf8A8aBbklSzrPtQcpaEJj70Al3CNPl3\nNBW1c5rZRvneqNO3OluMRH1qCpbQogxtn60EOIEehoM+4YWciN8UQncNlo+b61VAtl/xZNVG1bfK\nFEZqBhx8adqgUW7227daKjm6VAaoFBbKvxTVEpT+A/U0QF5gzjpQAbcUt1PZO5oDktyfzqojWkg/\nrQLKSgK365FASDqxiKKE4gldAN6dUdtqBRxS/l6d6ITWzPSaoq1IJHUGgNyiZ1UBkNyKgopzzR1x\nFBrW58gneM1FWKm07GgG66noaoXUsgfpREIJmRQWh2Z/Oirgr6xRHJd0YO9FULxny4ETNBVp4zvQ\nNNQR75qApSkq796KWe1D1I2iiBobKpO3aiobRpO2+1EXV/pPSgWdwc9KqmrZeB0FQO60nP3oF3ig\nDEz0oArSVQd/SggIxQQkrA/09O1AZt0jBoIU7NByXtPtQDdd1fN83SKDggTq+/8AeoPauOnRMkj+\nhrCs3UVYAM1oTbtO6lqUjB6+1A6gyCE9Kg4JTkq8wTPoYoKc5O20f4aCvxCvKJIH1oL/ABHy5nvV\nEl4JSROTsTQCLiU53NAYAFsTsf8AzRFLiwLidSBBiCfWlVYcNtiwgKjmAeb6ztQQqzYgICoogbQY\naO0g0FXLgFohGY26UAkGTkwo7VRV9MCPzoEHz5cEGgUMeUDJPzAVUO2zIS2EFW3aopkqhPaaICm5\nCTBOT1NBLt30GRRVEPA5O1EJcSVrb2nOM1VZaAdYHWqjaZuDyBP2qDipam6BVbpComO9ByVqwcnO\n/SgbZlIJWKKKzJUSOtEGaEg+u9Ak4wgGUKx1NBX8JyaoE4/pk7iNu1ELBeoSdqBlGIopkNleRioB\nu2/l3z1oF3GQIP4j6UCTjOJjeqgCLaFapoDISr6CqLCQmoAoJ5mBnrQaAWsdYooajk/nQU1BRiiK\nFR1aZoCtgaxO3WimVPISPTtUCzrpOEmqhJa3J60HJccM4+lFGbbcKT1negeaQpCesd6gMjUox09a\nCQkaie35UVVakTtg0RVQSvSE9KCq2SgSTJooCRq/6p2oL8r196C2oozO1BGrmKGNqBgNjr0qChal\nzBgdqCjjawkx8vaqABeZ60Bh/wAufuKAaB3GDQSvRiM0FyCpvy+8VB7WQXI0fwRHlVsVd/vWGl1t\n26GglKfN/MDI37RQA86EFO3YenrRFUN/wyrVCk1QC5WRKTj2NAi5caVRMx2qiqHVKMddxRDLbakn\nf/O1ARzVpOxP5feiqM6lrkZ07UQ8nUUjpE0Uyy5iDtEmagUvnQiCNz0FUDDxWdU+YUQK5cIQTGDu\nJgetUJ89Oo6TpBjSPSgMpJQDBkjfpQCdVrTvnY0ACw2lBCvmNAs4gI835VUTzNHrOSKKIHFLOpXy\njaoBqTzFEp6YqoEtK0DOIO9FQhwDy7E5miIKfJnNAm8mMxVF7d+QoK6VAdq5EeZUkdqCdLjx1DPS\ngdaRFuAoUBjBQNG3rRQhKXN8GiGtSQANqKC5yYIohBzzJISdNULIYcWYVkTtRDyOHJCPMZPSooiL\nQI6apoGW0wmD9RRQniNeBRHLAW3AFAgpsKkR7UCK2ylOdqorMgJnaiKlMCYxQShtPzdaC2v1oq8g\n+9AEyie5ojkfNlX1NAdKtIPX1oK805ET3NAIuH27UVUKBVM5oiIUPOKBi3dO3XrRTrbg6moC8xAb\nmM9qAAfMnpQV5oIz81BPfTiKDvMTnEUVwHUUAlvaFQrFANVwFnFBds6VY3oHg5iZqIGhwqcxRUvL\nhWKoryErzUBXGQhsAYmqFHFxPSgUQ4mfMetENodSQfSg9+sgIgJCgk+aOp3kntmuTYbnn/5kKiYG\nw+lVA5UJAGY+kUUElOP/ABVQo+vJxPcVQmQlR+XAyfaqjtYn+U5I+tBcP+fJ3GT+tARt3Wj0ONPr\nFQNWobAM7/rQHCo3MgUVLspbBSd9/tQZjieY4NZgAbVUFtSUq9IzQdxASgACURqz2oERyURt3JoL\ni4QSYyaDkupk9qBZaxzNKhRA1PoCfl/rVVmrdVz8bdaIYS6Q2Mz39KKaRAAO1EUuEFapO29FDQkS\nCrftRF3FaUgfailiUqkDM0QmzqSpSIqq0UMDlD+WetQalu02y3G/WiBlYSo5ie9Au5cAKlORtQQh\n89BA7mg5V0FeXvvVHJWs+WMnvmoLMtq1GdqCW1BEiKBnmJIz9KKkLj+5oiq3PL5T5qKSccg985NE\nVLywMdaCzKXVSo/agpepaS3CqoyuuKIPoJT/AEoCJQAnaooLiJcxMVR2lQ2PvREFLqht7UAnUZ9e\ntByFJGJoJW//ACigEXCpenee1BZEJP8AgoLLUY0zg0EMBQOTigbCjO9FEKlJxMVBRemJnaqK26kr\nXE570DoQiaioSQCRRAXHNKsd6KVclxf61UWQ1pUI2qKNyl643oGtKlIgUFmrZaBNBVclfpQcHBGK\nArjko9RQKOtpUjb7UQi6NJINUc3cAYqK949eFMoSnTG3fP67VhXOvFKEqVlJ+UjH1oF/jARAV/4q\nixd1RJmIFAs9pTOo4nFEBDmqffFUDfJVqJ8xmfSgAoHWB0O5HaiLJW6FhCVR094PUUVpNpKvaMmo\nH2gANvf6fpRVFnXI/DMUGfdqQhxRUr5p22qoC2/5gE51ddqC9ys8oAdNvqelBmXTym0qSMjH+RRF\n2w5y0ryO4H5VVQC4QYEVAs684Faik43NUAuLkFMpkDqP0ohdpxsrMkjEk+1UM2SgCdSsmoHC4NOn\naO9FWW6gNDUftQK85YxHlO1EWenSM+5FAFRCEY69KC7DaUgJI8x36UU84ICe29RByuWvXrQJrbW4\nSRIFUV5a5iIoOcbjyk0A0jzCNutUMNIOtKtx/LUDagQjAlRNAutQSr+tBbWhKMd80FA6NXfqKC+t\nByPrRSbyhqMQKIGFEn2qhxpyET1qBW6XzZHQUCQR54OPaqG0pgDtUHKV5imgnyz+poKpSDHagl2U\njbpVC6oUJjA3oF3dCdtvzogHORP9KokKMyN6gutYO1UES1/D1Cgvkok/eoLMOGY/Oirv69M79hQK\nKWVD1PSgatER5jv2oHQqAehqCoWNRzmgG6NQ2z1NFLtpJcic7RRDZttqKZZRGMUBkpSk439agklQ\nTPagXcqgY0bdaDlrRpiaAbS59B9xRC98iB5Y9xQZwGkn1qj6A8VhopJlRJIggj8vSubTOccfLQA3\nTJ+lULtkpcPbrmqh1F0qOhn8R6e1AB9Zc+ZQ3kzQD0aUg7jtVEKSDMGf9qgstaQlOicDB60AFyM7\netBpW74abTJ1A5IqCV3bmqOp+9BbnLJ+bB27VQpdtuLfGJRGRP6UFW5QCk5zEDfFAyCVDzYT6jtt\nQK31nIDjR90+nagCpWi3CpnUelUAS+AnVgzQddOamvU0QqUtco8wx2orPUAU+UY6VUEbWoZG9AdV\nwSoFWNQqCzZJVp6VQxPmiAZHX0qChVsjaY1GijuNJSgDSPVVAtqCnAnbpPpRGg4pQEbnEGaiulxJ\nGo4iiJS7FVXFxMFUZP2oBKUNVEBgJJ/KgYtbmcdaAvNCTMwntQJ3is95qi3PlsA1AJVykJUZ+1UK\no4l8wG3SaCqX+YTq60BEvQCDQMJfRyzNAsSrViiIK9Kx+dFMBxOFR6VBIMkntQVUIj1oCJWmR0oK\n3KiB3FAmtcIxuaoTfnptQLlMj5YiqgrTip7mgbQ2VjIoH22YTtvvUFLkJCY2oFEKKV7UDPOBSEqm\naKAAJJoJYe0riPegcL6dBAFRQteZoCJ83egulsM575zQED6NtzQdrIBV+e9ARpZ1A/rQG5k4oAvq\nSgTuelAuhepz070A7keUxgd6BZpZTtJoAXb6i7M/WqhXWVKor3Dq1EpgQkZxWFW1pLW3v/m9Akl1\nKDkT/wCaqCh1tOop3OY/WgWcdTrGFF1WfNtHSgLKwlO/06nv0oBN3Mrck5kyP7/SgK8UpEgCPw/3\noAagEDcgYzREs6grJxEp9TQPs27jisqCDuk/37UU60yQEp/l2nbPrUAnUecDsPrQUUeYoaYJOMjF\nUc7qSCpKPm2n/feoAvvwmVwqMFP0qjNeWkp8ghreD3qgaGxpnTHpQSA3pOTM570CV1cNpBAG+xoh\nVl3UrSoxO9UGaQEq2mgOu0Gk9fQUEssktyoRFQECTr0d9qBq2Yb+dfegPcltadp9qiskQHVHTjYG\nqh9sEtQMxsaCNepUHfpQSop5ZJHmoElXCwBn2FBKHtRziqGGwiPXPSoKakpdAHXrQNJbUqM4NAO4\ntV6iZ1CgSfPLGo+XoKozXnJUZMRkJoKtal+yqAyUFPzZ7UDSQVpHYdKDtUg+h+lBbc9KCCgK796A\nhTA/pQEbSQO9QWcQlKJGVH8qBd3WpvV+KgCt1xQAPSqgJC9zmiuQ0pxwJj6mgYcsuVAIz1NAL4Pz\nymgeS2NGN/WoLJJGKoVuPMoxsKALaFqWKB74SAM71Am+y4lU/MKoq22rfY0BkyVZFQMIslFYMeWi\nmA0UriKCFpJxEzQCFsoE/pQWKtIjpRF0K67j0oqvM0qiguFBzBz2oKo0hUYoKvKgD0NAndEgApHz\nURnLaOr+tUCU2qaK96hvW3q6p3O+O9YVS6VpRpRkTgmRQIvI1iR7Hp9zVRXX5lRhJ2/SgEoHmTOo\n9PaqHFOQjzbGoFFMK1k9TB/zNBCmXupkD/M0BXlaYTusHUfSaDrNDnMDmkkdJojZYXpMIOPTb6io\nphJwrWYzuMCfailL93+MUiEiB9aIAl5YAJ8sHBmKout7U1iRAoElKl3QkzG+PrQJKWiSk98kflVE\nhUJ0k+v1oFVhWvJ3zQKXSdOTkZ3qom3tzqSs7H5lUQ2PK55ZKPt9aKuy8orKM6ZmBUBluNhU4HpR\nS7b48yp868T/AGoh5pKkIM0URWxVHTHTJqDLcBbJz7CqCM3KeXpVjMhVEcl7KlTHQGiqXN2oIg56\nA0Gc5crO/SqgzTmqMZ6mgdQ+vnJHX0oHA0yVpcPzioGG1tqmNx1FAbmDEnFBS7RbPsRpB96KwLjh\nIUklIlzpVQrasqS5ocxQNFABjrQQkkCAfag5UJIxQXC0qWcUDQiQFD3qAnJJPmEaaAhQAJooLgMj\n+Y0RKm8Dv1oKrtkJn1oJt7NBMK2igbFqwBq29qAbzCFKHUUE/Dp07YoBr0BMDEUFDrDedqoUUgb0\nDLDaUo2kjrQMNkEx1qKCplb6tP26UE2/CXVKzAigdbs0plJHSgkfwz0oKulJymgoVaY/OgCt06oA\nigrhStIGaCqlBvp9KBY6sroJQvptQQpwJVqmDQcHFLOrv0oKKC1GD2xQAU0orgiiFnoTqBz61R6y\nzcRokAhO0g9t8HNYaXcJAJnygSZ+1Aqt06u/+qqgakoUrVMdSBQQNMSfzoLulSlJQTA/lkf5tQH5\nYKdUD/SAaAjKEqVA3A1R7UAXlJ0jprP50DKQnkgJzGaCqn9BIUQe+aA6rwqGqOmR0nv/ALUCiVuO\nrXrwJkp6xQXWW1BIUsJQjzEgbdhQSHFFCjr8qR1I26e9Bmi5UXFkmEr2G0fT2oKPp0KBGdWT6Cqi\nISRI96Koo6W/MnYb0CrgS9p6RsO9VFWnFpJR8pO2aAPxELz3OfSiGUK8oMg9aKKsJdgpk++M0APh\n3g4An8KsHaoNhkqJ0aYCRmoq7zitKQPp1oElt6lCTM1QFTSojoNqIWdUtKAD06UCz6pA/wBVUBTr\nKilZP1oHG2FAAnYfSg0G24OoSAdlRUFy5y42VPeg74hxvzdxkUEC5Kpjbr9aAxdUlqR9qCra1Zna\ngq/bJU2p1PzjagUbUlfmGBVAXlS4O1BC1pKetBe3I360DTbwDmTmoNPWnk5wrvQA1pOKAejMgdaB\nltqczgdKKEpsEyB5poiUA8yOnWgMhvUcH2oqXgGyB1NELKdJkR7UHJSk/NRXOk9NjQUbaBBVsOtB\nGEqxOOtEGGkkKwCaKMysI3HrQHQ/6RUFnVgJkjPSaAC1pU2f5qAC9jtVC2pcx+lBUyTKqCW1k7yI\n2oCKVjvRCz/yztGP96BNLi5JTkUEKDhPmqg7alBGelRVuZJz/vQQtOQOlEKXYSGyNicSao32kXC0\na2zq6EKhOk5xJ7gVhQl3cfNPt0mqBa9cZ679aouFESjYjY1BRSX0ARsnMiqBp5pAJPfSPSiGWXXC\nU7nT+npQHlUqVOkHG9RSyLtwOaSAqYk+lUaC1IQ0Sk5IwkGoF2glZGZjOfzoCLGltPQThRoJCm9A\nVqjFAF+4bjSfST1/zNEI3F8ArlhX8NR6/wBaqhso0grVBHSiDOPBe31G9FXbUmZ2IGI60FXH0CW8\nQroaBYMpLpUn6D9aIWfZlZz5oxVCzViFZH+RQaSLXSkCIO39qB1tJbbEx5eneoougK858oUMY/Ki\nKtvQkqJztQAh10HROlNFQ22SBHzUDC9LSBqz/eoAPMJWgwBVCSmGwZ6dqI5DEme/WqHktt6YMSkV\nFECITMzGJojPu1KK9sTuaolJK1480CaIKmBCqKaSAUioIUU7jYmJoLJUlKNJ64orNdCm1FKfeKqE\nw4eYdX2oDEJKdqDm9QSMxFAeFHKR96B5NwdISYBH51BbyQFdVbx3oDNo1K1HbegKpSQJ2PaghtKS\n4STvj0+9BBaKjIGKArSD+I+WgjlFRTnfrRQ3GWkqMHNEQoJCc0CroyP0qiQ4YgbVFVLgTNETbyVC\ncJOaB8N6gSZgUCa3FNriZTRXc7VEnB2oCodQUZ2AohO4ejKTAooQdCjGqKDnHNh94oiiX9IPTuaK\n5q58+ciiKOrSpefloKBwA5ECqODiCuoq+vybz6UQHWAZ+9UGbclOcxUA75pC2Z69utUb7Nu0lqXB\nq5hICDOI/FisKVfb5QBSrzDO0RHWqFrVbhcWCJCsiaAilFOR9P8AO1UX+JUVctUpJ+8/lFBZDC+a\nATpBEhNQPpbSUggZVgx39qBe40pMRCR+LfagzFOfxwoHKcqH6bVUaTTqHG5ODvpqC1ulCCFJHmPz\nTnNFEcSC519B/UUEOhJRp3SNv7UGa+lrlFUQVbRVGdrQp6FJn+WiG3VhDQA3FAIuj/tOwmaodBaC\nEhO22qf61Arc+VWoDH5ZqhhH8C3T5N9z1+lQZrr2pwEbTk1Qy0qNKh1IKv8AP6UGoy0y8vzK8mdP\ntUBuWgkdRGO+KBd+6Uk6Nk7R170A1BKZSr5hn70BmwUtBMwU0UC3VC1QM0FnSdMK3O1EDDoGRmgq\nWgTJ2+9FdPl0JG2aCyhKTAzRAlFyCT7CgUUhSl6jsOlVFgl1OUYSTt6UU02FaYjzgTQcXupgHbG1\nQUBKljOJk1Qd5SSANj96gz3nltvb+nbeqFl6JQsbnc0QRLoK9RyBt60UcAqEgdO1AdmcRsN6C60J\nVEb9KgYCVGOp6UEuXBRQCW6rrntVDDZwNWJ3qBpD7ewGO/r1oLpInzZT70E8xJASBttQBup06gI9\nKBVZOnJAooHMxmJqosHBG/rNBAKVqxvQOspycDb7VBKngkeUR3opS4/1COomgUcUQjQPm9P1qoqg\nqzucZ/vRQXrgr8ic9qIXN0UK9aCwugSDO+80BSnqo47Cgpv8plQoLNqnO/egItKVIOfaoAEacDfv\nVEqkI6Z/zFBVStITqTvtQMNFSs/nUUJ8BLmj0mqPVsWj7zYW0nVnqofMo9JNYVk8UUURmNEgE1UL\nc5PISd1n5/vVBLdtSgUkxqGPpQVLcK9f60DiEuIZ5iznp9KgbYu2lMbz/cUGdcXmtahHtjH9aqBo\nYbUQOg+Y0U4loHRpx7dKgM3rMjcem/50AllxL0ap7T3/AN6DnVLjVJHeNvegSdLfLVAnEj1+9UZh\nV/HlPzDpVQVR1GVHYfnQKPKcCtCTq09o+tA604hKBGZ26/eiCttkuBH4DjHaii3K1lHmTpSMD7UG\nY5oPyn2NEO26AYKtlQU9sUVpsJIXK0+QgbVAYoBQSBp04+/SgXeSdEr3EEH8poF31xGnJiCTvQVN\n2pYCQd9zPT1oCNKWhuPxDvQUUsqVJxFBCULKoG29RVlpWGyM96oAy9y3lIV12zVQ2lWnUJPvUAVb\nZOD1oAakYOon0qo4OpV5UnftQMhTafnJ/wCn/wA1FJPDTKuk5qizCQRKfUmgPzExtsJgGgzFOqWs\n6sZxNBIAgf5PeiJYaLilqjCYorSZUAgA/wDk0BtO507mJqAMnWIVHeqGkP480d9qgBdqUkiqE3Ll\nQMbCiGWHxASsH0FFPNIkQJSN6gN+LT03n3oCF0oT5cTuRQAef1T1x81ArdAFKYOe1Arp8xn6VRx1\nRjzUEJUuZkj1oDNvxIn0oD87UgZ2zNQQGw6nUcdKBR23PMViFDrVEJbKVSv5YoFbi3P/ADEmBtHW\ngXfbITr3HWgG1oUZmI3FA0HFdsd6AwQNGpQjGRQAEhfcUBdaxMUCbilpX5iYNARtYUoSYEUAHlEu\nzO3egctXUpFRU3OpWpZGIwYoPUMShEzg/kfSsqyOMuNrQIMnt/eqhBtxM7xVDTdwjoCDQUfWqCRI\nE5J70DLtwTaN9sx70Ec4oaTIgnoKDuYW1pPLlB3UelBKn5bKkZK85+00DFkrWjUnMbiMRUDAdSmY\n6delAmXXFO5mE/pQDvX0hrT+dAuhY2VufyqhJSk8xagYSPziqiji0wIO4oKW4gx+Gg0G0fw4xnag\nYt0pSTGVnFQC4ktaUQBgVRmIkwlXfNVGskNgpCN8faoreZYbNsNZ9qyE3pJ0hee3UVQstaQ2Eb4M\nmgTU8gFOr5weh6UE21vOSIHrQHU4EHf0FBASgqKiR5uu9Ay0u30lWMdOtQAef85n6VVIOOzdtwOs\n4qodX/EVANQBfWcoQI7iZqhLUD83mHvRHBOh0D+bJJoGlOJJ3hVFCceKgARnr9KgFrKXSlEmBJqi\nFvlLhSdyMwaBUvFDwgxH9aBidLIO0ZoiWVaUFSf8mimLcrUsbRO1BosuWxSQt9tpR/5aHDGqN49q\ngVUW9RHxDBkkQFH/AP5qi40agBdW+kbnWf7VBzwQ+paRc2+f9f8AtVGf8C4CVfFW0er3+1AawtnE\nqcUt5DqR8mheuPp0oNJT/KhCe2/61BfnIJKu3rQBN15wkKhMyKoEu4ClK8319aAL7/8AEwdvLQLq\nfJIOo+lAVlRMmcnegMh0AweuKAREJB6Gglu5SEnJkdaCGn1SRBjeooqFheqCQO570QBzKVD8R/Sq\nFy+QQDmP60AH3JmCTigFbFQUcf70D8pS2QTJ6e9AZqC3BOw61FLpHmMZ9qqI5nmA70BHA2WCVbja\noM8nrmKotzgViMmM0DKEKHm/Koqy7hop0OkjvVG6xcJW26pZ8yY0EVgZt0NR1zg9/wC1UKoUjmbl\nMbTVHteH2/gJnw0m94w8p24ceCW7WzJFwQB5itSvIBWdq48g8pCnXVIlu1KzywsyrT+EepitIhK1\nuFthJ8vRPuaII9dkqSkQdOJ70HLutSTrVvv7UAWdRJXHkSc/XaitaxfIlI3Pb+tQS+7okH5v5vWg\nVVclOZ1HbNVCt1coOCNJ3Ak4oEVXCgDCsnY+9UDCiBGn3OaArTalDUQeuIoGbZscxMiNWCYP0oGl\nJ+HvCwtQUkbrAUobewNQGZSpTqUow44dMnyj6q7UDV54c4062CEJKFGEu8xAb+5/WKUeey2pTah5\nkkgx5vlOYNVGhaXCEqSo4ExJnaorY+JtlpPmPUJPrHWgylXzQeCUjznef9ooOcukK1BQCAd49Pvv\nQZ4fCX0qkaQetUaCLxGnT1PSO9RSi3lKdGQkf5iqjk3CckGUzBH61BdlySVDyo2NAN+60qJBqhVV\nxkKJ23NBpc4kJny43HWoEXrg6y2VnOKqBEtNJmfMO1ABN2pS9Wx2E0DzSwBp6nrNFQ4/tqO0R/vQ\nDUtXMUR5dP8ASiFjcKU4T91daAaF+bzZ9aqj3FwoJSOhH0ogjDidIScA5qKZbuETpBx0FEBed1vM\nKB+Q9O52or0V60lPhS14ijV8Yu4cS7ca/MZCt0zOcTioDcUsG2eL8Jt2W/4N4y0VtJWClZzH4oGT\n1oD+KuGM8P4nw9LLIYt7vkodQhYIWNLWtKoUfxzvUxV2+HW7njV7h6rVfw4QtxFlzB5VqbTKtWqJ\n75qo8qy8UuFKnCtSRyxqOwST5R7VUGN35kHeMUBE3TaU5GeiKKSN0ELnVIFEUbuBk7ZwaBf4j+KQ\nTjOaouHkrHZKeveoojd6I/3oL8whOvOn6UFk3KdUHrtQV1tpcUFHbegIi4TEJ6dqDQ4S5br5zN4y\nPhS2V/ET/ESvpFQIXr1q28tNu4VtbNqVgkeoqjNdugVBPXvQX8nQ0ERoIjagOUrTMiEjMdaCouSo\nFHbc0BFOkJCU/WihuKSnRO5qAdxdAp0JP1qoXLgmB+dFXSQCImiDKfUlPkBM0CzhWoyYT/3CqPTX\nFxw5V0v4RIXYoI1rJIgn9aw02E8O4ShE3jSmlQCHFnSkg9RUDKPDvAXWBcMpDg/m1GPSpRCuH8GW\ny1bciFoWoIamZncwP60GVxxXh3h7jbTrKnHymV6TAQOgG9aQvw278MvuJGgsOKHzLJgHaCR+tA3c\nMeG2ApZOvloKlaSrHp7mgS4V+4uJLWOQLZCYLZWs6j6VdGk5wbhNuwt5XkbG6io9Ov8AapQvangl\n02gh0oW4fK2VGfL37UFLj4T4hVrbWTt0/BKUoJJMbn2FAfgnhy5vHluXzXwlulJPKM6+yeuJpuh1\n/wAKWiiXUt8ySdYzP2mpQFjw9wXnhDzMA4I8+P8AV1qjnOEWFoQ6wNLqQFIcSomJ26/lQF4Va8fv\nHW+Q48pC1wtxJ2j3NAF+94k0+tty5c8ilJjV9qD0/h66sOIMch3iDiLyToaQuCQESPfNTTBb6yU3\n4S+PU++bxxMqQTKfMYPl9qis674bZt8CtOII1C9cdIdJONOYxvVRsL4LcN23D3Q66k3jiEOISpIG\nlSZMDpUoy+P3VtwniZsheqMNhWlRSuCTnPsKuAyrbihat1NtvqLykCQ4z+ITEdMGZNBncabvGHXr\nW5LjTbiOqkLUEHG6Z3oPNL4TwZXkDz2pWyRn/wDWtIGngfCNTYWu41unQiN/riijMcL4JbcVtra5\nddl2YQcExgdBvQYviEcIYublqyae5jUJkrlOoHzE+kGmIDwvhS7pSlLuPhWWkBb6zgCfwCYlUZqj\n0nDOGeDnyhm2u3XnPllSlCT9gJrKn3PBnAvmUpyB0DlKFnvDnhpChqU5KVaAOYcnoIFUFe8OcKQ0\ntTiFhtKSonmH8ImiPMu8MsG7f4lx17SoaoCUn7d6oz23uHvvNtq5qWSYny6vSg32/DVkk+fn+4KK\nlFl8FsBhHxGVaSrywMUE/uCwUga03J9AU/2oLK4DYDBbuAOh1JP9KCivDvCwISLgr6ypP5YpQtdc\nH4RbguKFwhIHzqUnSO2Y60HWXDOCcQYUtl26W20dK9vLOe3U0Fbmw8P2pIXc3CSndJ0kgfagBdp8\nPWpZKF3NwpY1FMpT5SKoVccsHL5tNqlSEpUDK1a8KTkRAjScUG65w60Rwdu6BV8SVaVfyRJH3xUF\nrrhVszcWzaNeh0S7j2yMbZool/wiytl2wbCil5NupwqERzmkLP21GgMvgPBkeJFcPW8tFgEFYf3P\nyyMZ60HkeMAtr5aVRocUmR6YrWMmODM2jtv/ABbvlvIMcsxkU0aCuH2RbOm7xIkwCBNRSDvD7Fl9\nDZuvKsGXIEgjuO1BCbSyLYK7stmTCdI271UEHC7TBTdmI3096gH+7LVETcqgjUPL0qqsi0tFq8l1\nJGCNNQV+HbWpKOfCD+KI69qIl6xS2jUXdSU+gk/nRS55aWUOG4PLOyIzVQmL7eTmd46UD/Dku3TR\nVzAknKT6CBG/rQFXwC9V/ELsgicI/wB6isNx3lvKSk69JgKgitIZQ6wljUX5d25KUmfedqgasr1k\ntHUCpZlSV9Bp/wBqCPikl+SuSTkVFWN4EhQahQAjUZ3mgVTckOqcMyBhNUafDUt3TiA6BonUpPoN\n6yp7iI4HZnXykhagDb2cFRcT1cUr8I7UC9hxXw8+4lD1vofJjWrKT2lXSkC3HLYm+C+Hr5tk4nOi\nYQtPzp+hq4hZFqyUG5uHTyemn5pHvQJquLQ7BSQZ3I3HWY2qj73b/s84PeSU8LbbbRsYI1npEH86\n8/TrEX3hfgzwFq5aBAZMaCpSsn3NWoBbcLskNFtNoA02rSUjv60GP4p8PptuDfHWbHKc5itakyZS\nN59M1c1Nx82uAt5Lmx1w4pxQ80jolXROciujJZhbKAtS/wAMaQcHNUHev0LSQQM7T2ohRx2d/KOs\nd6A7vEOJO2jaHFLNoDEq+THSaBRT7hMN7zAj1qj7R+zHwoqz4Uu8ft3L254m2FsupUlOhhs+fTrM\n77nauXvWsxoL8KcTHEFcZdt3tSG8oBSfIv5Se5AFTpYRfumbd3mL126ZEFwRBOwkkVUZPicv8pm/\ntXPh1hwBTShHMx6ZGKuJryzfEOKC+cVcYISpLbQ2ztI61pGIOMcReuQxrUDOyMQPpmtA93xG4t9C\nWjpQ0IH+qagDYKcUwu4W4SpCk+TYgT6RVR7Dht6h7hqwWh/6fzh5SsrnIAB6zisqzn/FCk3KkC0Q\npnSkoj50nrqVnV6CkG/w4t8RtbS5a4cGrFC9F1cKcnWf5UD9TUCLvgy+d4g49bpQu2cOqFKlWnoI\nFOlj0CLO1YKS5baG0Kl7z7o0gHr3FQZ7pYu21XHD1BNuryp33Pf+lBkcPt7y557jakjQotBP80du\n1UN2Vg6tkKvv/tuhfI1eUp23GTt7UCPiocBuW2Lh0K5zk8p1nKvLjz+nSrg8gxbNu8wvPBlnzAmZ\nX/N7maqESo8xMy63/Is6Sr+tUev8E3ALb6lXOhlk8z4QhJ1deucAdKzo9G7xl1/hKb+1CU25IK3H\nAMJiSB3PSorxt1x25auEXCFjmQpJJExNaQKy4/fJLsypxeQ6rpO8g4g9KRCt5f3N1aoAbUssg63Q\nTgGqrJbuNKuYgnWnKT2oPc8Mubt7h6HXw/qUMwlwiO/sazurDSkvvyeXcLOofhXnapSD/um+WRpt\nbgpMmNK/elWId4Lxgb2d1/8AFdKkTZ2HHGHVqPDnF2xQeYtxGohQIj1jTM0pHlPERuWL923Uhdug\npJQ0sESg7EitYjL4c5dfEIatiVuPLH8JP4/SqPQcS8OMJbI5jiLnXlBEq06ZODB3xNSig8FXBs03\nTBcWQpUogBSikiAB0xUqmeGcBuE3SFXHBbtbIILml9KT+lKPRq4c1+6gytviikhev93m4a0BR/FG\nmJqDJ45d8d4c+3cWJvyhm1LLqn3AvlsqypKFJA0hMVRLviZbrbjFwq/Wi6ZRoW4+FRC0LlOpAzgU\nRo2C0X3iFtPxPE0cVWDpui82VgJScTonYUViueHbPjXEvg7V11l5t5arhy7WjSQJ1RpHzEiiK8d4\nFb8A8PPMuuIVxBT4QhxsBXOYHmDiHPwaVGCDvVweX4TxS4ZlpC4Q7hcCT6HrFa0aN89c26bdD+lP\nmxqA1f8AyqIiwteHXl2r4hl5DEJOhChA8vRw9z0iihl5thty1K8g/wAJah07K7UQncqSl8pbe5gU\nMKjTFUPss8YYt3FGydUpwR/ylEwQIIPt1qVYRulvNrQl0LZdjUsKEET8vloiH+IOvMs85etOxnfy\n4HrtirB1zdMuWwbIWFT5Cew6Ugpw7ht9fvaLZOAfO6r5E+5/pTdV6pduiytwyt915xvdUISgk9Tj\nAFZGg34lddR8Iy2pq7LZCDgjUBgAY1D2zQfPHuYblxtaShzXCwqdQVOd61UQjkBK+YhWtCZwesgb\nR60G1wPwnxm+UDCrRk5HMGSn/p/vWd9LHrmOF29otFnyE3FsGzquHUp1lwq2ncJjasq8rxvw4vha\njdMPEW0yQfmTP5GtVC3Dm7d5t1bt3y0pV5CU4Xn2MUBr5i3Q0o2vEuavq0EacdRMCgOzwtly1aub\nlIZa05uXVEfbqr2FRSNxc8GbSpFnYh2QR8S/qH/wQlQj3JmqjR4fxLw9/wAO33D1By14mTzbRwyt\nJViW9Y2BA/En60g87d3D0ckmGydekYBMRNXAK2SVPITMSRVH6+c8O+L+F8Ft2m3m7hSFwtehR1IJ\nxp7RXj6dmXf8C4griSn31BNo22mSlWkayfMJzqMbAVc1IPa+FWXOMtXDzDx4e62dVrqPNLg+VR6B\nIFKQS+4WL1N5aMsRwtDPIQheQtzVmD1MUV80u/2TceRcFNgptGtUoQASEIJ2Kj1HaunbHLVt/BHi\nDg9y25rs13MRz1MJUsj8WTiY6xUuaRfivh4cPt1vrsGXy2ku6Dy9UeuJ9aDyCvFXhpJ1t8OPUuQy\njPsSa3ylXY8Z8Nvrxu3TYLDZVjUhspTiY0gU5SsriVgsvvXDfJ0lYWgKGlRzMAgVcHueE/tDtrRj\nh9p+53Qpq0XYqcbdCgEuKlRQneZ71nfK0fxP+0Gy4jwR7hibK7QlxLSOaClKk8rvHtUzybr51fX3\nF795NzcalKt0kNJWiRoVgq07SB1roja4JeWFlZpQi0Nw4tgB1x4kpU8CSDKtfQ5ipoK3cPrvibRt\nDSn2xzVxrKD+JULEH0oEOH8N47avXaG0NhGSLhYZWpxRzG0jeg+jcF/ZRYXnBmnfEZZcunE6kMM6\nU6Sf5lpAmue+m8x4fxt4NseC8RDLNv8ADEIkpZlbSkdzOZrfnWdx89cu3mVkIeXyz8ikbfaujBJV\n0dchXlmZjrVG3w/jvEGeGKtviVNIAOhmfK5r9PzqD2HA/EF1cWbDQuRbpaSnmvFOnX5YDc9DPas7\nivUoRwp1k8S4s6m54LaKRPDrZWt110ny/Ex8rYj5evWsa0Fxji9lxd17ilu23aOqXyGbRCdLdzbo\nMF1S/lBT+E0zIbrG+P4U24pDKwpQwdCVKHrlAiqi/wC8rNSsqCEnGtaFJ37kiivm96yhq3UFPa3E\nuKKYz5ZlJHX3rbLP1uLgCCvqepFUEtHVpaftnm0lFyjlpcWJ5agZ1gig3vDlnwRlDrt3dQ/bOBbK\nvPLgHTREZ7zWdBF+L+LPOQ4kKtdK2xbFsJQNX4o74qwYRQtx1S1JIKpKEj3x0zQPnw+78At66cFm\n/CS2l5UKcT+IhG/tQKscOLVzDauaNkHICv8As+bNBv8ADfBnElunW2Cy4AoDRqCD6DrUqtVjg37R\n2s/vJyytkqwkvfgB2CemKdJHqbVxbatXx96+e7jv/wCorKqWPFnr61U+VOtJ1qQgapKgkxOaAgZS\n+dPNuHFEnqP6Gg8N4x8S3Vm/+6bK4KUaYvFK8rgUk/LMnetYjyPFHbu+u3b15zW8/wCZZ2z7Voex\n8H8Cat/D95xhPKfvHQW7RwydIwlzSnYq829Y1SJ4RdL4g+6u+GlxMagCDtB6be1aqPXeB+B3Fxa3\nFu1dNww5r1uJWSSvsfLtprO6rUTwy7VxxfCA80l9LfO5mg8vTgxOqZzUVbi3Cn+FhjmutufEr5YD\nKFKKf9SvMPLmoaw/2gW154dskpuOTdtXYWytKdSIx6956Vrz9TXhbrxfccR+FRftJdatEFm2CfIp\nKVEHJ6kadzW4jas+L8HdbN426+1xJuW03anV6UEjBKkeaIxWQK5u7cLVc8JeWpc+Yrz/AN2ojPmz\nVGbfr4rf8oPXaf4S9aENpTlXrAE/Wqg7PEeKMgyGjkxNu0fmMxMbTUVncXFxxFSC6AhSP+XpQlse\noMR7zVRo2l22xaNWaeSkI+Vy5ZkmTmVCoMXj9+1cvJS1yzoJC1tJKAr1O87VcGe00+5oUlJPmACo\nxM9Sao9nbcY4u6FvOXBStrywEp0kAbHEd4rKvMcWffev1uOu81QMB0gSUxiYrWIVL5UkNmCkHAig\nK66eUhZiYIB+u1B6fhfFblPDUl9VvbWicEso0KBPfJ371lUG/wCGrgNPtkqUNIGZT+Ix1igpdPcO\nUk8x9tCU4STI9jPftQefvkNsXLfJeFwl2F6k5Iz1NUW4bdFlbr6m0r5aVRMhR9ARTRs+Hr684nxJ\nTC3/AINBQVai6QDGw8xyfrU0egVxSxZYQt0OOul19jkB5QWVMK0pI1fzGsq7jzPDLjw04pCX03jj\nSStpxyUoUfwdJimDxdoB8HZoUT/H5qJHSDjPvWhNo3pcU6QonWhpKVZyT7elEfTeHeFfBnGvi0It\nLxV/w4stXCnLjBDon+FuEj6Vz2tvH+OLDgHC7h2xsLN5t61KBdPPOa0EOp1I0aa35ZeLUpufLv8A\nX+tbQZN6sMKY/wDtqypBg57idj6ipAWxdZ5iCGjrSRmR9/Tag/b3FW0upas1POanARykLyoAAnzV\n43YitV2lkuFtIC0hFtbL8ywpPRPSe5oCMMuNcNSglK3lJ/8AULQdJVnp/pG3rQV5hbtVOXDgYs20\n6pIxH/8Arjb9aozXvFZRZoasrA815RDfNGjVGy1DeKsR5TinH33V3ATcoc4gwysulpEhsgEhMqMd\nM9q1Er5A74l8U3KF6uJP/wARMLSY8yf5ftXWMBWXBLi+4fc3bA1tWxT8SSrzAHqEjJApQDjHCbmx\nLSVQee2HWXUrwpB6zVEcDPD0Pv8A7zfUhJEo0ySrvGDQaw4rw63cS5J5KMtrWjTGds1AZfH7J8la\nIUT/ANS4/wDxpCh/FvPrUghKW1iAQnTI7GaBW6s7lx3hhDoS3bLK3kpVCfQY3q1DFy29+9U3Ldyp\nNqg+ZGlUkT60EJsHby6fFshxxWFfMlAHvNKPROftb4+ri9nYENcPtbd1DVwWyXJQMKmf6VOF6P8A\n7QuN8M407/8A2u8SscgILrkoELPmjH4RWfORd18ZfhsuhhIDKVEAH5vfua6shJhMuGFCCfU1UEDi\nlLSv5uqIxEUHsvCvC2eIcJvbm6ulqda/5NvhI16cL1f0rGrj2i7dNv8AsrsrpK4vEtLS+jWCpRKi\nPl3rH9aaX7GLFhvwubriSyzca16G3oSOXsPmzmp+RfKvEF2rvhTjB4aubpTK0shHlXzNvJ1qAj1i\nXf2XMMPuauIfBoDjSlAucyJgzmaf0/j4M+FytPKcKQdKF5O3612ZCSFEpLaHDg49R2iqh1gIU06V\ntaA0kHQudaytcQjHTeoNe4srq30qFs+0yv8A5K3E6ARH+rNQbbXCOBsWjLtxe/vN94x8JY+QJJ6L\ncUJ+wqK43XFEqVbcN4WqyHyzbtfxT6l5yVfaKDU4B4G4hcrNzxw+RacsHzOK7ajmCPSm6R6qz8L8\nFsCF21sG3B1G/wBzNZql+Kr8S6y1whq1Zb25rjo1+4TECmDCe4R40cUVPqtluH8S1g+4rSLscC8a\noRqSiz74049dqlWAcG4V4lueHtuWb1uLZE+VfzeY5Ox60o9pwjhgsrUJcd59wrzPO7AnskdAKzVf\nDPHPF2+IeJ764Zy0HChtWjSo6cGR9K64yweZJBB9yf6VpH07g5ad8B8OtnEKaWxf6VpclBUXyfk7\njTG1c/6r3qv2acLQDpW4OwCyP61jpqBf/T2wSJDtwEKxpDh37zTpIFceDmENQ3c3KHTA1BwyBt/S\nr0QmfA/MLQevrlbjQlJ1YnqUzt6VeiPK/tZ4XeW3DOG3L1+5dBlZt0tuxsuVatQ6jTFXxqa+YbnM\neldGVgdKRH4sKTP60D/DXoacJdASnpmoGw6iQUuTq69zQE5gSk5hXtQD5rYUNTgSDkr60Cz3EtMp\nC9enHmyPpQXaLD+yG9QzEf3oG7a+u7TyWjgQgDzNY0n/ALciijvcSeuB/wCpSNKoKVIWQpJT+X3q\nBLiXCGlDn2l4m6bgS0oFLoExtGY61aMtqwu3Q6W2XHEMSVqSkwAPWgcVY2auGodS+pV2IlkohME9\nFmoJvBptnLchKHSgJwRkg5BoK8OsXmlsPpUhwtatTKTnPrEVRe+uDc60QlDSlAhuP61Amt11wtsI\nZKy3GmN8ew2qqPbWLrgUsuJT5V6m5yCQe4qIVtXg3dIbd1hSTEtR+tUat/xVtlBda5vx8w2+65rl\nPsRAxUUufEnFbhss3K+alQiY0wem29IJQ0DwZlLr3KLTrpEx82PUdaIdvrN5tuzUQG0vnngEiVLk\nAxkVFemHiP8AdodVY3QYuXtHPCIBWpv5dZ6xHWpBl+IvEaeP39r+8QnktgI0t4T2K9I67UwYXifg\nVhYOhfD7xFxbr2bmVpPuNx+lbzUYsQf1qoPYqi6ZI/nEg5FNV+wOI8SeDLt60jQytxthl1RUFpLi\nwhSz7JkivG7Nsp5baUPKdW4RCCRK4/mj9BUCF7xjhHDeWq+eSXDOhtR+UDYz61cwecu/FPAF3TVz\ndXPxhkpt2m55LcdVJn8X4jWuUrzvifx3w5nmN2qXG1LVpfdbIKig9G5OK1nlN1gr8bcPHD3LGx4R\ncrDrZStZW2FZGfX71rlK8Y3Y3JKlfDrAVmMq+8ehrTL1nhzhxY8LeJ7jlKSfhWmgnTCvnms7q4yv\nG1o/ao4M2+f4yeHtcxyCY/yauGlWvBHHFlCg83bqdHkBUUk4nJGPzq9JCN1wO7srhly4cU4FLUG1\nHUQrR8xTqzAPWKtI2OEOcP5zh4i4tFuE+RLcyVTtj0qA11d8N+MS5btFLATpAOSon8R3zQYl74nt\n2nihhlZCSNRJiYqxCd14ovHXVBthtsKERlR/OrAdV1eKHlWuVxqSjE/agDb8C45c3AWzZPqkzqKD\nHvmlInjSOL2lmUOtqaBhJUUFMZpg8vzHUq1I85zKhvmtIovmLbE/KPXP1oKMkpMHI7VRu8JvHGwW\nN0qyOydVZ0aoVeqcSlJ0QYQPU+8Csqop7iRP/OUANiRQGbuuI+Xl3EgbkpjagN8XehtPMdRpk+aI\nkkRB+lFF4XacQu0NWtnbOXLTZlCwNCAU/wAyzjrUD6eE8F4SrVe3QurwmfgrdyEicwpcffNB1v4t\nybc8PbQD8qbUBavQEnUaQPW3BL3izgU1ZLt//wCR+6JKwP8ATvmpR6S08H2yUJ+KXzWmtgdIGBuT\nUqhcU8VcA4S3y0vG6fGEMMQo49dhSDxHFfH/AB6/cUGnfgLaPlZMrPus5+1bzGWFd8QvXlx8Y+7J\nwpa1f3rQcY4XcBoOXd0bYEkJSSokf6iM4qUKrdvUSguqKNwNRA9/rQa/h0ccuVBVpculIErOs6Ux\n3/oKm6r6H4XtFW3B7ZxTpcS8gLUOiJ2rnq41Xnm2/MrCYkAgjeor5F+07gnC2rj96s3E3F4vzW2C\nCY+YREAV18s68fwpi0euWrYtKeduFJaQInSVnTIHetaj9EtN8N4Y8zwtdsm8trJkG2UqFkgJ0JJn\nZRIMVxbNJ4/YW+pd4nkWyGkBbjhOXVH5d6kKcHGPDK20cu8beUudCErkzG0A9KTVYTvF+GJ8Q27T\nly58G5bct61U4dCXdRWFqzhQA+1anxlr2lxwJdtzuYlKJJKlO40zA61Jql79jwXxPk2lym2vETzA\nhatQG41D70mlYPFPBv7Mk2brTlpZ2b68oXr0qjVhe+3SrlQxa+Ef2ZNNMJb4fZXOjDjs6jtMqM0+\nrQLTgv7N3nVabe2ZExIXpj+UjPWn1Kc/4W8E2jrS7n4UpVKQHnB5lJ6gd6fQoU+Ab+44jaW/DbV1\n2zLYRpgh0OJklEHpsfWrCsK5sPCVz46tuEfui0/daEw48HQEErRrlQ7pPliav8SvR/8ABH7PHOcp\nPB7UoQndMlOqc9YrP1XlP2h+FeAW9nY3fB+HM2/JWUXLbYOQseUmD0UPzq+dNeFcZaCDNolAT8y/\nMI962yBbXvDnTykWSC6cD55PsJqwWunQ06f/AE3ICRCkKCup9ag9R4fS2vhKAhsucydSGxGgg9QO\npHU1lWRxxItHXrVg/CJegLQtKjKD8wGMZFUePuHAt5ZGxUfStIs28trIWUTgUBbe+dauA4VBw6tR\nCxqSrpke1A8eLto5Zatm2TlOtM5HbNQW57ZI0LkTEfSil22mS2Fx5kjVI3naaAr9jw64aVLrqLpo\nYQQlTajvvIIoMtu3lKyl0YSVaVA5joKqPQW9oH2rBoBSk3TutZBwBoBVtPasqD8cu/VeN6VghXxN\nomTsjCkie6RMVRa88NeILm6cu7axWu3fUVoWnTCgczvSkML4HxW0aTdXPDlW6ER/FKUFM7QU6jg+\n1SjCuGHFvurbQAjUogCR9E4FaRFtYuXD3IUeSehWkx/hpQw5w160ug3KHFICVFQOPPMDMZEUo+2X\n/i/xCu3aQ4lHL57ag30GgyBXDMbrSvfHvidYWIZMiNUH7VOcarH4pxTiPGFW6rsgcpOhCkSBEyBp\n9O9XGWatM6mxy/IrMTJjG5qo8TxVaHr55zJKSQPea3iAtO+XSHDKvmb2GNv1qjXtGOJLgWrLygR5\nnEpUqB9NqlHu+A8LvUeFeM2zzdxbXHEggWy3lBI8kSZWU6ZO9Y3WsI+PH+HJubG//eDaGbBlDb1k\n28l11ak/ypTqBzvqNXymkrfxyi5W18OzeXSlnzoRgieg0iB7CkAuOcXdfuGg/Y3jF8+Am2+IVICP\n9KCJzVwHtPDPjC7gtWCtP8y0pQD96lwjcs/2ceKHkaX/AIZgKwZAcMekdadLy0rP9itiTzLxznOE\n+bSgI/qanZy3rX9lvhe3Kf8A0SVn/XJrPS8tdnw9wa08jbTTcYhCQDSkS6rhdqgqffbaA6uEJoPH\n+JfFv7PHbd23dvVXWIKbVJX+fy1rPOpuviPiRHCUXGrhTLotXvm+IISqfZNdsYYSFNeb+Ac7wqqh\nhpm3WhxzkEhsefz5E+lA5ZlpQS2pCwVmEQrf8qivRO3DNs2ht1JJAAR5t/yrKmuGW95xCfhrW4W1\n/NIifdVTRvWvhnhzbClcWu02i4BQhKgopHXV0NZqgO8R8G2KQux4cu+UwP8A3DgOkH1n+1BdV34h\n4+5y7Ju5atE4DbADTAx1Vgmr+g5w39mDJ89+6tKlbsoj81DrToj1dlwPgnBWJZZaYSBBcWYP1JrF\nqsTjX7QuGWWtmyQLx9OJb/5YP/V1itZ5SvJXfiK54whxPEb1dvbJH/tGsaz0HlzWojCVc2SFQy0s\neUpJkZP8wqopb2qrj/lpdCEzrV5cDeqDrd4UygISyp9yPO6VDf0HpQJpbtnVSvmkq6kpn86DS4P4\nXXxi8S0gOC3Bhbilgp1/yjuqOlTdWPqPCfDVh4f4I+xbA4bccUonUZ0ZzXPdaaHA7VLnh/h6RhHw\n7e/aKmqzHgosuuXNyVNNOOoSVDZKDAH0qo+V8aRYX3ELh95DrwRhHLXAUE/yiMTXRlqfs6sOEf8A\nFvDn02z+ppSnmwXNWhaUmJEZqejHqbDhNzxV9y/5iS2txZblUSgrJT1CgOwrKtC48M360RFvHzaF\nLURI26nIpVhFfh/jjOWkWKzugcwpgn1g0qRSz4Jx5NzzLiztFoVJUpl/zzER5gB0q0N3FitTYR+7\n1nbZbXQz/N6VAouy0JUTwt9ZOZBaVGZ/nGKqOTam7UArhbg0JybjRkdhCztvQJP8W4Qw4q1uGFM+\ncpWvkvBHbKtJEetAjc+I/DVoogJQ6PwhpDioAwBMCqglpxvw/wAZK220hpy2RrWXdTeCofLq32zS\nDyz/ABhqz8Tv8m3ZTauOIC1ha/KkjKkkH1rUGgzxrhTji7pkEMIUGgfLsMawd4zPepCvR2jIU0fh\nXwppzzqbQo5KtztHrUVu8a4w3xbgNzwo2awp1oBLilJgLTEKIgT5hWVeF/4AfcZIfcGlfzaW+/8A\n3CtdJCp/Z2m2fSq3d5hHVxMQfoo1eiIf8DXT7s3KirtpMSO1KNXgFpxbg/M5dum5S8U6vPp0xjAz\nnNQYviDh3Gb/AIm5dLsVeZIA5a0/h96oXt/CIeSVPJXbuKidSdWe+KUXuPAzStDQvpkqUqUFOk+x\nOZpRT/6f2h8q+JFI6EMk/wBanRFEeClfJ8fDaZ0qWyo/oetKGmfCNuyjQm+t31KgNI0FCpPzEFXa\nlGbxHwo9ZoUsPlEqCG2QkyqeozmrQlc24sbi5tefzijT/GHkBBEzBziaDPXaOMJc5ra3To5aPKQg\nFQ8q5SelWo2rC5Qjw5qQ4jnIQu2RJA0qWqCZI6t1FZgUi0Wxcoda1sHUpIKNtinG4Kao9K5xDg6E\nAurKCoJKPIdMdgY6Csha7u+EuN6WH9S8GdEbGgUYeY/A9qhQKRB7x/WqHH0sKvkuh4FSEn2lOYqB\nW7facabUtQbOAAUnKknP0qj7HxHwrxlxq3Wi2hAcSVqMJAHczXDp0g54Y1r5fkLmyQCO1UKv8I0t\nlaeXoCskK/WiMy84cEDlgaQrp0M+tUeLYS2xxdxv4Ru9Tq0Nt3CiEau50kTW2Xolr8QcOeW0pnh/\nCXEiS01bhbgG8511n4v0l4l4h4ss7sWl7xR95K20O/wzy0woSBpRFaxNwfgfB7PiF7w+3uW1uqul\nJDrjjijE9AKm6Yb8c+HuE2HCHnLNgNlHEl2jJH/8baOs9ZqedXXs/wBlnBVO+Crc8xSEuLcwiEn5\nu4zWfe/V8vGftN4Lcq8ecO4daIcuHXQkIQVEkic5ntW/G/E9YHw+/wCNcCunLa64ouyeZcXrYNwF\n6E//AGxoXqFX9svW8N/affNWqHVXdveAhRKXQhLhjb5SAKzyvRi5/bbwyxfCb22lrSFc61WlwTG0\ne9P8zp5ni3/9RTjqlt8M4XH8rt0v89KB/WtZ+NOnkeNftD8c8VtuYL0sFcANWidBOo+nmNazMRms\ncN8RXDlq7eNXSyVgFbuo7z/NVpy2LPgz6FtrUwtx4OOBTUeXH61KsJ31jcr5HO4S84nUoJcCjPsB\n6UowHuDvM8ON8WHUILpSlwxo3gCd5rVRv8N/Z/xtyFrYPnQPmEb96zvtc8th3wHxa0CCGS8rqhtP\nlTPZRNZ7WGrHgfFbRYWjha9fR5zSog+kzUpGmzwHxNey29qtsYIM/knFOlh9rwEyH9bgQp3q66ok\nE/8AQKz0R6i34LZ2yEhLaTtukR9BtWWjWoJ8qPKKqBc0mNOB+c0GJxnwpZcYe5ty69ERoSSE/wDx\n2rVSMofsy4OFz8U9H0p0Rb/6acHnFy8fXyinRA2/2YcLQ9JuXChPmSMHzevpTrSGLjwRYvgpN2+g\nYlKEpSD06ClITR+y7h41pTeOaSOqRq+9XpIat/2a8D0aHHnVp9wD96nSx6fhvDLKwaQ1bNpQloQn\nSNh6e/WsqpxxRHCL+D//AI7mf+2rgvwl5LPB7FCRMMNbj/SKmjxvGbDxVxBp60srZ4sOPvFagNAM\nr8olXTrWriQhbfsr44Rqu3W28bDJ/Kr0ctXhngz9yXrd23drXcpwj8Pzb06SPF+MOCXXAOJ6Pj1v\n/G6n0EakFPnyMGOta86m4Ra4zxZPl+Meif5z0FVBf3txRSR/6p7AP4jQSjjPF0j/AN28kxvrNAQc\nb4xJPxzudOJoJTxzjGqBerwo9B/agr/xH4hEH49RPm/l/tQcPFHiPRm9URE5CT/SgqvjHGP4iviS\nDqBA0p2PbFBRzjHF5CFPBcHOpCCM7dKDFfRfm9F86tDqivWU8tMeTppGIIFUfUf+GvDLyA7+57dO\ntIWgJ1AZE94rDRK6s/CvDGA/dWCGG1KDaVAuCVGY+U+lEX4fw7wRxZKvhLRTimo16XX2zJnuc0U4\nfAnBNAIubgKP4UvvzH3jFTrSYE54BtyfJxa5SpYGkc1ZI98VeiFVfs+UvDfHrlKsiQsqTP1inRGe\n7+zrjTZKU+IFuDvqP96dJHI/Z5xFJlzjjiv9AUUH2nUadECe8C8W1JSxd3Dv88vEgZ/CpMDbvToj\n0Pwt9a27QXzAG0JAVoS5kCM5mooIu1clxSHmwcCNBCsqE9MUEC+bdKUocSs4GSpBP1ITQee/aCtd\nxwdielwNWUqyEn+XIq+UeMs+LcZs/wD2t24hI/BOpP2VIraCL45xf4x6+Q7D76Uh06UlJgR8pBFR\nTLHiQx/HsLcqjSV25VbK0n/o8v5Ug17TjvhVxGhxty2XuQ6nmpnadSf7VItFCbC4uUt2TdrctqUN\nakLGoDr/AA1QaJVrm4suHsoVdldkFzAcDwGDG6UqTmO9QZDiGry6VcWDlu5a+XzhoL80ZydJ3qgP\nEuFrFjcLJszDc6UNQqB/KdRg1cHnrNrU6mVBKZEqVJGT1ArSPUptUMsuBVrb3SXvk5S3BpjE6Vnv\nmsKduOJcYup+JuX3SZ+ZxUf2or6b4O8Eceu7zh/GfhZtFWsIc5if5CB5ZnrXL1rWJP7MPEZ8K3dr\n8EoXLt3zUt60TojvNTtYb4hw4shph1H8ZthCCj+UhMb0o+cs2Clcet2lgqWbtpJSdp5grow+mePr\nPhiPEDrt9xi2snFIQPh4U49p2HkSNjNc/LevFftNsQjxVywqQ3bM5jTuMYFb8az6aHgXxDbu+KeE\n8LteEtJyAq7ccUtZKUmSlOEinrDDn7TGW0cISy1cc6eKXC3okJSsidGeo61nw16fQP2WWa0eBuGg\n7nWo/VZrPv8AaY+WftsevLfxnbP2bqmbm3Z1odSYUkzuK3+M9PmSSt967uL1fxFy6oKU84dSlScm\nTXVhd1GixIkStH/70Gfy1BMaehB+taQFDDi3iY6UGuH32FsraSAbfSoZ3KcjaoPVWP7SOOv3DPOt\nLdQQrUIBBkD1JqcjXtfF/FnLhuLRmZWoEn+asxV3PGHIFprt9fLLi18vO80ivI3XiG4vuHscN5A5\nNu8q4V1kqMpB9prUR7nhnjPxLepSm24cA1AGpWE49TWNzFr21q+ty0Qt4pDsDmIGwNc2xUujbBNA\nQawDjfrQd55lKdqI5RJ7zRVCler17VUQGXMQKAhSuO07xUFfh1ncmO1UWFqd501FcptQB6jpVFU2\n8pNEQGEwUz5vWgabY0400UVKRO0COtQY3iWRwHiMRPIUMVcRueE+FB7h9q64mQy2hI99IqelxtuW\nSln0FZUi7w64XOhHtVGdeeHXW0/EO9Ix9atR8o/ataXFxxq0SygqcbRpShP8pO5rr5Y15Rvg99vy\n1RqrVRdPCOI6SS0ry9KUSOF3wE/DqUCMYqUVVw/iGf8A0ywMQY3pRJ4dxEn/ANu58wzFKOTwy/1B\nPwjilAqxooBfuviWkTaOfKQfKatFlWN+EE8heopAmOozHvUAzZXsqXyFhMpM6aqOds7pEkNqwrzY\n7iivR8T4/e8P8N8Gft1tl8L5NxauQoqQEkZTuBjBFTBh8Q8UM8U4YzZcRtywsXDZN6yCpBbEglTZ\ng6hPQ5qwa/hDi3hSxU9ZOPlt642vSoll7QYSY+ZqQazqvX6XPhVOskqYbVpdebVzEg7ySmYqBFd7\n5NXPclSsLHQ1QdPFNKkNlwqO8qAzPpUgHd8YecQoJSjzeUTikCmt0v6FpbU2dzJmO1UO2l4hvy6S\nEZEDI+xqBK+4khtTZbDoMnMSPY9KBD9/a/LLuJMaO28mg5rjrihoSHcg7oED6/pQAPDrjiaHXeJk\nPMKVqCCNMaRCQNMbd6ULH9n/AAa4cITzrWAVHSdcDpIV1Par0PM8W8J31ldcgKS8FSphcEFSR+ih\n1q1GWeH8SZ+e0cSBOdJUIPrVC8DzER6Db8jQBdaAUDpOfTb1qh7h/iDj1inl2t86231bUdbZ90L1\nJqDVZ8Z3e95w6zu1K+ZwN8lRHqW8flUgrxHi/h7iFqvlN3fD7opMISUPNKxhKtlAGmYMvg9navoG\nu/atHtWG3kr0x31JCh9DWqNw8CuS0EcPbRxDl/M/bltefUJOufpWareRwvyqkDBTO1ZqxsfsTYuH\nP2nPNh1QbasrqE61adgBCdutT3/yZ+3zq9/eKbi4T8U/KXHBh5f8xzvW8TX6G/cN23wSxdUApZtm\n8kyZ0dTXn6+umZ8eG4MEWXie2uLq0F6PiUp5SlaRzFrhCsfy711ZjS/adwl5/wDaGo4CnEM5xnNZ\n8b8X0f8A2o8Dad8Xaz5nDatQJiNIjt6VPGmsX9mFm0f2j8OQlMOtB1bntoMfWt+v0y9R4v4VZvcK\n1XJKGhxW5LunMx27GufnW9eu/Zq1cp4F8OpYXbtf+1HUJVmCetT0j5l+2Tg9w54oWskGbcQegE/e\nt/jT0+aL4c2gOt5kkFJx0PWurKxs0PJDeiEQAM9s0RVPBGslAlRG0/3pQG24KoPuBSMzETtSrG5b\n+E23YXMlR06ZqdEel4b+zlDidUJHXep0vLdtf2btMkOlQUpI+X3rPS8iteE+HJTp+GSpzMD/AHNT\nfS5izPglTTanlNNW7W6ktgEn3NOjkp/6O0X5UOLR7iqjV8PXNhxe6ctGgpDzQlWrtU34Y9Pb8B0n\nR8x6VitQ7+4lpH/LJ/7hSrDiPC9y6hJ0AJOY1CjIT/h59hyVN/wx0ChRXJ4XaOuwhBT6k5mKA11w\nIo0aT5CB96AY4YlCUBY+aZiNqUCRZ26nICSlKcCIqoKjhLfLCtJOcbVKpa4sU8zQE+UHJO/0qii7\nRlJGkZV1NVARattvD1oGvhWFpUv8XrUCDjCjOM71RmcZaZHCrjmNSSiN6qNzgPHlWnB2W12yVJQA\nG0JMY9fWpuK10ceeW5AsUx+Lz+lZU1+9SlAUq1CR/wBVEYt/eXHFudaoRyUjAcBzPeqPiX7QuGcR\ns+Ojn3Be8nkXtXXzrG46xaQLdoqbJnMz+dUPItEqdP8ABUVD5vMKgKq1CYJbMfhOqoqvI5pCBbgj\n1VVQ6jgy0gQyNJ2TqqVY644eoJCUtpQ4k5XqORt9KVAfgVhIKmWlAYAlW+1UAXw9ICptmzoxpk/e\ngUTw1Dq4+HRp6oJMVQQ+F0rSpXwba0AThRj8yDUpFP8Ag2yS2XF2qUXJy2pKpj0zVoxOJeHLXmx8\nOlCu4O571ajBPAHrZRVahKsedl3KSPQzINKsV4VxXi/COJDkLNpcqPyJVqbX1SlQ223kVUewV+0R\nJea5nCkIf1RcclWhtaT00+bSr2xWYr1/B73wlfPqZKn2ytjnNjSnUk51JOFAxHcVBe44PbtJt3EM\n6mnvxK0kZ2IGCKlWOb8NjUCtCQXlcphSt9Z6SJpSKXfhdNoCq4QlkzpDgJOT6CanRuFXuAqYZW/q\nYFqmNWsLUSVkD16ntSjVt/ArpsF8QTcMKtcaRpVOTEQU7zUqs0eGbhNyWypuEGVJHT/TsOlWh88F\nui2nl6Y3SnEe5H6CpVijPBLgNmHAI+Ynv9qIV4h4XVcJSxdICkvZT5yDPRQI2q0jKP7P75sa2HXl\nNE6SC6lSZHQheadEJ8Q8GcUhTauH29zOwcKEn6ETFOjlhXHgZtGou2D1so7Bl9taQOvzEfpWuk5Y\nrvhi0TdFlq81ThCFoIVPqU4irWU/8F8dSHHGg06lI+XVCv8A8hFOljJuOC31ujW+zCPxQpPUYq0g\nLLBLUkeVOxx1qizbC0kFsaSdlAxUR//Z\n",
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 1,
       "text": [
        "<IPython.core.display.Image at 0x7f31a53e4810>"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Above:Dust storm hits Hooker, OK, June 4, 1937. \n",
      "\n",
      "To find out how much rainfall was there during the dust bowl years, we can use the [USGS/CIDA GeoDataPortal (GDP)](http://cida.usgs.gov/climate/gdp/) which can compute statistics of a gridded field within specified shapes, such as county outlines.  Hooker is in Texas County, Oklahoma, so here we use the GDP to compute a historical time series of mean precipitation in Texas County using the PRISM dataset. We then compare to climate forecast projections to see if similar droughts are predicted to occur in the future, and what the impact of different climate scenarios might be. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.dates as mdates\n",
      "import matplotlib.pyplot as plt\n",
      "import urllib\n",
      "import os\n",
      "from IPython.core.display import HTML\n",
      "import time\n",
      "import datetime\n",
      "import pandas as pd\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from lxml import etree"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "etree.__version__"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "u'3.4.4'"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cd /usgs/data2/notebook"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "/usgs/data2/notebook\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pyGDP\n",
      "import numpy as np\n",
      "import matplotlib.dates as mdates"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pyGDP.__version__"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "'1.3.3'"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "One way to interface with the GDP is with the interactive web interface, shown below.  In this interface, you can upload a shapefile or draw on the screen to define a polygon region, then you specify the statistics and datasets you want to use via dropdown menus.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "HTML('<iframe src=http://screencast.com/t/K7KTcaFrSUc width=800 height=600></iframe>')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<iframe src=http://screencast.com/t/K7KTcaFrSUc width=800 height=600></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<IPython.core.display.HTML at 0x7f318302ef50>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here we use the python interface to the GDP, called PyGDP, which allows for scripting.  You can get the code and documentation at https://github.com/USGS-CIDA/pyGDP."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create a pyGDP object\n",
      "myGDP = pyGDP.pyGDPwebProcessing()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:owslib:http://cida.usgs.gov/gdp/process/WebProcessingService?version=1.0.0&request=GetCapabilities&service=WPS\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:owslib:<wps100:Capabilities xmlns:ows110=\"http://www.opengis.net/ows/1.1\" xmlns:wps100=\"http://www.opengis.net/wps/1.0.0\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" service=\"WPS\" updateSequence=\"1\" version=\"1.0.0\" xsi:schemaLocation=\"http://www.opengis.net/wps/1.0.0 http://schemas.opengis.net/wps/1.0.0/wpsGetCapabilities_response.xsd\" xml:lang=\"en-US\">\n",
        "  <ows110:ServiceIdentification>\n",
        "    <ows110:Title>Geo Data Portal WPS Processing</ows110:Title>\n",
        "    <ows110:Abstract>Geo Data Portal WPS Processing</ows110:Abstract>\n",
        "    <ows110:Keywords>\n",
        "      <ows110:Keyword>WPS</ows110:Keyword>\n",
        "      <ows110:Keyword>AAFC</ows110:Keyword>\n",
        "      <ows110:Keyword>geospatial</ows110:Keyword>\n",
        "      <ows110:Keyword>geoprocessing</ows110:Keyword>\n",
        "    </ows110:Keywords>\n",
        "    <ows110:ServiceType>WPS</ows110:ServiceType>\n",
        "    <ows110:ServiceTypeVersion>1.0.0</ows110:ServiceTypeVersion>\n",
        "    <ows110:ServiceTypeVersion>0.4.0</ows110:ServiceTypeVersion>\n",
        "    <ows110:Fees>NONE</ows110:Fees>\n",
        "    <ows110:AccessConstraints>NONE</ows110:AccessConstraints>\n",
        "  </ows110:ServiceIdentification>\n",
        "  <ows110:ServiceProvider>\n",
        "    <ows110:ProviderName>CIDA</ows110:ProviderName>\n",
        "    <ows110:ProviderSite xlink:href=\"http://cida.usgs.gov/\" />\n",
        "    <ows110:ServiceContact>\n",
        "      <ows110:IndividualName>Tom Kunicki</ows110:IndividualName>\n",
        "      <ows110:PositionName>Technical Lead</ows110:PositionName>\n",
        "      <ows110:ContactInfo>\n",
        "        <ows110:Phone>\n",
        "          <ows110:Voice>608-821-3851</ows110:Voice>\n",
        "          <ows110:Facsimile />\n",
        "        </ows110:Phone>\n",
        "        <ows110:Address>\n",
        "          <ows110:DeliveryPoint>8505 Research Way</ows110:DeliveryPoint>\n",
        "          <ows110:City>Middleton</ows110:City>\n",
        "          <ows110:AdministrativeArea>WI</ows110:AdministrativeArea>\n",
        "          <ows110:PostalCode>53562</ows110:PostalCode>\n",
        "          <ows110:Country>US</ows110:Country>\n",
        "          <ows110:ElectronicMailAddress>tkunicki@usgs.gov</ows110:ElectronicMailAddress>\n",
        "        </ows110:Address>\n",
        "      </ows110:ContactInfo>\n",
        "    </ows110:ServiceContact>\n",
        "  </ows110:ServiceProvider>\n",
        "  <ows110:OperationsMetadata>\n",
        "    <ows110:Operation name=\"GetCapabilities\">\n",
        "      <ows110:DCP>\n",
        "        <ows110:HTTP>\n",
        "          <ows110:Get xlink:href=\"http://cida.usgs.gov:80/gdp/process/WebProcessingService?\" />\n",
        "        </ows110:HTTP>\n",
        "      </ows110:DCP>\n",
        "    </ows110:Operation>\n",
        "    <ows110:Operation name=\"DescribeProcess\">\n",
        "      <ows110:DCP>\n",
        "        <ows110:HTTP>\n",
        "          <ows110:Get xlink:href=\"http://cida.usgs.gov:80/gdp/process/WebProcessingService?\" />\n",
        "        </ows110:HTTP>\n",
        "      </ows110:DCP>\n",
        "    </ows110:Operation>\n",
        "    <ows110:Operation name=\"Execute\">\n",
        "      <ows110:DCP>\n",
        "        <ows110:HTTP>\n",
        "          <ows110:Get xlink:href=\"http://cida.usgs.gov:80/gdp/process/WebProcessingService?\" />\n",
        "          <ows110:Post xlink:href=\"http://cida.usgs.gov:80/gdp/process/WebProcessingService\" />\n",
        "        </ows110:HTTP>\n",
        "      </ows110:DCP>\n",
        "    </ows110:Operation>\n",
        "  </ows110:OperationsMetadata>\n",
        "  <wps100:ProcessOfferings>\n",
        "    <wps100:Process wps100:processVersion=\"1.0.0\">\n",
        "      <ows110:Identifier>gov.usgs.cida.gdp.wps.algorithm.FeatureCategoricalGridCoverageAlgorithm</ows110:Identifier>\n",
        "      <ows110:Title>Categorical Coverage Fraction</ows110:Title>\n",
        "    </wps100:Process>\n",
        "    <wps100:Process wps100:processVersion=\"1.0.0\">\n",
        "      <ows110:Identifier>gov.usgs.cida.gdp.wps.algorithm.FeatureCoverageOPeNDAPIntersectionAlgorithm</ows110:Identifier>\n",
        "      <ows110:Title>OPeNDAP Subset</ows110:Title>\n",
        "    </wps100:Process>\n",
        "    <wps100:Process wps100:processVersion=\"1.0.0\">\n",
        "      <ows110:Identifier>gov.usgs.cida.gdp.wps.algorithm.FeatureGridStatisticsAlgorithm</ows110:Identifier>\n",
        "      <ows110:Title>Area Grid Statistics (unweighted)</ows110:Title>\n",
        "    </wps100:Process>\n",
        "    <wps100:Process wps100:processVersion=\"1.0.0\">\n",
        "      <ows110:Identifier>gov.usgs.cida.gdp.wps.algorithm.FeatureWeightedGridStatisticsAlgorithm</ows110:Identifier>\n",
        "      <ows110:Title>Area Grid Statistics (weighted)</ows110:Title>\n",
        "    </wps100:Process>\n",
        "    <wps100:Process wps100:processVersion=\"1.0.0\">\n",
        "      <ows110:Identifier>gov.usgs.cida.gdp.wps.algorithm.FeatureCoverageIntersectionAlgorithm</ows110:Identifier>\n",
        "      <ows110:Title>WCS Subset</ows110:Title>\n",
        "    </wps100:Process>\n",
        "  </wps100:ProcessOfferings>\n",
        "  <wps100:Languages>\n",
        "    <wps100:Default>\n",
        "      <ows110:Language>en-US</ows110:Language>\n",
        "    </wps100:Default>\n",
        "    <wps100:Supported>\n",
        "      <ows110:Language>en-US</ows110:Language>\n",
        "    </wps100:Supported>\n",
        "  </wps100:Languages>\n",
        "</wps100:Capabilities>\n"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Let's see what shapefiles are already available on the GDP server\n",
      "# this changes with time, since uploaded shapefiles are kept for a few days\n",
      "shapefiles = myGDP.getShapefiles()\n",
      "print 'Available Shapefiles:'\n",
      "for s in shapefiles:\n",
      "    print s"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Available Shapefiles:\n",
        "upload:testUpload\n",
        "derivative:CONUS_States\n",
        "derivative:NCA_Regions\n",
        "sample:simplified_HUC8s\n",
        "derivative:Level_III_Ecoregions\n",
        "upload:CIDA_TEST_\n",
        "derivative:FWS_LCC\n",
        "sample:Alaska\n",
        "sample:CONUS_Climate_Divisions\n",
        "derivative:US_Counties\n",
        "sample:nps_boundary_2013\n",
        "upload:OKCNTYD\n",
        "upload:CLP\n",
        "upload:SPF\n",
        "sample:CONUS_states\n",
        "draw:american\n",
        "sample:FWS_LCC\n",
        "sample:CSC_Boundaries\n",
        "draw:clp\n",
        "derivative:wbdhu8_alb_simp\n",
        "upload:Merced_HRU_2014\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Is our shapefile there already?\n",
      "# If not, upload it. \n",
      "OKshapeFile = 'upload:OKCNTYD'\n",
      "if not OKshapeFile in shapefiles:\n",
      "    shpfile = myGDP.uploadShapeFile('OKCNTYD.zip')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Let's check the attributes of the shapefile\n",
      "attributes = myGDP.getAttributes(OKshapeFile)\n",
      "print \"Shapefile attributes:\"\n",
      "for a in attributes:\n",
      "    print a"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "DEBUG:owslib:Making request: http://cida.usgs.gov/gdp/geoserver/wfs?maxfeatures=1&typename=upload%3AOKCNTYD&version=1.1.0&request=GetFeature&service=WFS\n"
       ]
      },
      {
       "ename": "AttributeError",
       "evalue": "'Element' object has no attribute 'nsmap'",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-12-6af4325a2652>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m# Let's check the attributes of the shapefile\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mattributes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmyGDP\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetAttributes\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mOKshapeFile\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[1;34m\"Shapefile attributes:\"\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mattributes\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m     \u001b[1;32mprint\u001b[0m \u001b[0ma\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;32m/home/usgs/python/pyGDP/pyGDP.pyc\u001b[0m in \u001b[0;36mgetAttributes\u001b[1;34m(self, shapefile)\u001b[0m\n\u001b[0;32m    122\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    123\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mgetAttributes\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mshapefile\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 124\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mshapefile_value_handle\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetAttributes\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mshapefile\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwfsUrl\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    125\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    126\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mgetValues\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mshapefile\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mattribute\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgetTuples\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'false'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlimitFeatures\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;32m/home/usgs/python/pyGDP/pygdp/shapefile_value_handle.pyc\u001b[0m in \u001b[0;36mgetAttributes\u001b[1;34m(shapefile, WFS_URL)\u001b[0m\n\u001b[0;32m     28\u001b[0m     \u001b[0mgml\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0metree\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mparse\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfeature\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     29\u001b[0m     \u001b[0mgml_root\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mgml\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetroot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 30\u001b[1;33m     \u001b[0mname_spaces\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgml_root\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnsmap\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     31\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     32\u001b[0m     \u001b[0mattributes\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mAttributeError\u001b[0m: 'Element' object has no attribute 'nsmap'"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# In this particular example, we are interested in attribute = 'DESCRIP', \n",
      "# which provides the County names for Oklahoma\n",
      "user_attribute = 'DESCRIP'\n",
      "values = myGDP.getValues(OKshapeFile, user_attribute)\n",
      "print \"Shapefile attribute values:\"\n",
      "for v in values:\n",
      "    print v"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# we want Texas County, Oklahoma, which is where Hooker is located\n",
      "user_value = 'Texas'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Let's see what gridded datasets are available for the GDP to operate on\n",
      "dataSets = myGDP.getDataSetURI()\n",
      "print \"Available gridded datasets:\"\n",
      "for d in dataSets:\n",
      "    print d[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dataSets[0][0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df = pd.DataFrame(dataSets[1:],columns=['title','abstract','urls'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df.head()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print df['title']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df.ix[20].urls"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# If you choose a DAP URL, use the \"dods:\" prefix, even\n",
      "# if the list above has a \"http:\" prefix.\n",
      "# For example:  dods://cida.usgs.gov/qa/thredds/dodsC/prism\n",
      "# Let's see what data variables are in our dataset\n",
      "dataSetURI = 'dods://cida.usgs.gov/thredds/dodsC/prism'\n",
      "dataTypes = myGDP.getDataType(dataSetURI)\n",
      "print \"Available variables:\"\n",
      "for d in dataTypes:\n",
      "    print d"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Let's see what the available time range is for our data variable\n",
      "variable = 'ppt'  # precip\n",
      "timeRange = myGDP.getTimeRange(dataSetURI, variable)\n",
      "for t in timeRange:\n",
      "    print t"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "timeBegin = '1900-01-01T00:00:00Z'\n",
      "timeEnd   = '2012-08-01T00:00:00Z'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Once we have our shapefile, attribute, value, dataset, datatype, and timerange as inputs, we can go ahead\n",
      "# and submit our request.\n",
      "name1='gdp_texas_county_prism.csv'\n",
      "if not os.path.exists(name1):\n",
      "    url_csv = myGDP.submitFeatureWeightedGridStatistics(OKshapeFile, dataSetURI, variable, \n",
      "          timeBegin, timeEnd, user_attribute, user_value, delim='COMMA', stat='MEAN' )\n",
      "    f = urllib.urlretrieve(url_csv,name1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# load historical PRISM precip\n",
      "df1=pd.read_csv(name1,skiprows=3,parse_dates=True,index_col=0,\n",
      "                names=['date','observed precip'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df1.plot(figsize=(12,2),\n",
      "    title='Average Precip for Texas County, Oklahoma, calculated via GDP using PRISM data ');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "df1 = pd.stats.moments.rolling_mean(df1,36)\n",
      "df1.plot(figsize=(12,2),\n",
      "    title='Average Precip for Texas County, Oklahoma, calculated via GDP using PRISM data ');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "HTML('<iframe src=http://www.ipcc.ch/publications_and_data/ar4/wg1/en/spmsspm-projections-of.html width=900 height=350></iframe>')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#hayhoe_URI ='dods://cida-eros-thredds1.er.usgs.gov:8082/thredds/dodsC/dcp/conus_grid.w_meta.ncml'\n",
      "dataset ='dods://cida.usgs.gov/thredds/dodsC/maurer/maurer_brekke_w_meta.ncml'\n",
      "variable = 'sresa2_gfdl-cm2-1_1_Prcp'\n",
      "timeRange = myGDP.getTimeRange(dataset, variable)\n",
      "timeRange"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# retrieve the GFDL model A2 more \"Business-as-Usual\" scenario:\n",
      "time0=time.time();\n",
      "name2='sresa2_gfdl-cm2-1_1_Prcp.csv'\n",
      "if not os.path.exists(name2):\n",
      "    variable = 'sresa2_gfdl-cm2-1_1_Prcp'\n",
      "    result2 = myGDP.submitFeatureWeightedGridStatistics(OKshapeFile, dataset, variable,\n",
      "            timeRange[0],timeRange[1],user_attribute,user_value, delim='COMMA', stat='MEAN' )\n",
      "    f = urllib.urlretrieve(result2,name2)\n",
      "    print('elapsed time=%d s' % (time.time()-time0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# now retrieve the GFDL model B1 \"Eco-Friendly\" scenario:\n",
      "time0=time.time();\n",
      "name3='sresb1_gfdl-cm2-1_1_Prcp.csv'\n",
      "if not os.path.exists(name3):\n",
      "    variable = 'sresb1_gfdl-cm2-1_1_Prcp'\n",
      "    result3 = myGDP.submitFeatureWeightedGridStatistics(OKshapeFile, dataset, variable,\n",
      "            timeRange[0],timeRange[1],user_attribute,user_value, delim='COMMA', stat='MEAN' )\n",
      "    f = urllib.urlretrieve(result3,name3)\n",
      "    print('elapsed time=%d s' % (time.time()-time0))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Load the GDP result for:  \"Business-as-Usual\" scenario:\n",
      "# load historical PRISM precip\n",
      "df2=pd.read_csv(name2,skiprows=3,parse_dates=True,index_col=0,\n",
      "                names=['date','GFDL A2'])\n",
      "\n",
      "# Load the GDP result for:   \"Eco-Friendly\" scenario:\n",
      "df3=pd.read_csv(name3,skiprows=3,parse_dates=True,index_col=0,\n",
      "                names=['date','GFDL B1'])\n",
      "# convert mm/day to mm/month (approximate):\n",
      "\n",
      "\n",
      "ts_rng = pd.date_range(start='1/1/1900',end='1/1/2100',freq='30D')\n",
      "ts = pd.DataFrame(index=ts_rng)\n",
      "df2['GFDL B1'] = df3['GFDL B1']*30.\n",
      "df2['GFDL A2'] = df2['GFDL A2']*30.\n",
      "df2 = pd.stats.moments.rolling_mean(df2,36,center=True)\n",
      "df2 = pd.concat([df2,ts],axis=1).interpolate(limit=1)\n",
      "df2['OBS'] = pd.concat([df1,ts],axis=1).interpolate(limit=1)['observed precip']\n",
      "\n",
      "# interpolate \n",
      "ax=df2.plot(figsize=(12,2),legend=False,\n",
      "    title='Average Precip for Texas County, Oklahoma, calculated via GDP using PRISM data ');\n",
      "ax.legend(loc='upper right');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now just to show that we can access more than climate model time series, let's extract precipitation data from a dry winter (1936-1937) and a normal winter (2009-2010) for Texas County and look at the spatial patterns.  \n",
      "\n",
      "We'll use the netCDF4-Python library, which allows us to open OPeNDAP datasets just as if they were local NetCDF files. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import netCDF4\n",
      "url='http://cida.usgs.gov/thredds/dodsC/prism'\n",
      "box = [-102,36.5,-100.95,37]  # Bounding box for Texas County, Oklahoma\n",
      "#box = [-104,36.,-100,39.0]  # Bounding box for larger dust bowl region"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# define a mean precipitation function, here hard-wired for the PRISM data\n",
      "def mean_precip(nc,bbox=None,start=None,stop=None):\n",
      "    lon=nc.variables['lon'][:]\n",
      "    lat=nc.variables['lat'][:]\n",
      "    tindex0=netCDF4.date2index(start,nc.variables['time'],select='nearest')\n",
      "    tindex1=netCDF4.date2index(stop,nc.variables['time'],select='nearest')\n",
      "    bi=(lon>=box[0])&(lon<=box[2])\n",
      "    bj=(lat>=box[1])&(lat<=box[3])\n",
      "    p=nc.variables['ppt'][tindex0:tindex1,bj,bi]\n",
      "    latmin=np.min(lat[bj])\n",
      "    p=np.mean(p,axis=0)\n",
      "    lon=lon[bi]\n",
      "    lat=lat[bj]\n",
      "    return p,lon,lat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nc = netCDF4.Dataset(url)\n",
      "p,lon,lat = mean_precip(nc,bbox=box,start=datetime.datetime(1936,11,1,0,0),\n",
      "                        stop=datetime.datetime(1937,4,1,0,0))\n",
      "p2,lon,lat = mean_precip(nc,bbox=box,start=datetime.datetime(1940,11,1,0,0),\n",
      "                       stop=datetime.datetime(1941,4,1,0,0))\n",
      "latmin = np.min(lat)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "raw",
     "metadata": {},
     "source": [
      "# look at March 1935, just before black sunday on April 14, 1935\n",
      "nc = netCDF4.Dataset(url)\n",
      "p,lon,lat = mean_precip(nc,bbox=box,start=datetime.datetime(1935,3,1,0,0),\n",
      "                        stop=datetime.datetime(1935,4,1,0,0))\n",
      "p2,lon,lat = mean_precip(nc,bbox=box,start=datetime.datetime(2009,3,1,0,0),\n",
      "                       stop=datetime.datetime(2009,4,1,0,0))\n",
      "latmin = np.min(lat)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cartopy.crs as ccrs\n",
      "import cartopy.feature as cfeature\n",
      "states_provinces = cfeature.NaturalEarthFeature(\n",
      "        category='cultural',\n",
      "        name='admin_1_states_provinces_lines',\n",
      "        scale='50m',\n",
      "        facecolor='none')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(12,5)) \n",
      "ax = fig.add_axes([0.1, 0.15, 0.3, 0.8],projection=ccrs.PlateCarree())\n",
      "pc = ax.pcolormesh(lon, lat, p, cmap=plt.cm.jet_r,vmin=0,vmax=40)\n",
      "plt.title('Precip in Dust Bowl Region: Winter 1936-1937')\n",
      "ax.add_feature(states_provinces,edgecolor='gray')\n",
      "ax.text(-101,36.86,'Hooker')\n",
      "ax.plot(-101,36.86,'o')\n",
      "\n",
      "cb = plt.colorbar(pc,   orientation='horizontal')\n",
      "cb.set_label('Precip (mm/month)')\n",
      "\n",
      "ax2 = fig.add_axes([0.6, 0.15, 0.3, 0.8],projection=ccrs.PlateCarree())\n",
      "pc2 = ax2.pcolormesh(lon, lat, p2, cmap=plt.cm.jet_r,vmin=0,vmax=40)\n",
      "plt.title('Precip in Dust Bowl Region: Winter 1940-1941')\n",
      "\n",
      "ax2.add_feature(states_provinces,edgecolor='gray')\n",
      "ax2.text(-101,36.86,'Hooker')\n",
      "ax2.plot(-101,36.86,'o')\n",
      "\n",
      "cb2 = plt.colorbar(pc2,  orientation='horizontal')\n",
      "cb2.set_label('Precip (mm/month)')\n",
      "\n",
      "\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "From the above patterns, we can see that it's significantly drier in the northwestern part of the county in both years.  We can also see that the *maximum* precip in 1936-1937 is less than the *minimum* precipitation in 2009-2010.  We can see just how much each part of the county was drier by doing the different plot below. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Hopefully this demo inspires other investigation of historical and projected climate data using the GDP and Python.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}