Shapefile_Polygon_ERDDAP_Plot.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example downloads and plots Antarctic satellite data for a given polygon. \n",
    "\n",
    "1. loads a shapefile of multiple polygons and pulls out one polygon\n",
    "2. extracts satellite data from and ERDDAP server for the area\n",
    "3. plots the returned data in polar projection\n",
    "\n",
    "Uses python 3 with helpful geospatial packages:\n",
    "* fiona - opens and writes files\n",
    "* shapely - handles polygons\n",
    "* pyproj - converts between projections\n",
    "* cartopy - map based plotting\n",
    "\n",
    "Other examples in this series build on this one:\n",
    "1. Automating generation of the data request url\n",
    "2. Outputing data as geosjon for leaflet\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Run in python 3\n",
    "\n",
    "import sys\n",
    "sys.executable\n",
    "\n",
    "# Import modules for importing and working with shape files\n",
    "import fiona\n",
    "from shapely.geometry import shape\n",
    "from shapely.geometry import MultiPolygon\n",
    "from functools import partial\n",
    "import pyproj\n",
    "from shapely.ops import transform\n",
    "\n",
    "# Import modules for fetching and working with the satellite data\n",
    "import urllib3 \n",
    "from netCDF4 import Dataset as netcdf_dataset\n",
    "import numpy as np\n",
    "\n",
    "# Import modules for masking and plotting the data\n",
    "from matplotlib import path\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy\n",
    "import cartopy.crs as ccrs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use fiona to open a Shapefile that contains multipolygons, read the coverage area geometries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "    \n",
    "#Read in the geometries of each polygon\n",
    "Multi = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('/home/jenn/aerdData/CCAMLR_Statistical_Areas/statistical_areasPolygon.shp')])\n",
    "\n",
    "# Uncomment the following line to see the content of Multi\n",
    "#Multi.wkt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the first polygon as an example. Uses pyproj and shapely to reproject the shapefile polygon into EPSG:4326 and make it a shapely polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Pull out a polygon in EPSG:4326 projection\n",
    "\n",
    "i = 0\n",
    "\n",
    "for polygon in Multi:\n",
    "   \n",
    "    # work with just the first polygon\n",
    "    if i == 0:\n",
    "    \n",
    "        # The projections of the shapefile data, satellite data and our plot are all different\n",
    "        # this step transforms the shapefile data into the 4326 projection of the satellite data\n",
    "        # Requests to the satellite data server need to be in 4326\n",
    "\n",
    "        project = partial( \n",
    "            pyproj.transform,\n",
    "            pyproj.Proj(init='esri:102020'), # source coordinate system\n",
    "            pyproj.Proj(init='epsg:4326'))   # destination coordinate system\n",
    "        p4326 = transform(project, polygon)  # new shapely polygon in new projection\n",
    "    \n",
    "        i = i + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Satellite Data File Retrieved\n",
      "Query was for 3627 data points\n"
     ]
    }
   ],
   "source": [
    "   \n",
    "# Get satellite data from ERDDAP\n",
    "dataRequestUrl = 'http://coastwatch.pfeg.noaa.gov/erddap/griddap/jplMURSST41mday.nc?sst[(2017-05-16T00:00:00Z):40:(2017-05-16T00:00:00Z)][(-75.428900):40:(-60.000000)][(-105.000000):40:(-68.123990)]'\n",
    "\n",
    "# specify filename and local location for the satellite data download\n",
    "file = 'dataset.nc'\n",
    "\n",
    "http = urllib3.PoolManager()\n",
    "r = http.request('GET', dataRequestUrl, preload_content=False)\n",
    "\n",
    "# Opens the file and writes the data to it\n",
    "with open(file, 'wb') as out:\n",
    "    while True:\n",
    "        data = r.read(1024*1024)\n",
    "        if not data:\n",
    "            break\n",
    "        out.write(data)\n",
    "\n",
    "r.release_conn()\n",
    "\n",
    "try:\n",
    "    # load the data into python as a netcdf_dataset\n",
    "    data4326 = netcdf_dataset(file)\n",
    "except:\n",
    "    print('Data not returned. See examples on generating data queries for help.')\n",
    "else:\n",
    "    print ('Satellite Data File Retrieved')\n",
    "    \n",
    "# Pull out temperature data \n",
    "\n",
    "parameter = 'sst'\n",
    "\n",
    "if len(np.shape(data4326.variables[parameter])) == 3:\n",
    "    data = data4326.variables[parameter][0, :, :]\n",
    "elif len(np.shape(data4326.variables[parameter])) == 4:\n",
    "    data = data4326.variables[parameter][0, 0,:, :]   \n",
    "\n",
    "lats = data4326.variables['latitude'][:]\n",
    "lons = data4326.variables['longitude'][:]\n",
    "\n",
    "# How big is this dataset?\n",
    "print('Query was for '+ str(len(lats)*len(lons)) + ' data points')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Get the perimiter of the polygon and create a mask\n",
    "\n",
    "polyListx, polyListy = p4326.exterior.xy              # perimeter of polygon\n",
    "polyList = list(zip(list(polyListx),list(polyListy))) # formatted perimeter\n",
    "\n",
    "# Create a mask to hide data outside of the polygon\n",
    "p = path.Path(polyList)                             # path for mask\n",
    "X, Y = np.meshgrid(lons, lats)                      # create the grid\n",
    "points = np.array((X.flatten(), Y.flatten())).T     # break it down\n",
    "mask = p.contains_points(points).reshape(X.shape)   # calc and grid a mask based on polygon\n",
    "datamasked = np.ma.masked_where(mask != True, data) # create masked data array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO4AAADuCAYAAAA+7jsiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvVmsbNl53/dba+1513zu2DObpDg02aRpUoI8SJFlW46N\nOLYDw0iQhwD2i4EggO1AiA0k8ENgS4ERA0H8kOQhMZAgCYxEgWNbkKWYtiTLkiXOItnNJntgd9++\n06l5z3uvlYd9dvW5956p6pwazuX5Aw2yT1ftvatqfWt90///CWMMV7jCFS4X5LYf4ApXuMLyuDLc\nK1zhEuLKcK9whUuIK8O9whUuIa4M9wpXuIS4MtwrXOES4spwr3CFS4grw73CFS4hrgz3Cle4hLCW\nefFLL71k3nnnnXU9yxWucAXAGCNOe41YpuVRCGFOer0xhvl8TlEU9Ho9pLy8B3pZlkRRRFmWBEGA\n53kIcer3uVN4+PAh165d2/ZjLAVjDEmSEMcxjuMQhiFKqW0/1sooy5LJZILv+wRBcOrrhRCbNVyt\nNePxePFlX7ZF3iDPc6IowhhDGIa4rrvtR1oZl9FwGxhjSNOUOI6xLIswDLGspRzEnYExhul0ijGG\nbrd7om1s1HCNMYxGo8XJdNlgjCHLMuI4RkpJGIbYtr3txzo3LrPhHkaWZURRhBCCMAxxHGfbj7QS\n5vM5Wms6nc6xr9mo4U4mk8WueJnQ7OpRFOE4DkEQXNpd/Sg8LYbboCgKoihCa00QBLiue+k8u9Ns\n5ayGe+5VGkURwKUyWq01SZKQJAmu6zIYDC51PP6jAtu26fV6lGVJHMfM53PCMLxU+YdOp8NoNEIp\ndS7v9FyGm6YpWZbR7/fPc5mNoaoq4jgmyzJ8378y2EsKy7LodDqL33N/fx/f9/F9f+d/TyEEvV6P\n4XCIUmrlkGzlT1kUBfP5nF6vt/O7XZMcGI/HWJbF3t4eYRju/I98hZOhlKLdbjMYDAAYDoeLxOIu\nQ0pJr9djMplQVdVq11jlTVVVMZlMLkXJpyzLxe62t7eH7/s7v9FcYTk0CcW9vT2qqmI8HqO13vZj\nnYjGaxiPxyttNEtbnTGG8XhMp9PZ+UROmqaLZ71MMfgVVoMQgk6ng+/7DIdDiqLY9iOdiCYhuorx\nLm244/GYIAh2OiXfuMZJkjAYDJ6K0s4Vzg7P8+j1ekyn00XydFfh+z6WZTGfz5d639KGa1kWvu8v\n+7aNoaqqRdau3+/vvCt/hfXAsiwGg8GlcJ1brRZlWZIkyZnfs/SqbrVay75lY8iyjNFoRKvVunKN\nr7BwnV3X3WnXuck0x3F85vcsHaTuYmLncI90v9+/1L2tV7h4+L6PbdtL9QxvGo3xnhWX3o9sXGMh\nxJXRXuFYNK5zURRMJpOdLBkts3YvteHmec5oNCIMQ1qt1k56A1fYHQgh6Ha72LbNcDikLMttP9LK\n2O16zjEwxhBFEXme/8idssYYyrI804lRliV5np/6OiEElmX9yGx8QRBg2zbj8ZgwDHc62XocLpSP\nuwlorZlMJouumad9sZVlSVEUi3+AMxvZaDQ6Uzuq1pqyLBFCYNs2juNg2/ZTvyFqrZlOp0gpd2Yt\nbZyPuwk08Umr1bqU9MHToLV+xEirqlr0szbGtMziWpYddPj+eZ6jtcayLGzbXvyzC4v7ohFFEWma\n0u12t95UtDF20KbQfLm9Xm/rX+5FwBhDVVULI2nircZAfN/f+IknpcR13YV4wOFnTJKE2Wz2yDM6\njvNUnMoNx3c8Hl+aQ2HnT9xddGdWwTZOs3XwcS/aK9gl7EIY9lS4yo1rfFkTCPCheoMxZuPx46aI\n9Ifj8DzPF0Txy9hqejjx2e12N+5RXHpXuREMu4yu8eN6SZeBkHEeWJb1SCtsnufM5/NLqdslhKDV\nai1Kje12eyeffydXU5IkpGnKYDC4VG7X4wqFvV7vqYgBl4XjODiOs1DKnM/nl04p03Ec+v0+4/EY\nIcTOkWp2znDzPCeO40tltFpr4jgmTVM8z7tS1jiAZVl0u12qqiKKIqIoIgiCS8OJVkrR6/UYjUY7\n5/ntVIxbliXj8fjSNFU0CzLP851ckLsmFvf4BhcEwaXY4IqiYDqdboRtdumSU1prhsPhoiVtl9Go\nDVZVtYjhdslgG+ya4Ta4jKLnTc6i3++v9be+VMmpRlWj1WrttNE+Lfq+24YQYuGhZFnGeDxGKbXT\nmWjP8yjLkul0Srfb3fbjLG+4RlcIebG743Q6xXXdnSx8N2LpURQt6nu7urguG4QQeJ6H53mXIhPd\narWYTCZEUbR1vvfShps/+CHuzY9c2APssi5zmqbM5/Mf6QzxptBkopswZD6f0+l0dm6TvChd5PNi\nacPV6fTCbr6ruszGGGazGVVVXWWIN4zDoueTyQTP83ZqU78oXeTzYukVKe2L6WBqdtZd02W+0qza\nDeyyZtRF6CKf+xmWfYOwzx97NLrM3W53pwzjSrNqt7DLmlHn1UU+L5ZXeexcP9cNm3Gcu9QG2LjG\nURTR7/evssU7Bt/3F3KrywiqrRvn0UU+L5Z3ld3VhbaMMUwmk53SZb7SrLocOKwZta1T7ig0usgN\n5XFT2KifOpvNdkqX+Uqz6nKh0YxyHGenNKNardaiK2xT2JjhxnGM1pp2u72pWx6LRs51Pp/T7/d3\nsmZ4heMRBMEivlxGRHxdaDaUJEnIsmwj99yI4WZZtpAG2TaaGNsYc+UaX2LYts1gMCDLMqbT6dZd\n56ZMNJvNNuIJLG240Xf/Femd18/8RRVFwWw2o9vtbt0VzfOc4XCI7/uXWk3jCjWasoxSaidcZ6UU\n3W53I+Wr5ZNTlqa49wbZnddOfW0jBbINJYHHEUURs9mMXq+3k62VV1gdYRguXOc0Tbf6LLZt0263\n155AW7oeY4Uh0nFI77yOe/NjCOvozpGGOLDt3t7DOkKXieN7FmitT9zZq6o68RRSSj0130fjOk8m\nE/I836pH5bruovNrmbEiy2BpWl/+/u8AkO0Psfov4D336SNfOx6PF3WubaHh915mzaoGx6lCnuTJ\njMfjExdO0/XzNKk2HtaM2vbg9UbkcJlBeWun9dmdNtkHb+De+vgTp24cx0gpt2q0TRLqMvB7j8JJ\nqpDNEKvTTpSyLE/d8Q/fZzqdLlQbG1G7y6ba2GhGNUPN182fPQntdpvRaESWZRdeuVjZcKVtI5Si\nSmdYrcHi7009a29v70IecBVcFn5vg2asyGHZ08NTBdapsfy4ljJ8qNqYpukTWsq2be9Mx9tJ2AX+\nbNOyOZlMdsdwAUxVIp1HT9U4jgmCYKu79C7zew9Da02SJCRJsjhNgyDY+hyfx1UbjTGLkz9N04Xy\nx66Lv+0Cf7b5XRu5ngu77qpvNFrXpPpDpAOtNWmabvW03WV+b4OqqojjmCzL8H1/56mDjcph06ba\nPP/+/j6+7+P7/s4+/y7wZ8MwZDQanSpxtEzzxtKGO40LQs9ClCXCeVQcbT6fE4bh1nbhXeX3Nmjk\nSsuyJAiCS9tm2SiBhGFIkiQMh0Nc1yUIgp1Lbu0Cf1Ypheu6JElyZN6nLEtms9lSa2Fpw3VtyXie\no5KIcPDc4u9NxnNbLY1FUTCfz3ey5JPn+WKawa7KsqwCKSVhGBIEwSIZ1Ewx2KU4uGnU2KaCaBiG\ni+afZn1qrZnP5wu7WYZ4s4LhKmwJw1FC5OwhD3z3bZ62Db932+n/x3FYXO6yJMpWgRBi4TI3LYi7\nJqh3mD+7jc1dSonv+4scUKNyGYbhSjXnlbZFnSR0ent4N24xn8+ZzWYYY7aSvds1fu+P2viRx9Fk\nqBsPqNnQd8HLOMyf3YbyShAE3L17lziOzy2cv/S7jDEU8wjn1scXvZnw4VS9TUqM7Bq/N89z9vf3\nF/XTXZi3ui3Ytk2/36fT6ZCm6eJ72TYa/ux8Pt/ofZtmIPiwLfI83uHS7ywmE1TQx+rdqv+9KLAs\ni1u3bmHbNsPhcBHPrRu7xO9teqH7/T7tdnvnkjTbQjOGZFd6iaEuEzWZ8XWjUVdpDphbt25RluW5\ntaqWNtxyNkenM/L7bwK18TQtXb7vs7e3t5hKsE5u4q7we7XWjEajhSLklcEejaaXOEmSrdPwNsWf\nTZKE/f39RZ98Uw4Kw3BRtlwVSxuu/+wz2J2Q7P3vMH79tyGLHnEHhRC02+3FFzMajS7cRdoVfm9R\nFItMYafT2bls9q5BSrnI6g6Hw60pJMJ6+bNFUbC/v09RFAwGgycakjzPoyiKc913eZVHKVGeh3fj\nOtFshPjhV8nufv+J11mWRa/XIwxDJpMJs9nsQuLfpua1bX5vHMdMp9MrmuAKaDKpTR/vtnDR/Nmm\nujGbzeh0OnQ6nSPj2KbKcJ44e+XoOCsNXq+D129TzR4c+zrHcRYu5HA4JEmSld2kw8SBbbmkzTM0\nu+mPavLpvGjWRZMb2JbrfBH82UYKqemOGgwGp5b+XNddEDxWwUqGa4whzkpCz8JUFcI6OdXfDHlq\nVPqGwyF5ni99z23ze4uiWPw42z7xnwY0rrMQYpEn2AaavvbpdPkpHU3GXAjB3t7eUt7XeU7dlQw3\nzTWOpZACqjTDHjx7tptJSafTodvtEkUR4/H4zD9WkiTYtr21emCSJEwmEzqdzk5ksZ8WNG5j08+7\n7IZ+UQiCAK31me/fHEBNi+0qzUdNCXOVz7xSHTfOSjzLkN17gAwGWJ0bZPffJnnnm5STe6dew7Is\n+v0+vu8zGo0WU9pOvOdBl8mm0dSKsyw7kwt0hdXgui79fn/RtLEN17ndbp+qj9z0K8xms0US9jxh\n26qn7tKGm2QVnqMwcYwMegQ/9pMIKSmH71I8eJPsg++d+Vqu67K3t4cQgv39/WOlNqMo2goDxRjD\naDRaDKLapXbKpxHNvKbGODaNhs54VK25UdYYDoeL0tZFbOK2baOUWjpJt7zh5hWBq1BBQDXfZ/71\nXyZ5+2tUsyHScUAsd8mmrjUYDBYqjIcD9oYquA01jYbXu00ljx81NORz4Ny1zlXQarWeaCDKsoz9\n/X2MMezt7V14qBSG4dJextIp0cCtBcaEbePdvoXOc6r4Ae71moOb7Q8xxizt70sp6Xa7CzlXpRSt\nVutDF1lXJO99m3J0B+fmx3Bvf3zZR18K8/l8salcYfPYFo+2of6laYpt2wu63TpZRYfJ9mfF0ieu\n53z48EIIlOvi9HpIx0E6DkIKiuF7y152gcYNcV2X/f39heyH0RX5B2+gfIf87hsrX/8sSNN0oRR4\nhe2gaZCIomjjU/qCIOD+/fsLocFNDDVftptq6RP3tJPU7nZI3/oq5egO3oufQ9qr7Zae55Fl2aL+\na0f3QNanvfTWZ1DN3N5tioytgqqqnnC1jpVnTSZY7e2plJwVjRe2KR6tMWZBtwuCAM/zNkZeacj2\nZ8WFdw8oz8O7dZNiOmH+zX9O+MmfQoXLa8uWZUlZluzt7VGWJXdf/3Wk61EOJ4QvvnLRjw3sLq/3\ncRxWZiyKYqHM+PhGkyTJIxlLo0vSu2+SP3wXy/VpvfRZvN4NZDpFuQHS272w4DCPdp2DxvM8Zzab\nLTw+gOFwSBAEG1sLy4Rla2n7EVLi9HrkjEnf+zbhJ/7w0teYz+cL8oJlWXTDgKIqSMObWMFNHK0v\n9AvdNV7vYRw20sMKkI1Uq1IKU5UI9ehzl2VJt9tFxxPy+2+RP3ibtudjPX+DbB5RDd9ndO9t8vkQ\nU5Y4fojt2HjtAf7tj6H83QgVGh5ts6lepCdUVdUjfPLDv31DfF9GF/k8WGY9r3WF6rzEe/7Fpd9X\nFAVa60dcB3vwHPrOa/T2AsTBmMUgCB6RAlkVu8brhUcJ+VLKhUzrUWSGKhox/9b/h733PM6Nj5Dd\nfQOdRuTuNebv/R6mzFGBh3/7FuLA3fSkJP3gPULXZ/DiLbLhiKIyGFswHd1j+ME79J77GO0XXtmJ\nkMH3/UWfepN1Pg+aNsU8z2m1Wke6qUEQsL+/v9FT96xYeZLBWZDef4j3/KvYe4e0qZIZ2Qdv4L/w\nGYR1tJGMRqNTpV601gvF+uO++LNiFcX5daFpNkmSBMdxCMPw2Ngueetr5PffRAY9hM4wVYUBrDAA\nY3jwYMqN2wPUMZub0RpxzILM05zJvX3sFz5Pa3BjZ6RYL2JCRpIki96A06SEkyShKIoL2SzOgrVP\nMjgLrMAju/MaVvcmwrKpojHR6/8aIWH2zQ8IP/VHUf6jX0jT/nVacVtKSbvdXuzCcRzTbreXdnMb\nXu+mfpjj0AjJN/q7p8maGGPI77+Fe+MG5SxC+T5WK3zkv7vaxgqOrzkeZ7QAjufQCV2q/e8TxUNm\nrRu0uv0nPBxTlaTvfQcA/8VXl/nIK6Hb7S4UG5fdrJtSYzPd/iynqOd5xHG8yCPsCtZ64hpjKEZj\nqrxESIUpM6x2C7vdpphO0ZWi9crPPPKe4XC4UpyZZRnz+XxxSp3lR2nE3LaZQS7LkjiOyfN8Kde/\nSqZE3/4y/jO3j33NcJYxaK/uieiiQOc5OsvJo5jS65ErH0cafAusoEvx4G2EJamShPbnfg7prr9Z\npaoqRqMRvV7vTOukqirm8zlVVa1EUmlkfzfB/96JE1cIgTPoUyUJCIF0PzQQKwxJPrj7SLNGmqYo\npVZKDrmui+M4C53f05T2m6TENuVcoygiTdOVlP6y919DrbmjS9o20rYhDLG7Hao0RRdTUi2YRAY/\nmeK3A6wgIH5v9dr9sjjMo21aZo9C06aYpimtVmvlRo7m1C3LcmcSlxuJuJXvox4zIqEUQimK/foH\nb77k88SZy9AHm3ttI+nwuNzNsvFjdv8tysld7M7msr5CKawwxOn16Ay6XLvepfJDEmFjjMFqd5h/\n69c2Rg6wbXthUEehodsBS9PtjsJ5ie8Xja1uH86gT/rWVzBVie7cWjRcnxcNfbAZ+tTEw821m6FW\n24hri6JgMpmsPPpT5wnp21/Hu3njxBh13ZBS0A1t4qxiPC9ot1qUkzE6HpPffxvpBghlYw2eRdrr\noWIGQbCQDmo24MMtsxfZtOE4zkK8fBcYYmuNcc8CXRSk9x8SDz7OzY98ci0nYJqmzOdzPM9bSOn4\nvr9xbm+TLT6PbGs5uUfy1u/hXr926il93hj3rMhLzSwusOMpdlWgwtqFN5WmShJU0EH6HVTQw7n5\nMkJeXJKnIQQEQbAwrE6nsxbjajaFpkFDF1ndt39MdWQV7ESMexZI20YHISp6gJRHD8l+HDqLSd78\nPfyXv3imZIjnebiuSxzH3LtX84XXNSn8KDR1YiHEuWNq6XcwGsrZDHvLmfAGjiXptRymsotUEs/7\nsIvLaI3OMnQ5Jb97j2J0B+V3sPq3sXu3ViKkHIbv+9y5c2dRVVinF2XbNlJK5nfeQN/9HiZPQCq6\nP/7n13bP47C04WbDETovsDutY+uDy6JQNmo+QqfRqW13usiIvvvrYAqSH/wu4ad/+kz3aJg+Ded3\nOBwemWE0xoDRF3YqlGW5aO64CDqYdHzs6y+iZ3cu4OkuDkoKeqFNlFaMo4JOYKOkqMUFfR8FWEFA\nGceY9CHxG+8g3QCdzFHtAe4zn8Lu3Vzqnk0loZmesAllEs/zmL7zLqFnI1s+hu2ooSxtuNV8hnPr\nY5ST+5SzOXanjTxHcd4Yg0bS6rWZfeNXkF6ICrqU0weYIqP9hT+DqQpAUIzukL//Gsr3APuREZ9n\nQZ7ni9jncCzUqMobY5h/85+js/hCShtVVS3E7S7KdTPGUDx4B6e/G6ftYQghaPkWWVExiQp6LRv5\nWELSPmBcqTDEFAWy36ZKU5If/A7l3vN4L37+1LV0eLpd01feyL2uu9bqOA70nqfafw3hC7TZjkrl\n0oargoBi/4dIJ8Da+yjFw7dhNsfp9+rSwZIoKoNtSeygjRWG6KLAFDFGlwjbZfaNX1mcftK2cG9c\no4pjiumM1keX64Gez+cLV+qwQHeT4PAdG53OUX5AOb2Pc/2lpT9Pg8O9zxcZbxX77yIwtWjBjsK1\nFcbANCrohvaRhigtCw7ifCsMUb5Peu9d7P6zWN0bR173cLfc49PtGjL6umutUkpk2EekA8pojHPj\nbHprF42lDdcYBVRUyRTvhVex+7eYf/vLsGIZIC80jlUnpISUKNcF16WYRaAr/GduP/HDi06HMkrA\nnF0L97gase/7eJ5Xy5JMpgi7hchTdL76qIxGkfKie59NVZK+8w2c/u4rTHqOotSGeVrS9k/fuISU\nWK2Q/P5bTxhu07cdRdGi3Pf4599krdXzPNQLX9iqaODSKdzWq38S/6NfwnvhVVTnOtHrv1XHuiue\nKnn5oeE+8mBeTZ43x6i9CyXR+dEaVY/jtBpxozTY7/dJi5JJUmK81V3R2Wy2IAVcJNIffhPp2KhL\nIsDe8iy0hjg7XbHfGIMpK3hsPzosZ3TUVIBH7rehWqvjOFtTo2yw/CQDIbB7t3Bvfayei/ryHwTh\nkn5wl2IyxRySWzWnqMNrXZ/SUj75Q1i+B0KR3rtPeajIboyhmExA2se6VI+jkSE5Lf5RStEyCb6t\niHFXmj4YRRFa6wslLBhdEb/5FYrhezi97Y5dWRadwCLLNVlxsgyvznOqNMN77jPAh/mBKIoWQ8NO\nKxU6jnMukfGzwrbtjatyPI5zF03t3i1an/ljhJ/6abDbJHfv1T2uRUHy/vt1u+MxyEuNYx/9CNJ1\ncXr1qZc/3EcffFH5aIwxNuGnfgpxBmG6ZTuyVOcaNiX9fn8xfTCO4zN1BGVZtpae1uQHv0c1u4d3\n4/qClndZIETdqBGlJWV1/Caoswyr/wzCDRZTAXzfp9/vL+X6nkVi9bwQQmBZ1laN98K6HVTYI/jY\nT+A9/xmy/SH5eILVvUk+nhx78ualxj3CTYYDPaswRLoHMaIxFNMZOisIP/lHkM7Z3NBmiPBZGzuE\nrBdJcf/NBUunqqpTpy80WeqLJnqX0weU03u4e4NLZ7QNpBR0AptJVFDpozdAU1bkwn1kKsAqDTJN\nrXXdM4lc192qu3zhbUrOzY9idW5iNAQ/9pOo1jWy/SH6iFi1KDWWOnqR66IgvXcfnSTY3Q75eIIu\nIHzl30Gos8XTjYbQMtxN6YYoPyC/8xqzr/4T0h/824Xw9XHTF5oGi3VI3ug8ORDh2y0i97KwlKTl\nW0zjJ0+potRMckOezBkMBitNBTiMRmJ1nXAcZ6sDyy48/SaEwP/ol6AqEFIRfPwnyN5/jfTuG1iB\nj3RdpONQIbBUvRgPd88YY8hHY6poDsYgHZdiHuE+92ncWx9f6gdtMozLGJNz62PMH76De20PjCF7\n+AHw4fSFLMsWZO5Wq4UQYnGqryObWSVThLrcRtvAtRVZUce7rq2otCFKSyptCIQmvPX8hWx8lmXV\n/QEXLG90GEqpOqF2zs6vVbEezSkh4KB/U0iF9/wrONdfIn/4DtXsIcX4IXGhUZYiGRZ1Yf6A6VLM\nIkxZ1OUlZWENXsB95hNIZ/lMapZl/Cdf/QUimTHLPH775/7uqe9RfhvvhVcp7r2O3e2AEJiyQFj1\nKX+YPri/v4/v+yRJwt7eelQTq8k9lLe7NdtlEboW4yinKDVZoQk9C89RZKmknD5Ahf0L4fQ2md91\najLbtk2e51uZZ7U0yeCiaFv79z6g5Vo4rR5VPCa78xoA1t6LoEuk7WL1b58pAXUc/tNf/3neqEq0\nMMwyj0nsEQ993vnLP3/i+0yZM/3aP8O7eYNiMkV1njlS3UFrzf3796mqihs3blx4Y7tO58y+9WtH\n1rIffQ5DUWmKsv7f5jeaxnXrIYAUom50sSS2ElurA2dFxYNJSuBa9FvO4jl0WVJMplRxTPCJP4Td\nu3Wu++R5Tpqma+1dbpKRF3mPnSYZGGMwysbtXasfonMdq3P9Qu+hteamO2fflIxzn6SwSKb17vvS\n//zfQKwIfyj59i/+tSfeKywH5+ZHKcbvY/c6ZPfeRPptrM51qvmQ4sHbNTPE7yFUi0F/wHQ6wbJs\nWq3WhbXd6TJfxLZGa8r5HKvVQkiJNoY0q0iLChDYlsBWEt+te4QbNOygstKUlSHNK+ZVbdi+o/Ac\nuREjLivNPCkRQnCz5zNLHo11pWXh7g3ItFm5mecwbNte+/whx3HWnsE+Dlsx3DzP166mmOc5P937\nKj/r3efL80/zm3yULLdJ7wcYWy/q/F/8y/8tqgAj4N/+w7++eL/37KeYPXgHnaa416+R3/ku2Q+/\nhbAtrMBHuS7T2Qirekg5+QFWUWAGL7Af3yZod86dYAFQYR8V9Or6eFmi84KyKCj9Dnmp8R1Fr+U8\n0g98HCwlsdSHkygqbUiyiuEsrw3YVWe6zrLQ2hBlJUVpaPsW9kEVwbEUaV7hu48tQSUxxepdaw2E\nECil1tpJJYRASrkVPaqtuMrNMK11xgbfeuOL+PYYy4p4vejzfzz8CYZZyDffeRbx0EUHFTKR9L4j\nsRKDMKBSgywNv/l//+cAVPGE6LXfAKNRrovVadc9tkBVGSZxQb9V9+LqsqSczSjmEbnbpWrfonf7\nJfzg7CLX5WyfYv9dVNDDvvY8Qip0njD//X+BcALMs68yeeN3CT1Ja+/0tsez8HGbk7se5mbhuxez\nAI2pT/c4qwg8C89+9GTXxjCa5QzaziN/z8cTZOsW/gufPfczNOoY6xzaFkXRQnnlInBWV3krhvvw\n4cMTtYIuAl/7zk8Ttr/P17Jn+Cf7n2OYBbw/6zLebxG85qBykAV4I4NKDdoRCF3/f1UYjASZab78\nyz+PTmfk99+mePAW7rU9pOMwjQtcW+Lajy50YwxVHJNHMVEGzgufo3vt5pni3/m3fg0hSnSlQdiL\nenVZ5EymdRulvPca5OMzdVAtQ6Q3xjCNS4SAtm+d67fJC808LXEsSeAdf5LHaYkBQu/DE7GMIso4\np/25n1v5/otrHbCI+v3+ua+1qXvsbIxbVRVSrjeuKsuSt6o25Dd4ULaZFh7DNKDlZqjrmv2kT+/3\nFRioHACBSg1GgbYEqjDIvN6gfvZP/iLO/Qjevcs//t5/RfrWVyAMKbBpH9GLLITACkOsMMSezUjf\n/zpT9ePPuZfWAAAgAElEQVTYnneqxpUxGunYuP0W5XTG/Pf/Bc7H/zDzrKyld/IZ8egO3q2ztXou\ng6bDKckqRvOCTmAtynVnRVUZZmkdu3YDG3VMjb6BYwr2pxmuClGWRRVF5JMp3oufX/lzHIZlWYuZ\nSutab5u4x1HYeIEwy7K1p8/zPCdTOalxSI3NS+E+rirJSwuJofPMlPGrJabZthonQoCx6nhXaIPQ\noKK6ceSLP3mD9O2vYoD5OEKN90nefZf4vfdJ7z2oO8Qea8yw221sW9ImXrRPPj579TD8l79IOZ2B\n1rWqom3z8Hu/R7/fx3Vdsg/ewOq01tpB5buKdmAxiQqK8mx92toY5knJJC4IHIte6BxrtMYYdFlS\nZRn5gwe4ts/0gwck79+hTEtar/wM7o2XLuzzbIIQsI3e5Y2fuHmer33m7K+++Rdx3TGIgoGKmFgB\n1/2IN1KfOHXI5w7OAwtZgCzByNpQqQRGQulJZK6RpQEBYhohckWW5LSfvUkcV7S8Oq7VaYbOC4QT\nkt67j3vj+iIOBlCuQ37v+4TXnv+QPjgcHjl9wWoNsG9+lGz4LrLXJxU2bZmglKKKRpTju/i3z1cm\nOQtsJemFDuMop9dyHslSH4YxhrTQxGlJ4KpFvH8cjDFkDx4uGF+qvUf3o19kPpvS7nbO3Ma6DBrD\nXedh0bQ/bnJ8zcZP3Kqq1sqXNMbwETXjJfcuw6rFpAq47Yz5Qu8dHFWfiGJmYU9AVJC3BKUnSK4p\ntFUbMUDlSWRRIeMCfI/f+Vfv83/9L98hfjikmkeYqsJutfBuXMe9NqCKx9jXXqJ8jFZmtdtgCrI7\nry/og71ejyRJGI1GT4zB9J5/hUpLJpOYjq+QByIC+f57WGGwsX5lpQTtg/5ifYSHUJSa0bygrDT9\nloPvniEuPtCfan/h36PzxX+f1is/g+P5aKEQK45jPQ2bOHGVUk+0wa4bGzfcdccCRVEwsDOeUSUS\nw6zyiLVDXLmUWlIVCpkLrBSEBu1AfBu0DVZisCONOy7xHmRY0wwZZ5CkSN/jl/7PN/E+8ydov/RZ\nvvKrr5PevVd/HilBV6igi3mMASOEwGqFFPvvLv6mlKLX6y0UJw/TB42ByOoQiBJZlUjHp0qmlMP3\nKTdcM3QsSeAqplGxcO8rbZhEBVFa0gks2r59JC3zSBzoT82+9k+Zff2XmX/7X5C++x3M+H3ie2+t\nRZO5yaes07AOEkpru/5R2LrK40UjyzI6doorPQqjUKI2iMIoul7CPHHJrxXkY4fkOggDwQfQ+qBC\nWwIEaFeiS4XMFWgLMT0wmLLkr3zubzD+YMZLH2nzyc9dxxeC5O5dikLzD/7yP+A/+7t//IlnEpaF\nydMnNi3HcR6Rz/EsQWUEndsvY374FfIkxv/4TxJ9519hBR7OObuJVoHnKIpKk+QVxtSdT6FnPZFN\nPwuEEDh7g9pVNgZTVVSz95ClYfbwbSwB7s2XL/wzNKfuuhQrrgz3ApDnOW0nQQpJrF1SbePKAltU\nvPuwT5UriBWlDyqv3eXgQW3cstB1b7IAURm0o5AcGF6WofMcx7Mo4py3vjciiXLC6ZSqMvydv/Hr\ntLsOHEFhFErVG0I6e2LIWVMDdF2XD37rl5jkguc//mnaf+BPL16TVAVW+/rWGEJKCu6NUm503Ufa\nFFeBEOIRtRTl+1jaMBQR2Z3v4lx/8UJ1l4GFNO/TZLhPB+3kED64/+cY6jkPq5gfZrVwdWEUXx8/\nR1UodK4wjkZUYE9AJSxa7GRuEFX9j1EC7SrkPIU8B8vCCV2yqI6XykLziz//G/zWP/kO/8//+l1+\n/yv3+Nu//AtUafpEdrl2l1vEb/z2sc+tlCJ4+QvsuYbJO9/hwRtfQ2tdd+fY/kJIYJMoKs1onlNW\nhr2OWxvdGsIcKUUtfqck0Xf+Jcm736YY3bkwY2i6m9YFKeXSSinnxUZP3E3UugwQG8P/O3+Fa/a8\njnErl5v+jB+EGfGkBZZBFrWbrBLI2xKVGwgtZGFwJiXaltiTDO05yF4XGfi0Q8jGcxASjOZ7v7/P\nL/78bwDwq/ofYaoS5bepsgzrsU4a6dhUJ2gvlWVJcu9tWhb4UmLfeomH77+DYwqUtKjiuBbS2wD0\ngchbpQ0tr25TNMYwnOVra410LAlOB1Xm6Nn7pA/fQnodgk/84XOfwI307tOEjZ64mzDc2UFxVqGR\nQjMufO5lHWaFx0f39kEBpSB9pqJywSjIu4LKrUtBsjAYKZCFxiiJMAaUBN/Dv90nUz6qFS5ofg3+\ntPeXmH/rVzE6P1KmtkpS7L0XKCf3ye5+v6YuHkIcx/RffhXnmU/ivfAZyu//G/wH36bcf5PRLCIa\nz07V8DovjDHEaclonuNYkn7LWfQWCyHwnbq/eB1wbElRsRgs5t64jikiotd+kyo+H1lgG67surHR\nE3edxGaoF96DqsWDvFbEHxUhs9JDG8G8cBim9SmoYkkVaKIXK9rfV9gxaCWQZd3uKCuDEQKZFhgh\nQAhEnCBSgdZg8uIJd/infu4lTJHhP/vM0c9WVRQP3yG/9wZCSqp4TPDyFxf/PcsyXJFSRiPyuw9x\nBn2sIMAHOtowj3NGUUHbtxfGdJHIioooLXFs9UT/cAPXVsySgmANB7+tBPNDGfkmkVXOZkTf/jLh\nZ/4Yyt/cdMJdx1PlKmutcWWGErWkSKYttBGU+qAWWinwK6oKZCZxxgLtQFWBXZgP67gGVFKClIiD\n2FIoia4MZjp7Qs9ZCPirf/NLmKrEaH1kEkl5Lroosa/dqFUqjalfXxUUWUp59w2yYh+7Fdb820PX\nkFLQabkUB9Q4JQWhZx3bGLEMDtPtuuHxzRZQ13YrvR7Vh6OuJ4TA7nRqAsf47k4bbnOqb6rt8aky\nXGMM16wZrjXjTtkn0xZSGDQCKQyhkxO3UwpPwdsh1UHNX2ZQBAJVQLKnUJkkqAyi1KAUch4hkxid\nHjptDxmvMdD+4p9l9tV/is4y1BHZS+uQyqTVapGP7jP9vX8MUpJqgWU7tYrjCR6JrWr3NSsqxvMc\nz1EErlrpOz1Mt2v51pHa1kfBsSRFaXDszfXlCqnQ6XZ4r2fFU224TZZ0XTDGcEsZWkowc+4yDgO+\nG91mmntIYai0pBukTBOXONQ4I4kwoF2Q+qD1sarbIK0oR0ZZHePqCmFLdFkbrDlCqTC//3Y93OoM\nJQdp23g3rtcxqxDkBy6wOOMJ6toKx5LEWcVoni9VV63pdpo4Kwlci1ZrOcKHY0uysjpWVncdsFoh\n6b0f4r30B1ZeP+s2LCHE2kPBw9j4ibvuGPe7pYPKAyojybSNLTQ3/RnfHd1kFNVGVVYKlEE7dfeU\nzOt/7NhgxwYrqhClRhRlXZfNcgQF2tgIx0FUFeZQiPvql26S33kNp7+cnnJzulbanMqkeeK9Qiz0\nmqK0JMkqWv7JjJ68rN1ixxL022cj4D8OS0mSbD0JKoFAG/PEcwmloKoO5jqt5i43hrUuwvumM9dP\nlaustWZU+Ziywzv5Nd5O93iQtnFUiWcVeI5NktkUiVUTCpQBBLKEIgSoM8tGKPz36u4ekrRWm3Tt\nOjFVlE+cuK9+6SbKc8902h6F83wl6kCzuCg107jEVrVBH25DrHTN3jHGrETXe+RZ+ZBMddEQ4qCk\nLh7/u8Du95l/41cA6Pz4n1+6RLRuw9p05nrjrvI6JT6MMcyNjdYe72c9PkjqE7Asa9ZGy83QRlBp\nSR5bmFLVbCAfrLg+cf0HJUKDsRVimGKyDIoSYdno1BwoUH4Y3/7Z/+gT/Ln/+FOoYHsDoABsSzJo\nO6R57T77rqrrsUlJXmpCT63Upvg4Fsa1Bny4+J/cyUxZItwA59bHV6rrrtuwmhN9U9h4HXfdrrJ1\n0N54w5nhqYJ54TAvHCypGcc+o7ttijsBGIG5loEAZwrOrK7hyrJWv9C+g+m2Eb6PCHyk56C1eeK0\n/St/5y8Btazsqs981EJdFZ6j6Lcdkqzi3jil1Jp+y74Qo4Um+7seAzhuUzBaU0Zz2p/943i3P77i\ntddruFeu8jmgteamNUJaMwqjuO+28VTBG5PrTGKf6YMQNbUwlkHkAp1YFB2NPZcIDUUoMcqm+9q8\n3tJsBZYCIVB7HfTkycym/8Jnca6/yPybv7bS5zPmfK7y4yhKzSypZWOud10EgslB8mvZOHrTkIIj\nKYQL/1muvlzXfSJu2lV+qjqnjDGMtc+wauGKgs+FNZVunrnMZx4iVahMoBKJPZW4dy2csSTv1NQ+\nowAD8XNBTTZQAt1rY7oh0lEYbZ7I/P4J+RepojHS81b6bPUpc/4fvNKGaVwwb+h2gY2lJN3QJvAs\nJnHBPDmaW7sMLtpDOIx68R/5HzjvKb+JGHeTrvJTVw56uxigdchNa0pPRXwqvEvv2YTvtG6R3LL5\n4fvXcN51cIc1B1doECUgIfxAY6W16kXZslGxQBQaYQw6To9drnb/GdJ3vkGVpkvPrq2/j9VLFcYY\n4qw6kW5Xty/adfw7ywlca2U95UqbC2n8OAraGKwjnsmUJdL2zsWO2oSr/LgowjrxVMW49T0k2kgm\nVcCwanHNmqGN4O6sTdvO8FoZ2c2S+Uc02aCu2+oDxRFtCzDUkjWAqGqjNUJQGIXT9WuCwQE+8mN9\nfv4X/gimzHGf+zTlPD7qkU6FkuLYKXYnISvqRBRAv+WcGMcKIfAPpgeUulavOKum1GGUlTl2UNt5\nUVXm6Iz3wRiYXcZTnVVet6tsWRYhkIi6zjiuAiZlwJvzemLCd969jakkQhmMqzG5oAxrYzWqlrEB\niR3VKo9l6NTyNYUmcVw6qkDY1mIR/c2/90cZXKsVKooH79RB2gqf1VLiwCDO9jnLqo5jlTy9TfFx\nSClo+/YjrY4t/+ztk6U22Osy3GNO82I6w77+4rmuXZblWucIPdWG2yjLX/SMnQaWZTHMeyQIotKl\nZ8e8He+hjcBRFUIYdCkQqURogT0VyIwFxe9ALAOh68xy5da6U2XLxnkQo67bj2g+7T9IGFwPye58\nD0yOM+hjjKGYTimnU5Qf1FP/TntuJchLjcfJlrug21V1m+J5yAaWkvQO2icnUY5rn619siw1nn3x\ny+bEjU4IpHe+PuWyLNdailz39R/HRl1ly7LWGgcopbitEp5zhnSshAd5i5vejE9173GrNUPZGnS9\nOIQGlYGVgDMHO6rdZpUZtBK10Za1RGvl1T+IyYtHklOdrluTB+IxVuAjhCAbjXntt9/if/8fvnnm\ndLFjSfJCH7tj13HsIbpd27kwhpBrq1rVAhjO8hNpe5U2aMO5GjiOw2ku+HnVP9bdQ7DOUSdH4akz\n3I/aQz7rvceP+Xe55U6JK5vKCFp2xs3+lGAvxngaIyF+VjP5ZMXws5r8YEMvfbGYalCGivS6W0vZ\nZAWFsHBevI3q9UBIun2XIk4QfpsqSchGY5L9Cf/7//gtqupoltBREELUPcDFkzFnHccWaA2DtrOY\n/dPAVBXlfE5xjgFXQggCr45/87JWvSiqJ58lzSt8Zz2Lvzwp6aU1plp93aw7KQrrVy99HBt1lS3L\nIk3PP9DpJLRkgisyPFFw2x7zTjIg1xY33BnT3OOhaWGFBdXB6WtKgYgV6XWoPPDvCxSGypVUjsBK\nDXZSkd9sU2QZtiPJD3bu//6//h26A49PvLLHT/2pl/jNX/shv/QPv8M7P5gwn2RUcYxuhUcS6x+H\n7yhmSbkwzEfpdvaRi9poTXr/AcJy0ckMoRQqCFZepLJpnzyCPtjMAhq016MdXFX6WC/CmFqva9Uc\nySZOw027yluJcdeJwIqxjY0nC74xfZ67SYeWlVPZgjf39ygKhakElltSphZCNj3LNZIbkCcCOxK4\nM407LGqGEJBHOe61AJRESMFvf/k9AP71P/8h/9Pf+wpJXH+2X9X/CID4B79LFY+R3dPJB437mZcV\nWaHPRLczVYXRmuDlP0h+702Kyb1a7/mc81ob+mCa1/RB31FIWbdVruvkKiuD7x79Wd1Bj+zhHbL3\nXLznX1n+2ms23CbE2eQIko0abvPB1pldfvnZf8m3P/hDfNK9T9p5nQ+SH+fdWQ8hDHlm1yqPlaCK\nBSIXWJFEZgJZ1bQ+WYKs6hi48AV2qLCnBbI0FLmm5ShwnIOyUB0P/tLkfzvyWar5ELt9th5mYwxK\nCO7sp9zseWei20nbxm6FFA/fJfjYj1POHhJ/77ew2u26IeCcAmmeo3BtSZSW3B9l3OitLytbacNx\nnrJQCmdvQHbv+0i/g3Pt+aWuXZblWicZbNpNhi3IszZDktb1QS3LwtMKVxg8WdSqF0BZSV64sU9S\nODUf94MWohKIUtRZZX3Aw43BntdjNzG1TGuD0nNxjOafvfP3TzUqnSfoLEb0T8+GNnQ72xL02zZw\n9o2tFmSvP6Nq7SFsn2I0opzPyXNBbjKsMKxVFFdA0810reNRlIZxkZ9KH1wWWptTFSSlZeFe2yP5\nwe+igi4qOLtXUZblWsfebDoxBVsy3HV+UMuySKsAoy3GVYCjKtLSouOldJyMKK93XtXPMfddqsBQ\neWDPxWKOUN4ReENTn8KWACn41d/8W0A92/css2hMVYKu0GmKPKR+cRhH0e2MMQfZY3Wm3uLDpAch\nBFbnOvnd7y/+Vs5mVEmCe/36mWLtx5EXmkobeq36vQv6oCUIXevsUwxOQJpXeM5JUwwNOk0pZnNU\n2EN6yxnh05ZRhi0a7rqglOLLs88wL6+jjeC2P2E/DXBVSaklUeaQxg66lEgjQH9Yv0VDGdTG+82/\n/9ePvL7v+8zn8xMN1xhD8oPfRfoB6oid3hhDlFXkR7Qp1g0RNtO4oHfKEK0qTdFpRvjJTy7+5j7z\nCezBs6jWHsVwSGdvj/SH36KY3sEZ9JcKUbQ2zJKCXuvD09pu2ieLOvscuArPWU0+p0FaVPTCoz2C\nKsvIR2OEtHFufwrnxkeW/Azrzyiv2xU/Chs3XKXU2jPLw8JnLLpUAiojkBi0Ebw37RJNPJygoBRg\npEGlElnB6//lXzvTtW3bRmtNVVXH7+K6pIpGCKkwRYE45Kamea2m6Lvq2KkAjiUpHckkKuiGRxuv\nMYZ8PMH7yBcQ1ofXl47/yNQ7IQTuM58gGn9APhrj9HtnWsjNjKCjuqoaqVbXlsRpXa5aRrfqMPJS\nYyl57MldTGa4z3wa5+bLq/VWbyD+/JGKcdeJv/WZv80v/OCvYqycYR5iq4pZ5uFaFbZfIqThzf/w\nb618fd/3SZKE1jEusFA23Z/4D4jf/Cpl9ADHcSgO4lilBP2Wc6qLGbgWQlSM57XxHn59Iw1rihy7\nf7Qc7GFI26X1mZ8l+u6vn8l4a6PNafn2icYoD9olm4HWcQZtbzn6YJJV+O7RG6DRGp2nODdeWvnU\n3ERGedNDrWELhtv0dK57SrgwDpNSUhiJNgLPKviXP/v3LuT6nucxHA4Jw/DYz2CMQXohycN3mczu\nono92qGDvURSx3cUAhhHBb3QhrIgn0zRSYwKgrrAaTRn6aMRyiL81E8x//aXKWdz7M7RSbOy0kyi\ngk5wdv1mpQS90CEvNJO4qKf8eadPPNDaUGlz7OZgtEYo+1yTDMqyXOvc2hM9rzViK0O/mnmi60xQ\n/Rcf/+/onLOeeRyklNi2fWKSSmcx4ze/RaYFrU4Pt4ix1fLllDp+hOEsQ4736b70WaxWj+y9bwMx\npioR6mzfo1AWwY/9JPNv/RrKcx/JNB+mB3ZCe6kNpoFjS/qWTdLQBz0Lzz6+rJUWJyelLkIjpygK\ngsfGwVwktpGYgi0N/Vp3gqoZq7hOtkYQBCRJcuR/S9OU0TzBvv1jdG2D3+6is5xqxdjetRVdx6Cd\nkMjpY7wu4ad+ms4X/yzSWW4zUF4L78XPk+2PFn97nB64itE2EEIQNPTB8mT6YJprvJOoiJYFRlPO\nhys9S1VVdX38Kcsow1NquEIIHMchy7K13cO2baqqeiReL4qC4XBIlmX0+306g+sIy6KKJzi3P0G+\nP1p5M1G2TaBjWr7HbDZjPB6TlccTE06Cc+15qiInyUrG85w013RDh9A7w1T5M0LKeqJ9J7CI0pJJ\nVDzCOc5LjZLixFhfCIHyPcrx3ZWeIUmStY3WbHBluBeMJoG07nvEcYzWmul0ynQ6pd1u0+12UUph\n92/T/dKfQ0cjqvlDZNChilYTlRNKYYUB+fd+g16nhe95pJOHPHxwj8lkQpqmixPmKBhj6omAScJ4\nMmWGR5HVzRTH9UJfBBr6oOdIJlFOlJaL4WLBMUmpw9CVfmKm8FlgjCFN07VycGE7GWXYUoy7icyy\nbdsYY9aaPPA8j7t37xJFEd1ul/ZBq+GR0BXeC58lfv23MNocmxw6DkII7G6XYjwh+d6/wdp7HvHW\nVwmuv4Tde5U0TUmSBK31IvE3m80eeb+UEsuyaLVauHt7CBOvhaJ3FA5PX/hgmOJY8gRSgaGcz9FZ\ngc4yhLP8qZnnObZtr11VdBsZZdiS4W4iswwfnojt9sUPi8qyjPl8ThAEFEWBd4JYXPjKz6CCLkJZ\ntD7zs8y++SuowEcuuVM3rmMxT9FxLVBn4jGO4zyROW02rb29vSeeKx/doRjfxR30l/vQ54QQAs9R\nxFmFbQlG8/rEPxxTG2OooogySvFefBWrtYd0l08uJUmy1jZH2F5GGbY4kb7JLK8TnueRZdmFJqnK\nsmQ0GpEkCd1ul16vRxiGj5xuj8Nq7y0yv9ILcZ/5JOnde+gTwoVmY9NFQRlFlHFcl5gcp6bw2T7m\nhPc3J+zjRltO7pN8/9/iXttbefLCeTCL67p0N3RoeRbzuGQaF7VmdVWRPXhImeT4L38RZ+/5lYy2\nyT2sS2nl8H224SbDlk5cYFFOWecHP5ykOm+so7UmiiLyPKfdbj9ywgVBwHA4JM/zM9UMvec+jSkL\n0rs/AATScbBaIcqvVTSM1qT37oE2IBWqNaCa7YPWSNdFSAv72gsUD97G2ntuuQ8iBAKo0mxjE+4b\nJFnN8W3qtvaBmkeaVwznOWo2pnX9BfwXXz2XJ7aJpBTUXtc6a8QnQSxzGgkhzEWdXlVVMR6P2ds7\nXZPpPCiKgvl8Tr+/mlvYJDmiKCIIAvwD43ocVVUxGo0YDAZnjquMMZgyp5rvk7799ZpNZNko3wMV\n4L34eaTXQghB8vbXKYbvYfIU+/pLBB/94qmhxsOHD7l27doTfy9n+8Rv/DZOr720nOyqqCrDJM6P\nbfMs05Tpwxnq5Z+g3W6v3PtrjGF/f3+p3+E89zkqFDkPDsLIUy+4VVdZSklRrFd283Bv8bLI85zh\ncEhRFAwGA4IT1CWUUqe6zI9DCIG0Xez+M7Q+/+/S+Ym/gFAOVZzivfAqyv8w2eW/9Hnaf+DPACz4\nqKsuGKu9h/vspygm040oExpjmCYHo0SPeWZpWfiypNfrkSQJo9FopcrDJpJSUNfqXdfdSmIKtugq\nw4dNDOuORZr7HNdb/DiqqmI2q6VSut3umd153/dJ05Qsy5Y+MRph9Nbn/9Shf3/yNZ0f/wvnFk4D\ncG58hPzOa+g8X7vLnGQVthIntlA2iTqRR/R6PfI8ZzKZYNs2rVbrzIa4iaRUc591deadBVs7cWEz\nHU5QJ6nSND31PsYY5vM5o9EI3/fp9/tLx+CdTofZbLayJ3EaofwijLa5j9V/Br1mplZWVKSFJvSO\n/x6NMeSTCQiF9OrN1XEcBoMBtm0zHA6JD5JzJ6HxrNZ9EDSewLYSU7BlwxVCLIxq3fdpNonjkCQJ\n+/v7CCHY29tbOcZSStHr9ZhMJmsPA84LFfTQR6g5XhSyoiJKK3rHUBMblLMZpoTWq3/8EUKBEALf\n9xkMBlRVtUgAHodNJaXiOF5r//NZsFXDhQ9rrdu6T9OmmOc5g8HgRMbPWWFZ1sJ4T1poW4eUaxt2\nm+QV8YHRntTWaLSmjBPc5z6DtI9OlEkpFx1pURQxHo+PzFkkSbL2TiljzJkUUNaNrRvuJpNUxpjF\nfbTWTCYTZrPZYlFcZELDsiz6/T6z2WytPdPnQTUfIs8692QJJFlJmlf0WicbrS4K0vsPsLq3sHq3\nTr1u8536vs94PF7kIaA2WsdxnvqkVIOtJqcabCpJ1el0GI/HC/e81WqtdYdWStHv9xmNanLBuk+D\nZaHjCdK52O88SkuKUp/qHkMzE+hlvOc+vZQhuK6L4zjEccz+/v7i9xwMBud9/FOx7aRUg62fuLC5\nJFVVVaRpynw+Z29vbyOGJKWk3+8Tx/HaSQ/LoIpGlNFoJQG54zBPSsrKHCu3cxhG61rEbkkNqQZC\nCMIwZDAYLEpH6+7E24WkVIOdMNx1J6nKsmQ4HJIkCbdu3aqlZDaYOGqMN0kSohXZQRcJoyvSd7+N\nfcYpC6dez9SicvpArfJUozWGYjzB6t1cqaXxMNI0pdvtsre3x2w2YzKZrM2AkyTZelKqwU4YLqwn\nSdXQ7SaTCa1Wi16vh2VZdLtdptPpRieICyHo9/uL7Og2M87F8H10Msa6APJFfkCWl6IeX3ImIboo\nQpeG4OUvneveZVkuSCS2bTMYDHBdl9FoxHw+v1APzhizUn1+Xdj+mX+Aw0mq88a6xhiSJCGOY8Iw\nfIJup5QiCALm8/lG4xUhBJ1Oh6IomE6nSzcXXBR0Mke6zrlqwoc1obvB6QJx5tAmWUymhJ/+aYS1\n+u9sjGE6ndLpdB75bT3Pw3Vdoihif3//wvIYu5KUarAzhgsfkt/PY7h5njObzRYF/OOMIggCRqPR\nVnbR5nRIkmQhOreJ+mMDnc0fmfO7DIwxJFlFkle0/Ec1oY9DGUXkoxEgkK5bi8XL8y29OI6xLOvI\nJn8hBK1Wa6GBffhUXhW7kpRqsDOuMtTZwlWTVA3d7jCp/bSTrOly2qTL3EAIQRAEDAaDRU/0ulVB\nGlQrJqUat1hTj/w8k9HGMcV4SuuVP0bnC38Gq1vLyc6/8StL339xzQMlj9N41kqpxVqYzWYrh0e7\nlD1MSD8AABLPSURBVJRqsDtPwqNJqrOeQIfpdq1Wa6nT8zAxoHuGiXrrgJSSbrdLURRMJhMcxyEI\ngrUStE2ZI+TZ+3nLShOl1SOjUs70vjimGE0IP/1TqLAH1GQJ5+bLsOJm2bjI3W73zG7r4x6O7/sn\nEkYex6Y6spbB1mh9x6GqKiaTyak1ubPS7c6Cpjd523XW5jM1ZSPXdfn/2zu3J7et5Ix/uIM3DC8a\nS5a1q81m107FVXFeUpVy9iF/dlJ5yFNSSWWzFTtx7a7XK1uS7fFIQxIkAAI4AM7JA9UQhiI5IIkb\nR/i9jDTFmcFw8OH06dP9tWmaR4t4V1tfeP0MwfMvoHQ66xZCAJwx8CiGYGzd89vrg0EGi9dT9DqG\nAkNT1q2ISQIRx+uvid9kcAWgmDrUXm+dNV4skfgBup/8A9R+ceerlHQ61tWE6tHzPujLat/bRd62\nvsYJF1gLqd/v79yT0D62qOTOMb20ZZMkCcIwTI/IKOlyiIh3CRcAeBSCvfoWyfIVIASUwQMI00Kk\nduBePQNfXMPQZXQ6GsAFZE2F0ukgvJmu/290ofTHkDsXaZLLf/4FtH4PsR9ANvro/vrvb41DORVK\n6o3H45NFlO0AGwwGO8PgIAgQRVEp9kfbOGvhUmvcZvhKbzbnHJZlFbrnCIIAQRBgOBwW9j2LggpH\nqHTSNM2053SbPQ2xT7hCiLSbhjGGMAzTon7DMCABCL7/ChBrkcaLa8T2T1D6Y/T++h+3N8MvXyOa\nfg91+BDq8MNCVyghBGazGSzLKrTCLpvM7PV6tx7cQgjM5/PC77VdRFEEXdfPV7gA0myrYRjr6Xae\nV3qZ4mKxgKIouft264BETJVClGyRJCk9UlMUBYqiYD6fYzgcpgKl15NzBr1WVVWYpnlntCHiCEJw\nyFq1WXja15b1t9k8PiTjv9VqhTiOK8sm0979rIVL4Wu328VqtTo4oXAMQohUvFWFRkVBro4kTir0\nmEwmqUDpY1POIvNAfxOylS0Tzjlc10UURWlBUFV7W9pLX15e5hLuwet/EQUSeeCcI4oizOdzPHr0\nqBIbTEmS0qoqMjc/l5tckiSoqnorpAuCoBI3iLIQQsC27TSMLRtZltMCmR9++AG9Xq/0odgEFXjk\n5eBMTNm9s9l2uwcPHqDb7VZ2vgm8FS+wDl2qijBabsM5x3w+h2EYlT98GGOYTCawLKuU8sltHFoH\nfbBwoygqpWCB9rGz2eyWbUldRRKWZUFRFCwWi1a8FUOipe1RlcRxnOZSDMNIQ+XpdFpadxctTIes\n7AcLt4wunjAMMZ1OIYTAZDK5ddhNCYnlclnoz8wDHUnZtt2KtyIot1F1GSjwdj+drX/Otg+SW0rR\nDSLHWOEcLNwih2llpwKMRiP0+/2te0rKIpftTbUNymzP5/NaSiPfJ0i0ZRsc7MLzPBiGsTWHQ/tf\nKp8sqn3wWCucg4VbhNVMtt2u1+thOBzeGSZYlgXXdUtvlt4GVWa14i0PeohbllVL61wURQjD8M79\n9Gb7oOd5J0Vjx3YdHVUmtG+o8z6EEFitVpjNZukbkHeEgyzLtYXMwDrS6PV6mM/ntTw87jNRFMG2\nbVxcXNQy0mNXi+A+TNPEZDJJj3GOjQaPrYM+SrjHWM1QB0ySJBiPx0fVFlORQF0WMKZpot/vt+It\nEGquGA6HlRwzbsN13bQa7RCofXA0GiEIgoP3v6d0HR1Vx0VdPHlS2MdOBdjFYDBIM891jDiksKbO\nvdh9gKIv3/dTZ5I6iKIIURQdPVsKeOulzRg7yCDhFH/moyvqKUm1bwK64ziwbfvoqQDbII/dukJm\n4K3LfhiGlfbR3hco+uKcYzKZ1CbabVnkUzhk+gLn/CR/5pNKHh3HgSzL72zoyRStzDLF5XIJVVVr\nN+/KFqnvyorXxb4mgzpIkiRNMFZVuL+PMu+hu/rEqfd6c39bybS+fr+ftj0BuHXOVdRUgH0/e7Va\n1b7XpKesoignJSnuMxQWUyXUeDyuXbSMMcRxXNpZcXb6Av3uFJkFQQDO+Uk/++QmA8oI0jjLU719\nDoExls6+bcJK17QVpQkrbhMjEs45ZrMZRqNRZXkSeh9UVUUYhnjw4MHWPXAlKy5ZVgZBkGaLq8wM\n6roOVVUrmT2UB/I46vf7ab31+1pxRTXnruumvk9NEC2w3uKVbQ+0CUVmvu8jjmOEYXjSvXG0cIMg\nwHQ6BQB89NFHAFCLV/BgMEAYho2aEpB1mHzfwufsWX1TwuIs9DCtIzcSBAFM08SjR49uDZs7hoND\nZVryqYaYnlpxHKdeUVU/WcmpwDTN2pNVm2SPw/bZ8ZRBlaEyle65rrvVTaJu6JRDCFFYFvkQttkj\nxXGM5XKZ7oepV7qURvrpdLpzH+t5XrrPrRrq3aQztKbBGEsdFcgErmwRly1cmn4YBAEYY9A0Db1e\nr1ErLPC2MooM6euATCG2Hf+EYQjHcWCaJgaDQfUOGLTy9fv92krXmu5gkc0LlC3iMoS7TaymaULX\n9cbsYbNU6aCxC9/3wRjbawFMW4x+v1+PdU0cx7BtuzLLj02yT9cmJUS2QXaslNwrWsRFCpcxdjZi\nJap20NjGoQ6itbo8Vm2ytY2sG2STby6Cc56uxJzzVMSnhJ2nCjeKonS1OBexEpzzdBZynXmPrOlh\nHmoVbt0hM+G6LuI4Psj1vglsilhV1Vtmb9mP+8gj3E0HSPoYRVHq/tikYVd5IAcNasesi2MWsNp9\nlZtiMk5lZ8Ph8KxuPmKbsOjf1BtM/spky0r/ns1mGI/HW4WZ/dpNW1dZlqGqd8+5bSJJksC27dRm\ntS6O3TLWLlwg36a8ClarVWp23qQjiqLICjIrUhLupqBJ6PcNWiwGg0Gtc2xPiTjzCrfUvH2n00kd\n+Ot8I6nRYT6fYzQa3bubdpcQoyiq/aFZFbTCWZZV6/YMWC8UmqaVeh2l38F1jrLM0jpY3F/qdtDI\nknWJLJPShZsdZVk3rYPF/aMJDhpE0f29+6gkZux0OuCcN6Jm1zCM1Oi6bYA/bxhjqWibUK3leR50\nXa/kAVLZZu/i4gKu69YeMgPrJoCLiwvYtp1OwGs5L3zfh+M4hTmrnAq5RFZVnVWZcOt2adxE0zSM\nRiOsVivYtt2GzmdC1qyhyn7afRzjEnkqlaZX6VytKS14iqJgNBqlnsmneuS2lAd5cdMwNsuyGnM6\nsM9IvSwq/80ty4LneY3aX9KMGBqe3IbPzcL3/Vte3HUnobLQUPCqa6ErF64sy+n+so7G+12QRy55\nBLXhc/1sepgd48VdJtSOV0dVXi27ek3TMBwO07O3Jj1BVVXFaDRCGIa3JsY16Ya572QHTFuW1aj7\ngwiCAJ7n1VbQU9smQVVVDIdDLBaLo+07ymQzfG7iNd5HKCxWVbVxYTHh+z5Wq1WtVXil1irngYrC\nt3nPNoU4juE4Ttrj24RMZh6a4PKYF7JxIROEpiSeNqG697KcRRtRq5wHyuzO53MIIRo50oPC5yAI\n2vC5YLLG4U0Niwm6zibYATfisSbLcnqm2pSjom1kJ7S14fPpUFisKEpjw2LCcRxEUdSY9tDaV1yC\nxEsrb9PcGgnKPpumCcdxsFqt0Ov1Gn3TNQnyrHJdNxVsU8NiYrlcpkPrmiBaoEHCBdaiGI1GsG0b\nQojafILykM0+0/QC0zRPtpu5r2za4FQ58eJYsv5lTWuPrD05tY1zcGvcpAzPqFOpOzmVFeu52eDU\n5Q7ZCAeMU2iCF+6xNEXEdQiXrFvDMDw7sRJ1ukOevXAJ2l+ci1vjJiRimiVcpYirEu59ECtRtzvk\nvREucL5ujZtQT3IQBOnRl2mapZ0Llync+yRWognukPdKuMD5uzVuUoWIixYu2bIEQQBFUdJrvg9/\njya4Q3LmQzG651GAkRcakj2fz++FW6Msy+h2u+h2u6mIF4vFnbapefyUD0UIkVq27rOBJbE28QiH\n3bxE8OJLSIoK8+d/A3X4KPcDpSnukMGLL3O/9myEC9xft8asiIlN8TDGbgmLvm6XsOm9IVFuE2SS\nJGn/8eaDQtf10h4URSMEh//Nf8K4vIQQAv6z30LpjdD9+HNId1x7U9whw1ffIV68yv36sxIugLS1\ni1bec6kbPhQS376zzk1B0sBk+txisYAkSe+IMjsZoemizEM8+xGQFShv9qVKpwM2nWL19b+h88u/\ng6xvD33JaK7uDrVktUDw/H9gXubf1pydcAGk+ypaee+reO8ij/Amk0lFV1MPQgis/vQfkDJZekmS\noE8miOwFnC/+CcajX8F48umt0DnrDlnXWbuIGaL5FYLvv4J+YUE+YMU/S+ECSDOY2SzgfUiStByG\nJEkwnn6G+PU373xeHw2hDvpgN98Bigbz8ScQQsB13TTRWZRoOfMR2z8hdm6gDj+EPnmy87VCcPjP\n/hvR9CUUswOt34Fy4NHT2QoXWLs1TiYTuK6L2WzW+O6SlnLQLh6C/fiHW59LwhCJ74PHCXgUAW8S\ngK7rotPpYDwev/OgT7w52M1LJM5ryGYfSn8C/YO/gCTvjujCqz9BmzxB+NM3iF4/g9rrIXj2AxJv\nBqU7gmx0IWIGuTOArHcQL1+DXX0NHrroPH585x58F2ctXABpj+y59HO2FI+kGRA8AWcMsq6DxzHY\nzRT6o19D7w0hOkN4AYP0po82u7USPEG8eIXw6mvwlQ2l24HaMcBjF+x6CvbqW/T+6jeQ9bfnuok3\nh6Qa4IGD4PkXkHQT+uVTRNd/hjoYQDFNxO41ksWPEHECCA7BBQRPIOs6FNOAcfnuMLCE5z9qPVi4\nnPNGioIcE6hVrA2f3x9kzYD59DOw7/8P+oMJwpspjCefQn/0K3ieh3AVYDAY3MoaJ4EL/9lvkbgz\nyLoBpWNC//DtEZICQO33ES8dOF/8M4zHn0DpjRDPfwS7fgb94V8idl5DUlXwwIM++RkkzQBnDIph\n3PpZgnMkvg/FNCHtyMfECcfCy+/BdrBwm57N7XQ6MAwDnue14XNB8ChEvLgGDzxo48eQzT4ECyAZ\nzTET0C9/gXj6EsHVFYzHH4MPn2A6ne4Mi6PZD0Ds7w1XJUmCdmFB6XYQzV8guvkOsiLDuLwEm76E\npKjrDPbV19Af/hKc+ZC10bvfR5ah7ql5jhKOpRfB6uW/Tw8W7mAwaHw2V5blNnwuCJFE8L76V0AW\nkBUZ7OqPEDwBJBnGhx/D/NmndV8igLXIuh9/jigK4fkMUhhuvUeFEBBRAB44kA091x5T1jQY49uC\n1IQAm9uQjS4kTUP0+jkkST54zxrFHMtVhIueBlXJ/7UHC5fGd9DK2+Te0zZ8Pg0hOLw//jtkVYL2\nxq5FDIcAgMT3ES9fAWiGcIUQ8Px17fRmWBzZ1wi+/d36dTGDSCJImgF9eHzXmdrtQjFNQAhEjoNk\nZa8faAfAYg7XjzDs61Dkw+7Jo5YgTdMa6Y28CwqX4jhOfXpb7oZdfwvBXGijt/XhkiRBkiRwxqBa\nH9R8hWvCMMR0OoUkSRiPx+9UQEU3zyEbCvTRAMbDS2gXQ0iyBPnEmmRJliEpynovbF8BQkDknI0V\nRglcP8ZF73DRAidklbPeyHWXi+VBlmVYloUoiuA4Ttog3YbPu4mX11C770YognMkKx/mzx/XdGVr\nkiRJe7b3bd0S5wb66G2Bg2QNoA76hUVesqrC/OASzF4AOb5nwBKswgTDngb5CNECJ5rFkX3Lcrk8\nG+O07BiL2WyG1WrVzgvaQeLMIG8puo8cB4r1AZT+u4mYKqAiCtu20e129yZLeczAWQApk6CUJOno\n89NdSIoCY/JuEmwTP0zgswTD/vGiBQpweSR7VcdxzmrmThs+70fEEUTMbpUSAgCPY8Suh87Tz2q5\nrs2w+K5unvD730N501lWN6swRhi9WWlPvJ5CMktZb2TOeW1NyIfShs+7SVY2ZE2/dcMnvg82m8N8\n8ilko1p3iLxhcRY2+xHs9TOYDx9WcIX78YIYUcJx0dMKeYgUlhIme1XbtgHgbMQLvA2fKftMzdRN\neEpXTeI7iJc3CF98CS2TdU2CAGxmo/vJb6Ba1flYCSHgeR6CIMjdLyt4Av/P/4V4+QrGZAK55pMP\n14/ABXDRLUa0QMElj1nxcs4bba+6DSrecF03nXlqmuZ7U8AheAL3f/8FsqZDn4zWxx0AeBSBTWfo\nfvx5JaIl7+UgCMAYS43o89707OYFEm8K8+EHhe9lD0EIAcePIQEYdNRCF4LCH0WSJKXDvFzXrdTa\nsggofBZCIAzDdJbv+yBiSVagXf4C3HudilZwvi4hfPq3UC/KO/7ZFKumaTBNE4PB4OAbPp5fQclZ\nXFEWJFpZktDvFL/iH/odn0uS9LTwq2hpaSGe53nRQWZxLS0tzaBNn7a0nCGtcFtazpBWuC0tZ0gr\n3JaWM6QVbkvLGdIKt6XlDGmF29JyhrTCbWk5Q1rhtrScIf8PU3yAP66oW0IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff31140bba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make a colormesh plot in polar stereographic\n",
    "\n",
    "thisproj = ccrs.SouthPolarStereo() #define projection\n",
    "fig = plt.figure()    \n",
    "ax1 = plt.axes(projection = thisproj, ) # set projection\n",
    "ax1.set_global()\n",
    "\n",
    "# set type of plot (colormesh), set the data to plot and\n",
    "# indicate that it's in 4326 (Plate Carree) and needs to be projected\n",
    "dataplt = ax1.pcolormesh(X,Y,datamasked,transform=ccrs.PlateCarree()) \n",
    "\n",
    "# reduce equatorward extent of plot\n",
    "ax1.set_extent([-5500000, 5500000, -5500000, 5500000], thisproj)\n",
    "# add Land\n",
    "ax1.add_feature( cartopy.feature.LAND, zorder=1, edgecolor='none', facecolor='#fae5c9')\n",
    "# add coastline\n",
    "ax1.coastlines(color='#f0c7ab')\n",
    "# add lat lon rings\n",
    "ax1.gridlines(alpha='0.3')\n",
    "\n",
    "plt.show()\n",
    "plt.clf()\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done\n"
     ]
    }
   ],
   "source": [
    "print('done')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geo3jupyter",
   "language": "python",
   "name": "geo3"
  },
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}