From a56ef2c730b787a8b9b4f7c90bd3f49b749a19c2 Mon Sep 17 00:00:00 2001 From: yutsuo Date: Mon, 24 Jun 2019 12:37:54 -0300 Subject: [PATCH] final --- Arrays.ipynb | 52 ++ Colormaps_exampe.ipynb | 1182 ++++++++++++++++++++++++++++++++++++++++ Folium-Basic.ipynb | 369 ++++++++++--- Folium-Final.ipynb | 306 +++++++++++ Folium_2.ipynb | 338 ++++++++++-- pcolormesh.ipynb | 126 +++++ 6 files changed, 2268 insertions(+), 105 deletions(-) create mode 100644 Colormaps_exampe.ipynb create mode 100644 Folium-Final.ipynb create mode 100644 pcolormesh.ipynb diff --git a/Arrays.ipynb b/Arrays.ipynb index 33cf496..e08a0cd 100644 --- a/Arrays.ipynb +++ b/Arrays.ipynb @@ -594,6 +594,58 @@ "plt.ylim(-5,10)" ] }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np \n", + "array = np.arange(12).reshape(3,4) " + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0, 1, 2, 3],\n", + " [ 4, 5, 6, 7],\n", + " [ 8, 9, 10, 11]])" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "array" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1050000" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "700*1500" + ] + }, { "cell_type": "code", "execution_count": null, diff --git a/Colormaps_exampe.ipynb b/Colormaps_exampe.ipynb new file mode 100644 index 0000000..9550bad --- /dev/null +++ b/Colormaps_exampe.ipynb @@ -0,0 +1,1182 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.8.3\n" + ] + } + ], + "source": [ + "import os\n", + "import folium\n", + "\n", + "print(folium.__version__)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Using `folium.colormap`\n", + "\n", + "**A few examples of how to use `folium.colormap` in choropleths.**\n", + "\n", + "Let's load a GeoJSON file, and try to choropleth it." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "import pandas as pd\n", + "import requests\n", + "\n", + "\n", + "# url = 'https://raw.githubusercontent.com/python-visualization/folium/master/examples/data'\n", + "# us_states = f'{url}/us-states.json'\n", + "us_states = 'https://raw.githubusercontent.com/python-visualization/folium/master/examples/data/us-states.json'\n", + "# US_Unemployment_Oct2012 = f'{url}/US_Unemployment_Oct2012.csv'\n", + "US_Unemployment_Oct2012 = 'https://raw.githubusercontent.com/python-visualization/folium/master/examples/data/US_Unemployment_Oct2012.csv'\n", + "\n", + "geo_json_data = json.loads(requests.get(us_states).text)\n", + "unemployment = pd.read_csv(US_Unemployment_Oct2012)\n", + "\n", + "unemployment_dict = unemployment.set_index('State')['Unemployment']" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Self-defined\n", + "\n", + "You can build a choropleth in using a self-defined function.\n", + "It has to output an hexadecimal color string of the form `#RRGGBB` or `#RRGGBBAA`." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "def my_color_function(feature):\n", + " \"\"\"Maps low values to green and hugh values to red.\"\"\"\n", + " if unemployment_dict[feature['id']] > 6.5:\n", + " return '#ff0000'\n", + " else:\n", + " return '#008000'" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m = folium.Map([43, -100], tiles='cartodbpositron', zoom_start=4)\n", + "\n", + "folium.GeoJson(\n", + " geo_json_data,\n", + " style_function=lambda feature: {\n", + " 'fillColor': my_color_function(feature),\n", + " 'color': 'black',\n", + " 'weight': 2,\n", + " 'dashArray': '5, 5'\n", + " }\n", + ").add_to(m)\n", + "\n", + "m.save(os.path.join('results', 'Colormaps_0.html'))\n", + "\n", + "m" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## StepColormap\n", + "\n", + "But to help you define you colormap, we've embedded `StepColormap` in `folium.colormap`.\n", + "\n", + "You can simply define the colors you want, and the `index` (*thresholds*) that correspond." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "310" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import branca.colormap as cm\n", + "\n", + "\n", + "step = cm.StepColormap(\n", + " ['green', 'yellow', 'red'],\n", + " vmin=3, vmax=10,\n", + " index=[3, 4, 8, 10],\n", + " caption='step'\n", + ")\n", + "\n", + "step" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m = folium.Map([43, -100], tiles='cartodbpositron', zoom_start=4)\n", + "\n", + "folium.GeoJson(\n", + " geo_json_data,\n", + " style_function=lambda feature: {\n", + " 'fillColor': step(unemployment_dict[feature['id']]),\n", + " 'color': 'black',\n", + " 'weight': 2,\n", + " 'dashArray': '5, 5'\n", + " }\n", + ").add_to(m)\n", + "\n", + "m.save(os.path.join('results', 'Colormaps_1.html'))\n", + "\n", + "m" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you specify no index, colors will be set uniformely." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "0.01.0" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cm.StepColormap(['r', 'y', 'g', 'c', 'b', 'm'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## LinearColormap\n", + "\n", + "But sometimes, you would prefer to have a *continuous* set of colors. This can be done by `LinearColormap`." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "310" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "linear = cm.LinearColormap(\n", + " ['green', 'yellow', 'red'],\n", + " vmin=3, vmax=10\n", + ")\n", + "\n", + "linear" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m = folium.Map([43, -100], tiles='cartodbpositron', zoom_start=4)\n", + "\n", + "folium.GeoJson(\n", + " geo_json_data,\n", + " style_function=lambda feature: {\n", + " 'fillColor': linear(unemployment_dict[feature['id']]),\n", + " 'color': 'black',\n", + " 'weight': 2,\n", + " 'dashArray': '5, 5'\n", + " }\n", + ").add_to(m)\n", + "\n", + "m.save(os.path.join('results', 'Colormaps_2.html'))\n", + "\n", + "m" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Again, you can set the `index` if you want something irregular." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "0.01.0" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cm.LinearColormap(\n", + " ['red', 'orange', 'yellow', 'green'],\n", + " index=[0, 0.1, 0.9, 1.0]\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you want to transform a linear map into a *step* one, you can use the method `to_step`." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "3.010.0" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "linear.to_step(6)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can also use more sophisticated rules to create the thresholds." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "31100" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "linear.to_step(\n", + " n=6,\n", + " data=[30.6, 50, 51, 52, 53, 54, 55, 60, 70, 100],\n", + " method='quantiles',\n", + " round_method='int'\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And the opposite is also possible with `to_linear`." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "310" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "step.to_linear()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Build-in\n", + "\n", + "For convenience, we provide a (small) set of built-in linear colormaps, in `folium.colormap.linear`." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "0.01.0" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cm.linear.OrRd_09" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can also use them to generate regular `StepColormap`." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "0.01.0" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cm.linear.PuBuGn_09.to_step(12)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Of course, you may need to scale the colormaps to your bounds. This is doable with `.scale`." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "312" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cm.linear.YlGnBu_09.scale(3, 12)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "5100" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 24, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cm.linear.RdGy_11.to_step(10).scale(5, 100)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "At last, if you want to check them all, simply ask for `linear` in the notebook." + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Spectral_070.01.0
Greens_030.01.0
RdGy_110.01.0
PuOr_090.01.0
Blues_030.01.0
Reds_030.01.0
YlGnBu_050.01.0
Blues_040.01.0
Set3_060.01.0
YlOrRd_050.01.0
Pastel2_070.01.0
BrBG_030.01.0
Reds_070.01.0
PiYG_070.01.0
BrBG_070.01.0
Reds_090.01.0
Pastel1_030.01.0
RdPu_040.01.0
Paired_030.01.0
PuOr_040.01.0
Set1_060.01.0
RdGy_090.01.0
BuPu_090.01.0
RdPu_090.01.0
Oranges_080.01.0
GnBu_070.01.0
YlOrBr_060.01.0
PRGn_110.01.0
BrBG_090.01.0
Blues_070.01.0
viridis0.01.0
PuBu_040.01.0
PuBu_080.01.0
PuOr_100.01.0
RdBu_100.01.0
YlOrBr_080.01.0
Set3_030.01.0
Greys_070.01.0
GnBu_080.01.0
RdBu_060.01.0
Reds_060.01.0
Purples_080.01.0
Purples_070.01.0
OrRd_060.01.0
BuPu_070.01.0
Accent_060.01.0
RdYlGn_090.01.0
PuBuGn_040.01.0
PiYG_030.01.0
PuOr_070.01.0
Greys_030.01.0
Spectral_060.01.0
YlGn_070.01.0
YlGn_060.01.0
Greys_040.01.0
RdBu_080.01.0
Pastel1_040.01.0
RdBu_110.01.0
PuOr_050.01.0
Spectral_110.01.0
Reds_080.01.0
PuBuGn_050.01.0
RdBu_040.01.0
PuBuGn_060.01.0
Set1_080.01.0
PuRd_090.01.0
Set2_040.01.0
RdPu_050.01.0
PuBuGn_030.01.0
RdYlGn_100.01.0
PuBu_070.01.0
YlOrRd_040.01.0
Dark2_070.01.0
Accent_050.01.0
PuRd_040.01.0
Purples_090.01.0
PiYG_110.01.0
PRGn_070.01.0
PuOr_110.01.0
Paired_110.01.0
RdYlBu_100.01.0
YlGnBu_040.01.0
Set1_050.01.0
RdYlGn_080.01.0
BuGn_040.01.0
RdBu_090.01.0
BrBG_050.01.0
BuGn_050.01.0
YlGnBu_060.01.0
RdYlBu_040.01.0
PiYG_060.01.0
YlGn_040.01.0
Greys_060.01.0
RdYlGn_040.01.0
Set2_070.01.0
RdYlGn_110.01.0
PiYG_050.01.0
PiYG_040.01.0
Paired_080.01.0
RdBu_030.01.0
RdPu_060.01.0
Greens_070.01.0
RdBu_070.01.0
Oranges_040.01.0
Accent_080.01.0
RdGy_080.01.0
BuPu_080.01.0
Set1_040.01.0
RdGy_070.01.0
YlOrBr_070.01.0
Spectral_100.01.0
RdYlBu_110.01.0
BuPu_030.01.0
Spectral_030.01.0
RdYlBu_090.01.0
Set1_070.01.0
Blues_060.01.0
Pastel1_060.01.0
RdYlGn_030.01.0
PuBu_030.01.0
RdYlBu_030.01.0
Greys_080.01.0
Blues_080.01.0
YlGnBu_090.01.0
Paired_060.01.0
Greys_090.01.0
Pastel1_080.01.0
Set2_050.01.0
YlGn_050.01.0
Spectral_080.01.0
Purples_060.01.0
YlOrBr_050.01.0
PiYG_090.01.0
Pastel2_060.01.0
BuPu_050.01.0
PuRd_080.01.0
RdYlBu_050.01.0
YlGnBu_030.01.0
Blues_050.01.0
Pastel1_070.01.0
Oranges_070.01.0
YlOrBr_090.01.0
PuOr_080.01.0
Pastel1_090.01.0
BuGn_060.01.0
Spectral_040.01.0
PRGn_040.01.0
BuGn_080.01.0
YlGnBu_070.01.0
PuBuGn_080.01.0
GnBu_050.01.0
Paired_040.01.0
PuBuGn_090.01.0
Accent_070.01.0
Oranges_090.01.0
RdGy_100.01.0
Accent_030.01.0
RdYlGn_070.01.0
Dark2_080.01.0
Pastel2_050.01.0
PuBuGn_070.01.0
Paired_100.01.0
YlOrBr_030.01.0
RdGy_060.01.0
Spectral_090.01.0
BrBG_080.01.0
YlOrRd_080.01.0
Dark2_060.01.0
Set2_030.01.0
Pastel1_050.01.0
YlGn_080.01.0
Spectral_050.01.0
Set1_030.01.0
RdYlBu_080.01.0
BuPu_060.01.0
Pastel2_040.01.0
Set3_080.01.0
Oranges_060.01.0
Oranges_050.01.0
Set3_110.01.0
PRGn_100.01.0
RdGy_040.01.0
BrBG_100.01.0
GnBu_030.01.0
RdPu_070.01.0
PuRd_030.01.0
RdYlBu_070.01.0
GnBu_040.01.0
Accent_040.01.0
PuOr_060.01.0
RdPu_030.01.0
Purples_040.01.0
Set2_060.01.0
Paired_120.01.0
PRGn_080.01.0
Set3_070.01.0
Pastel2_030.01.0
BrBG_110.01.0
Set3_050.01.0
Set3_120.01.0
Paired_090.01.0
Set1_090.01.0
OrRd_030.01.0
PRGn_060.01.0
Dark2_040.01.0
Set3_100.01.0
Greens_040.01.0
GnBu_090.01.0
RdYlGn_050.01.0
RdGy_050.01.0
Reds_040.01.0
YlGn_090.01.0
Dark2_050.01.0
Paired_070.01.0
BuGn_090.01.0
Greens_060.01.0
OrRd_050.01.0
OrRd_040.01.0
Set3_040.01.0
PuBu_090.01.0
Greens_050.01.0
Set2_080.01.0
YlOrRd_030.01.0
Greys_050.01.0
Greens_080.01.0
RdYlGn_060.01.0
YlOrRd_060.01.0
YlOrRd_090.01.0
PRGn_050.01.0
PuRd_050.01.0
BrBG_060.01.0
Pastel2_080.01.0
PuOr_030.01.0
PiYG_100.01.0
Purples_030.01.0
YlOrBr_040.01.0
Dark2_030.01.0
RdGy_030.01.0
GnBu_060.01.0
RdBu_050.01.0
RdYlBu_060.01.0
YlGn_030.01.0
RdPu_080.01.0
Purples_050.01.0
PuBu_050.01.0
Blues_090.01.0
Reds_050.01.0
BuPu_040.01.0
YlGnBu_080.01.0
PRGn_090.01.0
PuRd_060.01.0
BuGn_030.01.0
OrRd_080.01.0
Set3_090.01.0
Greens_090.01.0
PRGn_030.01.0
Oranges_030.01.0
PuRd_070.01.0
Paired_050.01.0
OrRd_070.01.0
BuGn_070.01.0
BrBG_040.01.0
OrRd_090.01.0
PiYG_080.01.0
PuBu_060.01.0
YlOrRd_070.01.0
\n", + " " + ], + "text/plain": [ + "" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "cm.linear" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Draw a `ColorMap` on a map\n", + "\n", + "By the way, a ColorMap is also a Folium `Element` that you can draw on a map." + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m = folium.Map(tiles='cartodbpositron')\n", + "\n", + "colormap = cm.linear.Set1_09.scale(0, 35).to_step(10)\n", + "colormap.caption = 'A colormap caption'\n", + "m.add_child(colormap)\n", + "\n", + "m.save(os.path.join('results', 'Colormaps_3.html'))\n", + "\n", + "m" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/Folium-Basic.ipynb b/Folium-Basic.ipynb index 26b6f8c..171d06a 100644 --- a/Folium-Basic.ipynb +++ b/Folium-Basic.ipynb @@ -22,7 +22,8 @@ "metadata": {}, "outputs": [], "source": [ - "image = '/notebooks/resources/gpm/gpm_1d.20190531.tif'" + "image = '/notebooks/resources/gpm/gpm_1d.20190530.tif'\n", + "shapefile = '/notebooks/resources/centro-oeste.geojson'" ] }, { @@ -37,9 +38,7 @@ "10.0 -180.0\n", "-60.0 -180.0\n", "-60.0 -30.0\n", - "10.0 -30.0\n", - "ext = [[10.0, -180.0], [-60.0, -180.0], [-60.0, -30.0], [10.0, -30.0]]\n", - "geo_ext = [[10.0, -180.0], [-60.0, -180.0], [-60.0, -30.0], [10.0, -30.0]]\n" + "10.0 -30.0\n" ] } ], @@ -98,7 +97,6 @@ "rows = ds.RasterYSize\n", "\n", "ext=GetExtent(gt,cols,rows)\n", - "print(\"ext = \" + str(ext))\n", "\n", "src_srs=osr.SpatialReference()\n", "src_srs.ImportFromWkt(ds.GetProjection())\n", @@ -106,28 +104,39 @@ "# tgt_srs.ImportFromEPSG(3857)\n", "tgt_srs = src_srs.CloneGeogCS()\n", "\n", - "geo_ext=ReprojectCoords(ext,src_srs,tgt_srs)\n", - "print(\"geo_ext = \" + str(geo_ext))" + "geo_ext=ReprojectCoords(ext,src_srs,tgt_srs)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "010000" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "# # Choose colormap\n", - "# cmap = matplotlib.cm.cool\n", + "# create custom colormap\n", "\n", - "# # Get the colormap colors\n", - "# my_cmap = cmap(np.arange(cmap.N))\n", - "\n", - "\n", - "# # Set alpha\n", - "# my_cmap[:,-1] = np.linspace(0, 1, cmap.N)\n", + "colormap = branca.colormap.StepColormap(\n", + " ['#2b83ba', '#64abb0', '#9dd3a7', '#c7e9ad', '#edf8b9', '#ffedaa', '#fec980', '#f99e59', '#e85b3a', '#d7191c'],\n", + " vmin=0, vmax=10000,\n", + " index=[0, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000],\n", + " caption='rainfall (mm)'\n", + ")\n", "\n", - "# # Create new colormap\n", - "# my_cmap = ListedColormap(my_cmap)" + "colormap" ] }, { @@ -136,71 +145,71 @@ "metadata": {}, "outputs": [], "source": [ - "# # Choose colormap\n", - "# cmap = matplotlib.cm.cool\n", + "# recreate custom colormap setting transparency on the first color\n", "\n", - "# min_x = 0\n", - "# max_x = 256\n", + "change = list()\n", + "alphacolor = list()\n", + "for item in colormap.colors:\n", + " change.append(list(item))\n", "\n", - "# def my_cmap(x):\n", - "# x = np.clip((x - min_x)/(max_x - min_x), 0, 1)\n", - "# c = cmap(x)\n", - "# # multiplying alpha by 10 to make a sharper transition to the colors in the lower range\n", - "# return c[:3] + (np.clip(x*10, 0, 1),) " + "change[0][3] = 0\n", + "\n", + "for item in change:\n", + " alphacolor.append(tuple(item))" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ + "# set colormap\n", "colormap = branca.colormap.StepColormap(\n", - " ['#64abb0','#9dd3a7', '#c7e9ad', '#edf8b9', '#ffedaa', '#fec980', '#f99e59', '#e85b3a', '#d7191c'],\n", - " vmin=0, vmax=10000,\n", - " index=[0, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000],\n", - " caption='step'\n", - ")\n", - "\n", - "my_cmap = colormap\n", - "\n", - "# # Get the colormap colors\n", - "# my_cmap = colormap(np.arange(colormap.N))\n", - "\n", - "\n", - "# # Set alpha\n", - "# my_cmap[:,-1] = np.linspace(0, 1, colormap.N)\n", - "\n", - "# # Create new colormap\n", - "# my_cmap = ListedColormap(my_cmap)" + " alphacolor\n", + " , vmin=0, vmax=10000\n", + " , index=[0, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000]\n", + " , caption='rainfall (mm)'\n", + ")" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# set scale\n", + "colormap_scale = colormap.to_linear().scale(0, 1000).to_step(10)\n", + "colormap_scale.caption = 'rainfall(mm)'" + ] + }, + { + "cell_type": "code", + "execution_count": 59, "metadata": { "scrolled": false }, "outputs": [ { - "ename": "IndexError", - "evalue": "tuple index out of range", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 15\u001b[0m \u001b[0;34m,\u001b[0m \u001b[0mbounds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mext\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mext\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 16\u001b[0m \u001b[0;34m,\u001b[0m \u001b[0mmercator_project\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 17\u001b[0;31m \u001b[0;34m,\u001b[0m \u001b[0mcolormap\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmy_cmap\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 18\u001b[0m )\n\u001b[1;32m 19\u001b[0m )\n", - "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/folium/raster_layers.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, image, bounds, origin, colormap, mercator_project, pixelated, name, overlay, control, show, **kwargs)\u001b[0m\n\u001b[1;32m 274\u001b[0m origin=origin)\n\u001b[1;32m 275\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 276\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0murl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mimage_to_url\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimage\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morigin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0morigin\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolormap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcolormap\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 277\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 278\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbounds\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjson\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloads\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjson\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdumps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbounds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/folium/utilities.py\u001b[0m in \u001b[0;36mimage_to_url\u001b[0;34m(image, colormap, origin)\u001b[0m\n\u001b[1;32m 109\u001b[0m \u001b[0murl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'data:image/{};base64,{}'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileformat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb64encoded\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 110\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0;34m'ndarray'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mimage\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__class__\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 111\u001b[0;31m \u001b[0mimg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mwrite_png\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimage\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morigin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0morigin\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolormap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcolormap\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 112\u001b[0m \u001b[0mb64encoded\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbase64\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mb64encode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mimg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdecode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'utf-8'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 113\u001b[0m \u001b[0murl\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'data:image/png;base64,{}'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb64encoded\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/folium/utilities.py\u001b[0m in \u001b[0;36mwrite_png\u001b[0;34m(data, origin, colormap)\u001b[0m\n\u001b[1;32m 172\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mnblayers\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 173\u001b[0m \u001b[0marr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcolormap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mravel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 174\u001b[0;31m \u001b[0mnblayers\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0marr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 175\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mnblayers\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 176\u001b[0m raise ValueError('colormap must provide colors of r'\n", - "\u001b[0;31mIndexError\u001b[0m: tuple index out of range" - ] + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 59, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ + "# draw map\n", "m = folium.Map(\n", - " location = [-22, -114]\n", - " , zoom_start = 2\n", + " location = [-16, -63]\n", + " , zoom_start = 5\n", " , control_scale = True \n", " , tiles = 'Stamen Terrain'\n", ")\n", @@ -211,16 +220,252 @@ "# [[lat_min, lon_min], [lat_max, lon_max]]\n", "m.add_child(raster_layers.ImageOverlay(\n", " data\n", - " , opacity = 0.7\n", + " , opacity = 1\n", " , bounds = [ext[2], ext[0]]\n", " , mercator_project = True\n", - " , colormap = my_cmap\n", + " , colormap = colormap.rgba_floats_tuple\n", ")\n", " )\n", "\n", + "m.add_child(colormap_scale)\n", + "\n", + "from rasterstats import zonal_stats\n", + "\n", + "calcs = zonal_stats(\n", + " shapefile, \n", + " image, \n", + " stats=\"count min mean max median\",\n", + " nodata=-999\n", + ")\n", + "\n", + "media = str(\"%.2f\" % calcs[0]['mean'])\n", + "\n", + "text = 'Average Rainfall: ' + media + ' mm'\n", + "\n", + "folium.Marker(\n", + " [-15, -53]\n", + " , popup = text\n", + " , tooltip = text\n", + ").add_to(m)\n", + "\n", + "from folium import GeoJson\n", + "\n", + "folium.GeoJson(shapefile, name='geojson').add_to(m)\n", + "\n", "m" ] }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# m.save('results/map.html')" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(700, 1500)" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# data.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Driver: GTiff/GeoTIFF\r\n", + "Files: /notebooks/resources/gpm/gpm_1d.20190530.tif\r\n", + "Size is 1500, 700\r\n", + "Coordinate System is:\r\n", + "GEOGCS[\"WGS 84\",\r\n", + " DATUM[\"WGS_1984\",\r\n", + " SPHEROID[\"WGS 84\",6378137,298.2572326660159,\r\n", + " AUTHORITY[\"EPSG\",\"7030\"]],\r\n", + " AUTHORITY[\"EPSG\",\"6326\"]],\r\n", + " PRIMEM[\"Greenwich\",0],\r\n", + " UNIT[\"degree\",0.0174532925199433],\r\n", + " AUTHORITY[\"EPSG\",\"4326\"]]\r\n", + "Origin = (-180.000000000000000,10.000000000000000)\r\n", + "Pixel Size = (0.100000000000000,-0.100000000000000)\r\n", + "Metadata:\r\n", + " AREA_OR_POINT=Area\r\n", + " TIFFTAG_DATETIME=2019:05:31 13:35:19\r\n", + " TIFFTAG_DOCUMENTNAME=/NRTPUB/imerg/gis/2019/05/3B-HHR-L.MS.MRG.3IMERG.20190530-S233000-E235959.1410.V06B.1day.tif\r\n", + " TIFFTAG_IMAGEDESCRIPTION=DOI=10.5067/GPM/IMERG/3B-HH-L/06 DOIauthority=http://dx.doi.org/ DOIshortName=3IMERGHH_LATE Unit=0.1(mm) ScaleFactor=10\r\n", + " TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)\r\n", + " TIFFTAG_SOFTWARE=IDL 8.7.2, Harris Geospatial Solutions, Inc.\r\n", + " TIFFTAG_XRESOLUTION=100\r\n", + " TIFFTAG_YRESOLUTION=100\r\n", + "Image Structure Metadata:\r\n", + " COMPRESSION=LZW\r\n", + " INTERLEAVE=BAND\r\n", + "Corner Coordinates:\r\n", + "Upper Left (-180.0000000, 10.0000000) (180d 0' 0.00\"W, 10d 0' 0.00\"N)\r\n", + "Lower Left (-180.0000000, -60.0000000) (180d 0' 0.00\"W, 60d 0' 0.00\"S)\r\n", + "Upper Right ( -30.0000000, 10.0000000) ( 30d 0' 0.00\"W, 10d 0' 0.00\"N)\r\n", + "Lower Right ( -30.0000000, -60.0000000) ( 30d 0' 0.00\"W, 60d 0' 0.00\"S)\r\n", + "Center (-105.0000000, -25.0000000) (105d 0' 0.00\"W, 25d 0' 0.00\"S)\r\n", + "Band 1 Block=1500x2 Type=UInt16, ColorInterp=Gray\r\n" + ] + } + ], + "source": [ + "# !gdalinfo '/notebooks/resources/gpm/gpm_1d.20190530.tif'" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "79" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# data[0][0]" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "29.003405714285716" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# np.mean(data)" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'0.80'" + ] + }, + "execution_count": 56, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# from rasterstats import zonal_stats\n", + "# calcs = zonal_stats(\n", + "# shapefile, \n", + "# image, \n", + "# stats=\"count min mean max median\",\n", + "# nodata=-999\n", + "# )\n", + "\n", + "# type(calcs)\n", + "# calcs[0]['mean']\n", + "\n", + "# print(round(calcs[0]['mean']))\n", + "# print(\"%.2f\" % calcs[0]['mean'])\n", + "\n", + "# str(\"%.2f\" % calcs[0]['mean'])" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# from folium import IFrame\n", + "\n", + "# text = 'your text here'\n", + "\n", + "# iframe = folium.IFrame(text, width=700, height=450)\n", + "# popup = folium.Popup(iframe, max_width=3000)\n", + "\n", + "# Text = folium.Marker(location=[40.738, -73.98], popup=popup,\n", + "# icon=folium.Icon(icon_color='green'))\n", + "# m.add_child(Text)" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# from folium import GeoJson\n", + "\n", + "# folium.GeoJson(shapefile, name='geojson').add_to(m)\n", + "\n", + "# m" + ] + }, { "cell_type": "code", "execution_count": null, diff --git a/Folium-Final.ipynb b/Folium-Final.ipynb new file mode 100644 index 0000000..b98b21f --- /dev/null +++ b/Folium-Final.ipynb @@ -0,0 +1,306 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import folium\n", + "import matplotlib\n", + "import numpy as np\n", + "import branca\n", + "from matplotlib.pyplot import imread\n", + "from matplotlib.colors import Normalize\n", + "from matplotlib.colors import ListedColormap\n", + "from folium import raster_layers" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "image = '/notebooks/resources/gpm/gpm_1d.20190424.tif'\n", + "shapefile = '/notebooks/resources/centro-oeste.geojson'\n", + "date = str(image[38:40]) + '-' + str(image[36:38]) + '-' + str(image[32:36])" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "10.0 -180.0\n", + "-60.0 -180.0\n", + "-60.0 -30.0\n", + "10.0 -30.0\n" + ] + } + ], + "source": [ + "from osgeo import gdal,ogr,osr\n", + "\n", + "def GetExtent(gt,cols,rows):\n", + " ''' Return list of corner coordinates from a geotransform\n", + "\n", + " @type gt: C{tuple/list}\n", + " @param gt: geotransform\n", + " @type cols: C{int}\n", + " @param cols: number of columns in the dataset\n", + " @type rows: C{int}\n", + " @param rows: number of rows in the dataset\n", + " @rtype: C{[float,...,float]}\n", + " @return: coordinates of each corner\n", + " '''\n", + " ext=[]\n", + " xarr=[0,cols]\n", + " yarr=[0,rows]\n", + "\n", + " for px in xarr:\n", + " for py in yarr:\n", + " x=gt[0]+(px*gt[1])+(py*gt[2])\n", + " y=gt[3]+(px*gt[4])+(py*gt[5])\n", + " ext.append([y,x])\n", + " print (y,x)\n", + " yarr.reverse()\n", + " return ext\n", + "\n", + "def ReprojectCoords(coords,src_srs,tgt_srs):\n", + " ''' Reproject a list of x,y coordinates.\n", + "\n", + " @type geom: C{tuple/list}\n", + " @param geom: List of [[x,y],...[x,y]] coordinates\n", + " @type src_srs: C{osr.SpatialReference}\n", + " @param src_srs: OSR SpatialReference object\n", + " @type tgt_srs: C{osr.SpatialReference}\n", + " @param tgt_srs: OSR SpatialReference object\n", + " @rtype: C{tuple/list}\n", + " @return: List of transformed [[x,y],...[x,y]] coordinates\n", + " '''\n", + " trans_coords=[]\n", + " transform = osr.CoordinateTransformation( src_srs, tgt_srs)\n", + " for x,y in coords:\n", + " x,y,z = transform.TransformPoint(x,y)\n", + " trans_coords.append([x,y])\n", + " return trans_coords\n", + "\n", + "raster=image\n", + "ds=gdal.Open(raster)\n", + "\n", + "gt=ds.GetGeoTransform()\n", + "cols = ds.RasterXSize\n", + "rows = ds.RasterYSize\n", + "\n", + "ext=GetExtent(gt,cols,rows)\n", + "\n", + "src_srs=osr.SpatialReference()\n", + "src_srs.ImportFromWkt(ds.GetProjection())\n", + "# tgt_srs=osr.SpatialReference()\n", + "# tgt_srs.ImportFromEPSG(3857)\n", + "tgt_srs = src_srs.CloneGeogCS()\n", + "\n", + "geo_ext=ReprojectCoords(ext,src_srs,tgt_srs)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "010000" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# create custom colormap\n", + "\n", + "colormap = branca.colormap.StepColormap(\n", + " ['#2b83ba', '#64abb0', '#9dd3a7', '#c7e9ad', '#edf8b9', '#ffedaa', '#fec980', '#f99e59', '#e85b3a', '#d7191c'],\n", + " vmin=0, vmax=10000,\n", + " index=[0, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000],\n", + " caption='rainfall (mm)'\n", + ")\n", + "\n", + "colormap" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# recreate custom colormap setting transparency on the first color\n", + "\n", + "change = list()\n", + "alphacolor = list()\n", + "for item in colormap.colors:\n", + " change.append(list(item))\n", + "\n", + "change[0][3] = 0\n", + "\n", + "for item in change:\n", + " alphacolor.append(tuple(item))" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# set colormap\n", + "colormap = branca.colormap.StepColormap(\n", + " alphacolor\n", + " , vmin=0, vmax=10000\n", + " , index=[0, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000]\n", + " , caption='rainfall(mm)'\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# set scale\n", + "colormap_scale = colormap.to_linear().scale(0, 1000).to_step(10)\n", + "colormap_scale.caption = 'rainfall(mm) ' + date" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# draw map\n", + "m = folium.Map(\n", + " location = [-16, -63]\n", + " , zoom_start = 5\n", + " , control_scale = True \n", + " , tiles = 'Stamen Terrain'\n", + ")\n", + "\n", + "data = matplotlib.pyplot.imread(image)\n", + "\n", + "# Image bounds on the map in the form\n", + "# [[lat_min, lon_min], [lat_max, lon_max]]\n", + "m.add_child(raster_layers.ImageOverlay(\n", + " data\n", + " , opacity = 1\n", + " , bounds = [ext[2], ext[0]]\n", + " , mercator_project = True\n", + " , colormap = colormap.rgba_floats_tuple\n", + ")\n", + " )\n", + "\n", + "m.add_child(colormap_scale)\n", + "\n", + "from rasterstats import zonal_stats\n", + "\n", + "calcs = zonal_stats(\n", + " shapefile, \n", + " image, \n", + " stats=\"count min mean max median\",\n", + " nodata=-999\n", + ")\n", + "\n", + "average = str(\"%.2f\" % calcs[0]['mean'])\n", + "\n", + "text = 'Average Rainfall: ' + average + ' mm (' + date +')'\n", + "\n", + "folium.Marker(\n", + " [-15, -53]\n", + " , popup = text\n", + " , tooltip = text\n", + ").add_to(m)\n", + "\n", + "from folium import GeoJson\n", + "\n", + "folium.GeoJson(shapefile, name='geojson').add_to(m)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Folium_2.ipynb b/Folium_2.ipynb index 6f73ef3..a8b584c 100644 --- a/Folium_2.ipynb +++ b/Folium_2.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -11,7 +11,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -29,7 +29,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -38,9 +38,22 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "10.0 -180.0\n", + "-60.0 -180.0\n", + "-60.0 -30.0\n", + "10.0 -30.0\n", + "ext = [[10.0, -180.0], [-60.0, -180.0], [-60.0, -30.0], [10.0, -30.0]]\n", + "geo_ext = [[10.0, -180.0], [-60.0, -180.0], [-60.0, -30.0], [10.0, -30.0]]\n" + ] + } + ], "source": [ "from osgeo import gdal,ogr,osr\n", "\n", @@ -110,16 +123,64 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Driver: GTiff/GeoTIFF\r\n", + "Files: /notebooks/resources/gpm/gpm_1d.20190531.tif\r\n", + " /notebooks/resources/gpm/gpm_1d.20190531.tif.aux.xml\r\n", + "Size is 1500, 700\r\n", + "Coordinate System is:\r\n", + "GEOGCS[\"WGS 84\",\r\n", + " DATUM[\"WGS_1984\",\r\n", + " SPHEROID[\"WGS 84\",6378137,298.2572326660159,\r\n", + " AUTHORITY[\"EPSG\",\"7030\"]],\r\n", + " AUTHORITY[\"EPSG\",\"6326\"]],\r\n", + " PRIMEM[\"Greenwich\",0],\r\n", + " UNIT[\"degree\",0.0174532925199433],\r\n", + " AUTHORITY[\"EPSG\",\"4326\"]]\r\n", + "Origin = (-180.000000000000000,10.000000000000000)\r\n", + "Pixel Size = (0.100000000000000,-0.100000000000000)\r\n", + "Metadata:\r\n", + " AREA_OR_POINT=Area\r\n", + " TIFFTAG_DATETIME=2019:06:01 13:35:13\r\n", + " TIFFTAG_DOCUMENTNAME=/NRTPUB/imerg/gis/2019/05/3B-HHR-L.MS.MRG.3IMERG.20190531-S233000-E235959.1410.V06B.1day.tif\r\n", + " TIFFTAG_IMAGEDESCRIPTION=DOI=10.5067/GPM/IMERG/3B-HH-L/06 DOIauthority=http://dx.doi.org/ DOIshortName=3IMERGHH_LATE Unit=0.1(mm) ScaleFactor=10\r\n", + " TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)\r\n", + " TIFFTAG_SOFTWARE=IDL 8.7.2, Harris Geospatial Solutions, Inc.\r\n", + " TIFFTAG_XRESOLUTION=100\r\n", + " TIFFTAG_YRESOLUTION=100\r\n", + "Image Structure Metadata:\r\n", + " COMPRESSION=LZW\r\n", + " INTERLEAVE=BAND\r\n", + "Corner Coordinates:\r\n", + "Upper Left (-180.0000000, 10.0000000) (180d 0' 0.00\"W, 10d 0' 0.00\"N)\r\n", + "Lower Left (-180.0000000, -60.0000000) (180d 0' 0.00\"W, 60d 0' 0.00\"S)\r\n", + "Upper Right ( -30.0000000, 10.0000000) ( 30d 0' 0.00\"W, 10d 0' 0.00\"N)\r\n", + "Lower Right ( -30.0000000, -60.0000000) ( 30d 0' 0.00\"W, 60d 0' 0.00\"S)\r\n", + "Center (-105.0000000, -25.0000000) (105d 0' 0.00\"W, 25d 0' 0.00\"S)\r\n", + "Band 1 Block=1500x2 Type=UInt16, ColorInterp=Gray\r\n", + " Min=0.000 Max=2237.000 \r\n", + " Minimum=0.000, Maximum=2237.000, Mean=34.214, StdDev=111.334\r\n", + " Metadata:\r\n", + " STATISTICS_MAXIMUM=2237\r\n", + " STATISTICS_MEAN=34.214233333333\r\n", + " STATISTICS_MINIMUM=0\r\n", + " STATISTICS_STDDEV=111.33382267193\r\n" + ] + } + ], "source": [ "!gdalinfo '/notebooks/resources/gpm/gpm_1d.20190531.tif'" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -128,7 +189,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -149,11 +210,25 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "metadata": { "scrolled": false }, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "m = folium.Map(\n", " location = [-22, -114]\n", @@ -195,18 +270,76 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "{'_children': OrderedDict([('stamenterrain',\n", + " ),\n", + " ('image_overlay_92155413a6b84d5387a381702772f37a',\n", + " ),\n", + " ('marker_ac7d26db7daf4a5789ade1e0b5287481',\n", + " ),\n", + " ('marker_29b98858dab84a42ba7fd08322bdfefe',\n", + " )]),\n", + " '_env': ,\n", + " '_id': '5302fb64bbe34b86b83a04708129f218',\n", + " '_name': 'Map',\n", + " '_parent': ,\n", + " '_png_image': None,\n", + " 'control_scale': True,\n", + " 'crs': 'EPSG3857',\n", + " 'global_switches': ,\n", + " 'height': (100.0, '%'),\n", + " 'left': (0.0, '%'),\n", + " 'location': [-22, -114],\n", + " 'max_bounds': False,\n", + " 'max_lat': 90,\n", + " 'max_lon': 180,\n", + " 'min_lat': -90,\n", + " 'min_lon': -180,\n", + " 'no_wrap': False,\n", + " 'objects_to_stay_in_front': [],\n", + " 'png_enabled': False,\n", + " 'position': 'relative',\n", + " 'top': (0.0, '%'),\n", + " 'width': (100.0, '%'),\n", + " 'world_copy_jump': False,\n", + " 'zoom_control': True,\n", + " 'zoom_start': 2}" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "vars(m)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 2 3 3 ... 0 0 0]\n", + " [ 4 6 3 ... 0 0 0]\n", + " [10 9 3 ... 0 0 0]\n", + " ...\n", + " [ 0 0 0 ... 0 0 0]\n", + " [ 0 0 0 ... 0 0 0]\n", + " [ 0 0 0 ... 1 0 0]]\n", + "(700, 1500)\n" + ] + } + ], "source": [ "print(data)\n", "print(data.shape)\n", @@ -217,9 +350,54 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "\n", + "## ds ##:\n", + "\n", + " >\n", + "\n", + "\n", + "## data ##:\n", + "\n", + "[[ 2 3 3 ... 0 0 0]\n", + " [ 4 6 3 ... 0 0 0]\n", + " [10 9 3 ... 0 0 0]\n", + " ...\n", + " [ 0 0 0 ... 0 0 0]\n", + " [ 0 0 0 ... 0 0 0]\n", + " [ 0 0 0 ... 1 0 0]]\n", + "\n", + "\n", + "## gt ##:\n", + "\n", + "(-180.0, 0.1, 0.0, 10.0, 0.0, -0.1)\n", + "\n", + "\n", + "## proj ##:\n", + "\n", + "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.2572326660159,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]\n", + "\n", + "\n", + "## inproj ##:\n", + "\n", + "GEOGCS[\"WGS 84\",\n", + " DATUM[\"WGS_1984\",\n", + " SPHEROID[\"WGS 84\",6378137,298.2572326660159,\n", + " AUTHORITY[\"EPSG\",\"7030\"]],\n", + " AUTHORITY[\"EPSG\",\"6326\"]],\n", + " PRIMEM[\"Greenwich\",0],\n", + " UNIT[\"degree\",0.0174532925199433],\n", + " AUTHORITY[\"EPSG\",\"4326\"]]\n" + ] + } + ], "source": [ "# First: read the geotiff image with GDAL.\n", "from osgeo import gdal, osr\n", @@ -243,9 +421,17 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(0, 0, 0, 0)\n" + ] + } + ], "source": [ "q = lambda x: (0, 0, 0, 0)\n", "print(q(0))" @@ -253,9 +439,18 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(0.72312, 0.11873, 1.0, 1.0)\n", + "(0.72312, 0.11873, 1.0)\n" + ] + } + ], "source": [ "cmap = colorcet.cm.bmw\n", "\n", @@ -267,9 +462,32 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+0XWWd3/H359x7k/BTIkGgJAprTB1TVMA06KKjFGVNQAudgtPAaGWWTpxZRPFXXTCdotLVNXXa6rQl/rgVKjoqOvhjRRubsgZYjrMUEyDgJIjNoqME6ED4JSgkufd++8feNx4P5569z7l7n2ffcz+vtfZa55y97/M8B777m32e/eznUURgZmaLSyt1A8zMbPic/M3MFiEnfzOzRcjJ38xsEXLyNzNbhJz8zcwWodLJX9KYpLskfbvLvqWSviJpj6TbJZ1cZSPN6uTYtsWonyv/K4B759j3DuCJiHgp8AngY/NtmNkQObZt0SmV/CWtBN4EfHaOQy4Ebshf3wS8QZLm3zyzejm2bbEaL3ncnwMfAo6aY/9JwAMAETEl6SngWGBf+0GSNgIbAWiNv1rLjhmgyfPzquPKfuXq3f3IwST1Hvvs00nqBdjHgX0Rcdxc+1tHrwymnissJ559bFtErK+0cZmaYnt5DU3t7egl4jdWzfmfujZ37fl/Q6+z3RlrTk5S7x133DFnbJeNa6g1tnsqzISS3gw8EhF3SDp7PpVFxCQwCdA64rhYsuai+RQ3kNv+cMXQ65x13LUPJan3op23JKkX4DP89Kc9D5h6jvGXXVBYzsGd/6Py/3H1xfaLYtmpw4/tf7Jqgm98YtPQ6z3igv8w9Drb7dhxfZJ6Jc0d2yXjGuqJ7TLKXAafBVwg6XxgGXC0pL+IiLe2HfMgsArYK2kceAHwWOWttdEjodZYqtod21aPtHFdSmGff0RcFRErI+JkYANwS8fJAbAFeHv++uL8GM8YZyWI1viSwq0Ojm2rT7m4riu2yxi4A1zSNcCOiNgCXAd8QdIe4HGyE8msWAOvkBzbNm8NjOtOfSX/iLgNuC1/fXXb588Bb6myYbY4CNBY+pPEsW1Vakpc95Ju6IsZgESr4VdIZn1bAHHt5G/JNf3nsdkgmh7XTv6W1gLoGzXr2wKIa0/sZkkJ0RqfKNxKlSWtl3RfPg/PlV32v0TSX0m6R9Jt+dO9ZpUrG9cpY9vJ39LKr5CKtuJiNAZsBs4D1gCXSFrTcdh/Aj4fEa8ErgH+tOJvY5YpGdcpY9vJ35Kr4gQB1gF7IuL+iDgA3Eg2L0+7NcDs4863dtlvVpmqkj81xbaTv6UlobGxwg1YIWlH27axo6RDc/Dk9uaftbsb+Bf5698BjpJ0bB1fyxa5knGdMrZ9w9eSEqVHReyLiLXzrO6DwLWSLgO+SzZ1w/Q8yzR7nj7iGhLFtpO/paUWY9U84j47B8+slflnh0TEQ+RXR5KOBC6KiCerqNzs16jFeHVTN9QS207+lpYqGw+9HVgt6RSyE2MDcOmvVSWtAB6PiBngKiDNdJA28lRdXENNse0+f0tKVDPaJyKmgE3ANrJVub4aEbskXSNpdm7ds4H7JP0EOB749/V8K7O++vx7qiu2feVvyVV1hRQRW4GtHZ+1z9NzE9lqXGb1qvbKv5bYdvK3tBbAk5Bm/Wt+XDv5W2LNP0nM+idQs+Payd+SkkRrIt2CFmZ1kMT4kqWpm9GTk7+l5W4fG0ULIK6d/C25pp8kZoNQq9mDKZ38LblWS6mbYFYpqflxXfhPk6Rlkn4o6W5JuyR9tMsxl0l6VNLOfHtnPc21USMJtYq3mup2bFtNRKjclkqZK//9wDkR8YykCeB7kr4TET/oOO4rEbGp+ibaqBsbS/bz2LFttZBgyUSzuzMLk39EBPBM/nYi36LORtkiImq7si/i2LY6pYrrskpdckkak7QTeAS4OSJu73LYRfkqMjdJWtVlv9nzZLMfpun2Ace21aRkd2bKfyBKJf+ImI6I08hmk1sn6dSOQ74FnJyvInMzcEO3ciRtnJ2zOqaem0+7bWSIloq3utQR20w9W1t7bWEQlIrrOmO7SF+drfkUobcC6zs+fywi9udvPwu8eo6/n4yItRGxVuPLBmmvjRqlvfKfVWVsM35YvY21BWHBX/lLOk7SMfnrw4BzgR93HHNi29sLyGaeMysl4Wgfx7bVIoCZklsqZUb7nAjckC8i3CKbTvTbkq4BdkTEFuA9+dSiU8DjwGV1NdhGiwRj48mufhzbVouWYOmSBf6QV0TcA5ze5fP26USvIltAwKxvStTv6di22ihtl04ZfsLXkpLU+CchzQbR9Lhu9u8SWxSq6vOXtF7SfZL2SLqyy/4XS7pV0l350M3zK/8yZvxqeocyW7nyqo9tJ39Lrorkn/fbbwbOA9YAl0ha03HYn5D1659Otg7qJyv+KmaHVDXap67YdrePpSWqGuu8DtgTEfcDSLoRuBDY3XZMAEfnr18APFRFxWadApiubihPLbHt5G9JCdEaL/UDdIWkHW3vJyNisu39ScADbe/3Amd2lPER4H9LejdwBPDG/ltsVqwlsWxp6fSaJLad/C2t8lPf7ouItfOs7RLgcxHxnyW9FviCpFMjIuVwaxtRfdzwTRLbTv6WXEVDPR8E2ufdWZl/1u4d5E/wRsT3JS0DVpDN62NWGQlU3R3VWmLbN3wtqWxit+KthO3AakmnSFpCdtNrS8cxPwPeACDp5cAy4NHKvoxZmwqnd6gltn3lb2lVtOJRRExJ2gRsA8aA6yNiV8fTuh8A/ruk95HdILssn9bZrFIBTFcUWnXFtpO/JSZaFS3mEhFbga0dn7U/rbsbOKuSysx6aEksW1Jdeq0jtp38LamFsNapWd8WQFw7+VtyTZ8DxaxfwsnfrCcJxhp+kpj1T42Payd/S67pJ4lZvxbCRY2TvyWlBXCFZNavKkf71MXJ35KSYEm56R3MFowxwVHLmp1em906G3kSjPvK30aMpMbHtZO/JSWa3zdqNoimx7WTv6Ul9/nb6Mlu+Da7O7OwdZKWSfqhpLsl7ZL00S7HLJX0lXyVmdslnVxHY230ZFf+rcKtlrod21aT2V+0ZbZUylz57wfOiYhnJE0A35P0nYj4Qdsx7wCeiIiXStoAfAz4lzW010ZQwhPAsW21CGBqZoGP9sknB3omfzuRb53f6kKyxQQAbgKulSRPmmVFWlKy0T51xnaKyI+A+OWTw684sbds/pvUTXieljQao33yNSTvAF4KbI6I2zsOObTSTD4D3VPAscC+jnI2AhsBli4/ntf+3qXza/0Ajvqnhw+9zll7n/hCknrf9M8niw+qy0d+u/CQsWrm8x9IHbF95LEncMkVf1B305/nt44PJk5cVXxgxZ7+0rvQ4ccMvV6AS7/ROa19MwgYr3BC/zqUal1ETEfEaWSLCKyTdOoglUXEZESsjYi1E0ekCRZrltknIVP1i9YR28uOXl5tI23Bkcr196fs8+/rn6aIeBK4lXzFmDaHVpqRNE62gPBjVTTQRl8TThDHtlVtwSd/ScdJOiZ/fRhwLvDjjsO2AG/PX18M3OL+fitj9iGvoq1cWVov6b58ZM6VXfZ/QtLOfPuJpKcc21aX6ZkotZUxQGwX3vwp0+d/InBD3jfaAr4aEd/uWEXmOrIFg/cAj5MtM2ZWSFRzwzePz81kCXwvsF3SlnyRCwAi4n1tx78bOAe41bFtVWtVOL3DgLF9elG5ZUb73NOtoI5VZJ4D3lJUllmnCmc/XAfsiYj7s3J1I9lInd1zHH8J8OGIuLlzh2Pb5kvVPrw4UGwXFdrssUg28vqY3mGFpB1t7ycjon0Y06FRObm9wJld65ReApwC3NJfa83K6XPakiSx7eRvaZW/8t8XEWsrqnUDcFNETFdUntmv6+8XbZLYdvK3pCqcz//QqJzcyvyzbjYAl1dRqVk3FU9YWEtsO/lbchWdJNuB1ZJOITsxNgDPe4pQ0m8Cy4HvV1GpWTcBTFU3KKyW2Hbyt6RaFS3mkj99uwnYBowB10fEro6RO5CdODd6uKbVaUziqKXVpNe6YtvJ39KqcK3TiNgKbO347OqO9x+ppDKzXipew7eO2Hbyt6SEks7tY1YHkXbOqjKc/C25VsNPErNBNH2NIid/Syq7QkrdCrNq+crfrIig1fRLJLM+jcRiLmZ1EjDR8LVOzfrVkjhiyVjqZvTk5G9JudvHRpG7fcyKSO72sdGzALoznfwtKeHRPjZ6srhO3YrenPwtOXf72KiJoPRCLak4+VtSEkyM+YavjZZWSxy+pNnptdmts5Hnbh8bRQthIIOTvyXX9JPEbBBN7/Mvs4D7Kkm3StotaZekK7occ3a+GPbsAsJXdyvLrJMQLRVvtdTt2LaazP6iTRXbZZS58p8CPhARd0o6CrhD0s3tiwfn/joi3lx9E22kVTz7YZ8c21aPtHFdSpkF3B8GHs5fPy3pXrI1JedaPNistJRD4hzbVpuA6ekRGu0j6WTgdOD2LrtfK+lu4CHggxGxq8vfbwQ2Aixdfny/bbUR1JTpHaqM7SOPPaG+htqC0JI4vOHTO5Q+6yQdCXwNeG9E/Lxj953ASyLiVcB/A77ZrYyImIyItRGxduKIYwZts40SwVireCtVlLRe0n2S9ki6co5jfretj/9L+WeVxvayo5eXa7CNtJbKbWUMGtu9lLrylzRBdnJ8MSK+3rm//YSJiK2SPilpRUTsK1O+LV5VDfWUNAZsBs4F9gLbJW1p77+XtBq4CjgrIp6Q9CLHttVC1Q1hHjS2i8otM9pHwHXAvRHx8TmOOSE/Dknr8nIfK/5aZtlKXkVbCeuAPRFxf0QcAG4ELuw45g+AzRHxRP7+URzbVoPZid1SxXZEPFJUaJkr/7OAtwE/krQz/+yPgRfnlXwauBj4I0lTwLPABi+QbWX0ceW/QtKOtveTETHZ9v4k4IG293uBMzvK+IcAkv6GbCHsv8SxbTUIYLp8mNQR2x+JiP/Vq9Iyo32+R3aO9jrmWuDaorLMOmXTO5RK/vsiYu08qxsHVgNnAyuB7wLLI+LJuf7AsW2DaEkcNlH6hm8tsS3pFb1iO/0wC1v0pOKthAeBVW3vV+aftdsLbImIgxHxf4GfkJ0wZpUrE9cpY9vJ35JrocKthO3AakmnSFoCbAC2dBzzTbIrIyStIPupfH9138QsI8rFdcrYdvK3pEQ1V0cRMQVsArYB9wJfjYhdkq6RdEF+2DbgMUm7gVuBfx0Rvnlrtajqyr+u2PbEbpZcVU/4RsRWYGvHZ1e3vQ7g/flmVqsqn1yvI7ad/C2t8v2eZguKF3Mx60GUHutstmAIWNLwRYqc/C25hk9+aNa/PqZuSMXJ35Jr+Dli1jfR/Lh28rekvIyjjaqmx7WTvyXX8HPEbCBNj2snf0uu2bfFzAYz49E+ZnPTAljuzqxfEiwZb/ZljZO/Jdf0n8dmg2h6WDv5W1LZHChmo6fpFzVO/pacmn6WmPVpds6qJnPyt7QWwMMwZgPxDV+zuWXL3aVuhVm1hBj39A5mvbnbx0ZR08O6zALuqyTdKmm3pF2SruhyjCT9V0l7JN0j6Yx6mmujJnvCt3irpW7HttVIJbdUyvwumQI+EBFrgNcAl0ta03HMeWRLhq0GNgKfqrSVNtKqOkEkrZd0X56or+yy/zJJj0ramS/Y/rs4tq0OJRdyKfvroN/YlvTOojLLLOD+MPBw/vppSfeSrSa/u+2wC4HP5wsK/EDSMZJOzP/WrAdVMgeKpDFgM3Au2Xqm2yVtiYjdHYd+JSI2df69Y9uqNLuMYyVlzTO259JXn7+kk4HTgds7dp0EPND2fm/+2a+dIJI2kl09sXT58f1UbaOqusVc1gF7IuJ+AEk3kiXuzhOkezMqjO0jjj2B/VPTfX+B+ZqegakE9X7rzoe44+8fKD6wBr84eDRHjC9JUnexykb7zCu251I6+Us6Evga8N6I+PkglUXEJDAJ8PJXnhZ/+junDlLMvPziuMOGXues8546L0m9t131W0nqBTj6I733KwLNlEpYKyTtaHs/mcfTrG5J+swu5Vwk6XXAT4D3RcQDVcf2K047I65+08sHKWZejpgQM4cNPxF+8ZFjUaKBLb931km8Zd2Lk9T9tXf32hu0yif/2mK7V6Wlkr+kCbKT44sR8fUuhzwIrGp7vzL/zKyQYqbMYfsiYu08q/oW8OWI2C/pXcANkn4bx7ZVLYBycQ01xTZwTq8/KDPaR8B1wL0R8fE5DtsC/Kt8ZMRrgKfcJ2rlRHaSFG3FCpN0RDwWEfvzt58FXo1j2+oSUW4rNmhs91Tmyv8s4G3Aj/IREgB/DLw4r/TTZKvKnw/sAX4J/H6Jcs0y5U6AItuB1ZJOITsxNgCXth/QcaP2ArKfz45tq0llff6DxPa9RYWWGe3zPQpG2+UjIS4vKsvseSL6+Xnco5iYkrQJ2AaMAddHxC5J1wA7ImIL8B5JF5ANX34cuCgiflxQrmPbBhComouaQWP7sqJy/YSvJVeyz79QRGwlu1Jv/+zqttdXAVdVUplZLxHE9MEKi6s+tp38LbGAmanUjTCrVnVDmGvj5G9p9TcqwmxhCKq6l1UbJ39LLGDGyd9GkJO/WW9V9fmbNYmqG+1TCyd/S8/J30ZNzBAHD6RuRU9O/pZWBJSb3sFs4ZDQ2FjqVvTk5G/JudvHRo2Cysb518XJ3xKr5iEvs2Zpflw7+Vt6DT9JzAbjK3+zuVU0vYNZ47jbx2xuwn3+NoJiBg4+l7oVPTn5W2IB0x7tYyNGgrGJ1K3oycnf0vL0Djai/JCXWQF3+9jI8dw+ZkV8w9dGlJO/WQEnfxtJzU7+hWv4mtVqdnqHoq0ESesl3Sdpj6Qrexx3kaSQNN9Fs826i2niuV+U2sqoI7bLLOB+vaRHJP3tHPvPlvSUpJ35dnW348y6C2LqYOFWRNIYsBk4D1gDXCJpTZfjjgKuAG4H/q1j22qhFkwsLbcVFTVYbBcqc+X/OWB9wTF/HRGn5ds1ZSo2A7JfxtVc+a8D9kTE/RFxALgRuLDLcf8O+BjwHPBtHNtWh9mHF8tsxQaJ7UKFyT8ivku2ILBZ5YIgpqcLN2CFpB1t28aOok4CHmh7vzf/7BBJZwCrIuJ/5h/dhWPb6jIT5bZ6YrtQVTd8XyvpbuAh4IMRsavbQfmX2ghwwkkrK6raFrSg7Epe+yJi4D56SS3g48Blff5p37H9D1auGrSZNkrKj/ZJEttVJP87gZdExDOSzge+CazudmBETAKTAC9/5WnNvhVuQ1LZfP4PAu1Zd2X+2ayjgFOB25StrH0CsAV4V48yB4rtV5x2hmN7kYuZaeLZn1dV3ECxLemCiNgxV6HzTv4R8fO211slfVLSiojYN9+ybRGIKHVDt4TtwGpJp5CdGBuAS39VTTwFrJh9L+k24IPAnHHq2LaBqQVLDq+qtIFiu1fihwqGeko6Qfk/N5LW5WU+Nt9ybbGI7CqpYCssJWIK2ARsA+4FvhoRuyRdI+mCQVrm2Lb5iIhSW4lyKo9tKHHlL+nLwNlkNyX2Ah8GJvJGfRq4GPgjSVPAs8CGKPONzOBXo32qKCpiK7C147OuwzMj4mzHttWqwocX+43tMmUWJv+IuKRg/7XAtWUqM3u+KHvDt/qaHdtWp4ZfJnh6B0srmB3KaTY6FsAiRU7+llhlo33MmmNmmvjFk6lb0ZOTv6VV3Wgfs+ZojaHDjk7dip6c/C0xX/nbKHK3j1lvFY72MWuUhg8Mc/K3pIIgEo32MauNV/IyK+Arf7MknPwtrQji4IHUrTCrXNN/0Tr5W2LpHvIyq0/Q9Ke8nPwtPXf72Chyn79ZDxGlJm4zW2ganvud/C29pveNmg2m2dnfyd/SiiCmnfxttMTUFNNPNnv2byd/SyoimDk4lboZZpXS2Dito45N3Yye5r2Yi9m8BMT0TOFWhqT1ku6TtEfSlV32/6GkH0naKel7ktZU/n3MmH3Gq5rFXKCe2Hbyt+SqSP6SxoDNwHnAGuCSLifAlyLiFRFxGvBnZItem9UiSm5F6optd/tYUhHBTDXz+a8D9kTE/QCSbgQuBHa31dW+ovYRNP2OnC1s1Q33qSW2nfwtuZKjfVZIal+QejIiJtvenwQ80PZ+L3BmZyGSLgfeDywBzum/tWYl9PeMV5LYLrOG7/XAm4FHIuLULvsF/BfgfOCXwGURcWdRuWZAP6N99kXE2vlXF5uBzZIuBbZKegGObatYTB9k6olHyx5eR2z/CfD2XseXufL/HNk6pp+fY/95wOp8OxP4FF3+VTLrpsLRPg8Cq9rer8w/m8uNwCTwOhzbVjGNTdA65riqihsktj9VVGjhDd+I+C7weI9DLgQ+H5kfAMdIOrGoXLNZM9MzhVsJ24HVkk6RtATYAGxpP0DS6ra3bwJ+jGPb6lLVHd/BYvv/FBVaRZ9/t/6ok4CHOw+UtBHYCHDCSSsrqNoWvHyo57yLiZiStAnYBowB10fELknXADsiYguwSdIbgYPAExT8LGYesf3Uc8N/dmHs4DQHDzwz9HpnpqfIesiGL2ZmmJlu6HMiFd3vrSm2h3vDN7+JMQnwm684LQ5OD3+wxYPvuWTodc66Ztv9Sepd8qHvJqm3lAqf8I2IrcDWjs+ubnt9ReffSDq5oroPxfapp50RRy4Z/liK/Xt/xtgLDx96vTe8/giWrTx56PUCMDHGeCvNPzzDNEhsF6kiQvvtjzI7JGj03D6ObRvIzNRBDjz+SOpm9FRF8p/9yXEj2c2wpyLieT+LzbqKYOZAQ3+2O7ZtQBqfYHz5i1I3o6cyQz2/DJxNNhZ1L/BhYAIgIj5N9lPkfGAP2XC436+rsTaCAmYSXfk7tm0xK0z+EdGzkzyyySkur6xFtqgE6Wb1dGzbYuYnfC2tgKhmegezhmn2jWgnf0ssmnzD12xkOflbWhWN8zdrkjh4kAOP/n3qZvTk5G9JRQTTzR3tYzYQjU8wvuKE1M3oycnfEnO3j1kKTv6Wlrt9zJJw8re0AiLBNB9m9fNoH7M5BVF21k4zq5CTv6UVEDO+8rfRElNTHHi09GIuSTj5W1IRMH3AD3nZaNH4BBPHebSP2dwi3Odvo0dAojUOyipcycusbjPTUbiVIWm9pPsk7ZF0ZZf975e0W9I9kv5K0ksq/zJmh6jkVqKkGmLbyd/Syod6Fm1FJI0Bm8nW3V0DXCJpTcdhdwFrI+KVwE3An1X8bcwqV1dsu9vHkgpgppobvuuAPRFxP0A+B/+FwO5DdUXc2nb8D4C3VlGxWaeZg1Psf3RfVcXVEttO/pZWRFU3fLutt3tmj+PfAXyniorNOrUmJlj6ouOrKq6W2Hbyt6Si/ENeKyTtaHs/ma+b2zdJbwXWAq8f5O/NSil/vzdJbDv5W1rlk/++iFjbY3+p9XYlvRH4N8DrI2J/P00164tK31JNEttO/pZYZU/4bgdWSzqF7MTYAFzafoCk04HPAOsjotmra9sIqGyoZy2xXeqfphLDjC6T9Kiknfn2zjLlms0+4Vu0FRYTMQVsArYB9wJfjYhdkq6RdEF+2H8EjgT+Mo/TLY5tq4eycf5ltgKDxnZRuWUWcJ8dZnQu2Y2G7ZK2RMTujkO/EhGbCr+JWZuA0uP4C8uK2Eq26Hr7Z1e3vX5j+748tn+CY9sqNjM1xf59lY326Tu2yyjT7VM4zMhsYBHMpJvewbFttWhNjFc52qcWZbp9ug0zOqnLcRflT5fdJGlVl/1mzxNR3RO+A3BsW42qe8K3DlU94fst4OT86bKbgRu6HSRpo6QdknY8+fhjFVVtC13MzBRuCfUd2088Vt3PfVuoBK2SWyJlkn/hMKOIeKxtaNFngVd3KygiJiNibUSsPeaFxw7SXhs1UXzVX+OVfy2xvfzYFbU01haahX/lf2iYkaQlZMOMfu1OsqQT295eQHZH2qxYPs6/aKuJY9vqo1a5LZHCG74RMSVpdpjRGHD97DAjYEdEbAHekw85mgIeBy6rsc02QoJ0a/g6tq0uM9PTHHj8idTN6KnUQ14lhhldBVxVbdNsUYhg+kC6Pn3HttWhNT7OsuOa3f3nJ3wtqQiYCS/mYiMo4c3cMpz8LblpJ38bNQtgJS8nf0sqAK/iaKNHSW/mluHkb8n5yt9GTUxPs//pp1M3oycnf0tqJuBANSt5mTWGxsdY+sLlqZvRk5O/JeduHxs9Qi13+5jNKQh3+9joER7tY9aLb/jayPINX7PenPxt9DR/tE+zW2cjLyIb7VO0lVFiVa7XSbpT0pSkiyv/Mma5mJlh6tlfltrKqCO2feVvSQXVjPYpueLcz8jm5vngvCs066E1NsaSFxxdSVl1xbaTvyVVYZ9/4apcEfF3+b6kCwTYIiCgutE+tcS2k78lV7JbZ4WkHW3vJyNisu19t1W5zqygeWYDaZWf3iFJbDv5W1JZn3+pQ/dFxNqam2NWEfVz5Z8ktp38LbmKxvkXrsplNjTVdvvUEttO/pZUABV1wB9alYvsxNgAXFpN0WZ9imBmaqqq0mqJbSd/SyqISkb7lFmVS9I/Br4BLAf+maSPRsQ/mnflZh3UajFx+LJKyqortp38LalstE81w31KrMq1newns1ntqpzbp47YdvK3tMrf8DVbONT8id1Kta7E02VLJX0l33+7pJOrbqiNptkr/yqe8B2EY9vqolar1JZKYc1tT5edB6wBLpG0puOwdwBPRMRLgU8AH6u6oTa6pqN4q4Nj2+qi/Mp/QSd/2p4ui4gDwOzTZe0uBG7IX98EvEFq+AKW1ggzZNM7FG01cWxbbcbGxkptqZTp8y/zdNmhY/I7008BxwL72g+StBHYmL/d/7rfWPG3gzS6AivoaNtI13v08nR1w8t67dzHgW2f4acrSpRTR9tri+2XHX9UithO9f84Zd0pv/OcsX3HHXdsk1QmriFR+4d6wzd/ZHkSQNKOVE9spqp7sX7nXvsjYv2w2lKnJsT2Yo2vlN95rn0LIa7LdPuUebrs0DGSxoEXAI9V0UCzGjm2bdEqk/wPPV0maQnZ02VbOo7ZArw9f30xcEuE1+azxnNs26JV2O1T5uky4DrgC5L2AI9bfypoAAACNElEQVSTnURFJosPqU2quv2dG2QEY9vxtXjqnjf5IsbMbPFp9iNoZmZWCyd/M7NFKEnyL3qkvsZ6r5f0iKShjsGWtErSrZJ2S9ol6Yoh1btM0g8l3Z3X+9Fh1NtW/5ikuyR9e5j1puK4Hk5c53U7tudp6Mm/5CP1dfkckGL87RTwgYhYA7wGuHxI33k/cE5EvAo4DVgv6TVDqHfWFcC9Q6wvGcf1UOMaHNvzluLKv8wj9bWIiO+SjdgYqoh4OCLuzF8/TRY0Jw2h3oiIZ/K3E/k2lDv8klYCbwI+O4z6GsBxPaS4zutzbM9TiuTf7ZH6oQRME+SzQp4O3D6k+sYk7QQeAW6OiKHUC/w58CEqW6ir8RzXQ4zrvE7H9jz4hu8QSToS+Brw3oj4+TDqjIjpiDiN7OnVdZJOrbtOSW8GHomIO+quy9JLEdfg2J6vFMl/US60LWmC7AT5YkR8fdj1R8STwK0Mp2/4LOACSX9H1v1xjqS/GEK9KTmuE8Q1OLYHlSL5l3mkfqTkUwBfB9wbER8fYr3HSTomf30YcC7w47rrjYirImJlRJxM9v/3loh4a931Jua4Hm7dju15Gnryj4gpYPaR+nuBr0bErmHULenLwPeBl0naK+kdw6iX7GrhbWRXCTvz7fwh1HsicKuke8iS080RsWCHpjWZ43qocQ2O7Xnz9A5mZouQb/iamS1CTv5mZouQk7+Z2SLk5G9mtgg5+ZuZLUJO/mZmi5CTv5nZIvT/AWShD23TQIGcAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "import numpy as np\n", "import matplotlib.pylab as pl\n", @@ -302,9 +520,45 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[0.0000e+00 0.0000e+00 0.0000e+00 1.0000e+00]\n", + " [2.7065e-02 2.1430e-05 0.0000e+00 1.0000e+00]\n", + " [5.2054e-02 7.4728e-05 0.0000e+00 1.0000e+00]\n", + " ...\n", + " [1.0000e+00 9.9953e-01 8.7115e-01 1.0000e+00]\n", + " [1.0000e+00 9.9989e-01 9.3683e-01 1.0000e+00]\n", + " [1.0000e+00 1.0000e+00 1.0000e+00 1.0000e+00]]\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+0nVWd3/H35978UkBAkgKTBIIlajOoMFyDLrpmWChtQIXV4kwDSyuztLEuGfFXZ8F0ikqnq+NMq+OUjM4tUvHHiA46NmqclNawkFYwN/xyEkQzoBJAIQQDiNzk3vvtH+e5mcPJuefZ99znOfvk3M9rrWet82PfvfeB7/nmOfvZz96KCMzMbH4Zyt0BMzPrPSd/M7N5yMnfzGwecvI3M5uHnPzNzOYhJ38zs3koOflLGpZ0l6RvtHlvsaQvSdol6Q5Jq6rspFmdHNs2H83mzP8K4L4Z3ns78GREnAp8HPjoXDtm1kOObZt3kpK/pBXAG4DrZihyEXBD8fgm4HWSNPfumdXLsW3z1YLEcn8G/D5w1AzvLwceAoiICUn7gOOAPc2FJG0ANgAcccQLz3z5y1/STZ/n5Ofb/67nbU47/swzs7S7ffv2LO0W9kTEspneXLduXezZs2emtw/avn37lohYV2nPGgYmtsf3PQ0s7Hm7i48+uudt9oPt27fPGNupcV3UU1dsd1Sa/CW9EXgsIrZLOmcujUXEKDAKMDLyihgb2zSX6rrysQW9/1JOe//YWJZ2M5+o/qTTm3v27GEs4b+LpKWV9egf6hyo2H7wW7cAJ/S83VPOP7/nbfYDSTPGdmpcF/VUHtspUs78zwYulHQBsAR4kaTPR8Rbmso8DKwEdktaABwNPFF5b20ATQK/yNW4Y9tqMgk8mbsTHZUm/4i4CrgKoDg7+mDLlwNgE/A24LvAm4Fvh1eMsyTjwANZWnZsW33GIfLEdarUMf9DSLoGGIuITcCngc9J2gXsBdZX1D8beBMQe3N34nkc2zZ3EzTCpX/NKvlHxC3ALcXjq5tefw747So7ZvPFJP0wiuLYtmpN0A9x3UnXZ/5m1ej/MySz2Zvsu1+0rZz8LbP+G/Yxm7v+j2snf8us/8+QzGYtnPzNOov9MPXTSqqStA74BDAMXBcRf9zy/snA9cAyGmNNb4mI3ZU0bvY81cU11BPbTv6W2STE03OuRdIwsBE4D9gNbJO0KSJ2NhX7L8BnI+IGSecC/xl465wbNzvEVCVxDfXFtpO/ZTYFU7+soqK1wK6IxuRqSTfSWJen+QuyBnh/8Xgr8LUqGjY7REzB1LNV1VZLbHs9f8srpiB+VX7AUkljTceGlpoOrsFT2F281uwe4F8Wj/8FcJSk4+r4WDbfJcZ1xtj2mb9lNgXxXErBPRExMsfGPghcK+ky4FYaSzdMzrFOszamphN7iiyx7eRvmU3B5DNVVDS9Bs+0FcVrB0XEIxRnR5KOBC6OiGwLC9kAiymYrGQ4E2qKbSd/yysCpg5UUdM2YLWkU2h8MdYDlzYXKFZP3BsRUzTW9Lm+iobNDhUQlcQ11BTbHvO3zAJisvwoqyViArgc2EJjV64vR8QOSddIurAodg5wv6QfAscD/6mez2QWMDWVdpTVVFNs+8zf8gqSvgBJVUVsBja3vNa8Ts9NNHbjMqtX0Bj6qaq6GmLbyd8yi0q/JGZ9o6KTmro4+Vt+ff4lMZu9gKn+nkjm5G95BTDp5G8DJoA+D2snf8ssYMobY9kA6vO4dvK3/Lwrog2iPo9rJ3/L6zD4eWw2a4dBXJfO85e0RNL3JN0jaYekj7Qpc5mkxyXdXRzvqKe7NpAiyo8aOLatPsUF35Qjk5Qz/3Hg3Ih4RtJC4DZJ34qI21vKfSkiLq++izbQApjM9vPYsW31CGCqv++hLU3+ERHA9OIrC4ujvwez7PCS6cKYY9tq1ecXfJP+aZI0LOlu4DHg5oi4o02xiyXdK+kmSSvbvG/WXiQcNXFsW21S4jrjvw9JyT8iJiPidBqrya2VdFpLka8DqyLilcDNwA3t6pG0YXrN6scf7+/9La1HopjqWXbU1rxj22owfcE35chkVoNSxRKhW4F1La8/ERHjxdPrgDNn+PvRiBiJiJFly17cTX9tEPXBF8SxbZU73JO/pGWSjikev4DGPpI/aClzYtPTC2msPGeWJt9sH8e21aRYqjzlyCRlts+JwA3FJsJDNJYT/Yaka4CxiNgEvKdYWnSCxs7xl9XVYRsweWf7OLatHiGY6u/bqFJm+9wLnNHm9eblRK+isYGA2exl+unr2LZa9flNXv39T5PND31+G7xZV/o8rvv7LgQbfKmzIhJIWifpfkm7JF3Z5v2TJG2VdFcxdfOCij6F2fMFMKW0I0Edse3kb/lVkPyLcfuNwPnAGuASSWtaiv0hjXH9M2jsg/oX1XwAszYqmudfV2x72Mfyq+bn8VpgV0Q8ACDpRuAiYGdzS8CLisdHA49U0bDZISJgcn9VtdUS207+llcAaWtbLZU01vR8NCJGm54vBx5qer4bOKuljg8D/0vS7wFHAK+fbXfN0ghiUWLZ/Vli28nf8ksb098TESNzbOkS4DMR8V8lvRb4nKTTIryJsNUgPaqyxLaTv+U1vbzD3D0MNK+7s6J4rdnbKe7gjYjvSloCLKWxro9Zdapdz7+W2PYFX8uvmsWvtgGrJZ0iaRGNi16bWsr8FHgdgKR/AiwBHp/7BzBro7qF3WqJbZ/5W34VnCFFxISky4EtwDBwfUTsaLlb9wPAf5f0Phpfu8uKZZ3NKhYwOV5eLKWmmmLbyd/yqnB5h4jYDGxuea35bt2dwNmVNGbWSQhicWLh8vV96ohtJ3/Lz5dbbRAl3sCVi5O/5ZV5Qwuz2vR5XDv5W34+87dBU+1sn1o4+Vt+ff4lMeuKh33MOgj6fvVDs1mLKZh8NncvOnLyt/wmcnfArGpDEC9MLPtUrT2ZiZO/5eULvjaQ1Jju2cec/C0/j/nbIOrzuHbyt/z6/EtiNmuHwS/a0rV9JC2R9D1J90jaIekjbcoslvSlYpeZOyStqqOzNoCmL/iWHTVwbFutqlvbpxYpZ/7jwLkR8YykhcBtkr4VEbc3lXk78GREnCppPfBR4F/V0F8bRPnO/B3bVo+YhP2/zN2LjkqTf7E40DPF04XF0frv1UU0NhMAuAm4VpK8aJaVSt/Mpfqm64ztDFsEREzx5JZv9rzdU84/v+dtPs+OV+Ztv61h0JGJZX9Ra09mkjTmX+whuR04FdgYEXe0FDm400yxAt0+4DhgT0s9G4ANACcdD2xdPafOd2Ptwp43edAfDee5+h9jqTsKVU8jCVvZZRzzryW2T1gEY+vr7vohXvD/HuMFp13V83a3f+C9vPiCN/W8XYDlR72TRUekTqnssT6f7ZO0nn9ETEbE6TQ2EVgr6bRuGouI0YgYiYiRZUd3U4MNpAo2cO9WLbF9jOdRGGlxnfHEZ1abuUTEL4CtFDvGNDm404ykBTQ2EH6iig7agJteAyXzF8SxbZUKimWdE45MUmb7LJN0TPH4BcB5wA9aim0C3lY8fjPwbY/3W7KKZkRIWifp/mJmzpVt3v+4pLuL44eS9jm2rRYxBfufTjsSdBHbpRcSUn6fngjcUIyNDgFfjohvtOwi82kaGwbvAvbS2GbMrFxQyfIORXxupJHAdwPbJG0qNrloNBXxvqbyvwecC2x1bFv1hkAvSizb+Ydkl7F9RlmrKbN97m1XUcsuMs8Bv11Wl1lb1WzgvhbYFREPAEi6kcZMnZ0zlL8E+FBE3Nz6hmPbKlHd78OuYrusUm/gbnml3wizVNJY07GhpaaDs3IKu4vXDiHpZOAU4NvVfRCzZonj/Y0x/yyx7WkJll/aBd09ETFSUYvrgZsiItMdBjbwZreZS5bY9pm/5VfNbJ+Ds3IKK4rX2lkPfLG7zpolqm62Ty2x7TN/y6u67e62AaslnULji7EeuLS1kKSXA8cC362kVbN2YhKe21dVbbXEtpO/5VXRbJ/i7tvLgS3AMHB9ROxombkDjS/OjZ6uafUahuFjE8v+vOO7dcW2k7/lV1EajojNwOaW165uef7haloz6506YtvJ3/Lzev42iLyBu1kH1Y35m/WPwyCunfwtvz7/kpjNnvfwNSvn5G+DZmoCnn0ydy86cvK3vCqa7WPWV7QAFh2XWHimKfv1cvK3vA6DsVGzrviCr1kJJ38bRH1+J4mTv+Xn5G+DJvCZv1lHHvaxQTQ1Cc/094ZvTv6Wly/42iAaWgBLliUW/kmtXZmJk7/l5zN/GzTTe/j2MSd/y8/J3wZRn8d1ygbuKyVtlbRT0g5JV7Qpc06xGfb0BsJXt6vL7BABTCYcNXBsW62mlHZkknLmPwF8ICLulHQUsF3Szc2bBxe+ExFvrL6LNvDynSE5tq0egzDsExGPAo8Wj5+WdB+N/SNn2jzYLF3G2T6ObavN5ATseyx3Lzqa1Zi/pFXAGcAdbd5+raR7gEeAD0bEjjZ/vwHYAHDS8bPtqg2kPpntU2lsn7Covo7a4WF4IXHECYmFd9XalZkk7+Er6UjgK8B7I+KplrfvBE6OiFcB/w34Wrs6ImI0IkYiYmTZ0d122QZNTJYfKSStk3S/pF2SrpyhzO80jfH/VfFatbF9jOdRzHvTN3lVNObfbWx3kpT8JS2k8eX4QkR8tfX9iHgqIp4pHm8GFkpamlK3zXPTwz5z3MBd0jCwETgfWANcImlNS5nVwFXA2RHx68B7HdtWm4qSf7exXVZvymwfAZ8G7ouIj81Q5oSiHJLWFvX29+1t1j8qSP7AWmBXRDwQEfuBG4GLWsr8G2BjREyvtfs4jm2rSyQe5WYd2xFResEh5ffp2cBbge9Lurt47Q+Ak4pGPgW8GXiXpAngV8B6b5BtSaanepZbKmms6floRIw2PV8OPNT0fDdwVksdLwWQ9H9pbIT91zi2rQ4TE/Bk543Zm9QR2x+OiL/t1GjKbJ/bgI6/TSLiWuDasrrM2kq74LsnIkbm2NICYDVwDrACuBU4NiJ+MdMfOLatK8ML4ahfSyy8r5bYlvSKTrGdfMHXrBbV3eT1MLCy6fkKDt0lYzewKSIORMSDwA9pfGHMqlXtBd9aYtvJ3/KrZsx/G7Ba0imSFgHrgU0tZb5G48yI4qLtS4EH5tx/szZiKu1IUEtse06a5RXpUzk7VhMxIelyYAuNMc/rI2KHpGuAsYjYVLz3zyTtpPF74t9FhC/eWj0qWrqhrth28re80i/4llfVmIq5ueW1q5seB/D+4jCrT8XLO9QR207+ll+fr35oNmsTB2DPo7l70ZGTv+UVEPtzd8KsYsOLiKNXJBbOM/Lo5G95eRtHG1SH+6qeZrWrab1+s6z6/KTGyd+yiopm+5j1kwiIjBu1pHDyt7wqnO1j1j/kYR+zMj7zt4FzYD/xaOtNuP3Fyd/y8mwfG0QLFhHHnZRYOHkBuEo5+VteHvO3QeVhH7MOPOZvgyjkC75mZXzmbwPJZ/5mM/NUTxtEsX8/U7sfKi+YkZO/5RUw5Qu+NmgWLiKOPzmxcJ5ZQU7+ll3imuZmh5n+HvZJ2cB9paStknZK2iHpijZlJOnPJe2SdK+k36inuzZwimGfsqMOjm2rT+OCb8qRS8pOXhPAByJiDfAa4N2S1rSUOZ/GlmGrgQ3AJyvtpQ2sAKamyo8UktZJur9I1Fe2ef8ySY9LurvYsP13cGxbTaZQ0pFitrEt6R1ldaZs4P4o8Gjx+GlJ99HYTX5nU7GLgM8WGwrcLukYSScWf2s2s6hm2EfSMLAROI/GfqbbJG2KiJ0tRb8UEZcf0g3HtlUoAqKi2T5zje2ZzGrMX9Iq4Azgjpa3lgPNl7Z3F6897wsiaQONsydOOn42Ldugmj7zr8BaYFdEPAAg6UYaibv1C9JWpbH9j4aIJ34y6w8wV/HEPiZ/2Lq1aw/aXXALPPS1nrcLEL/2M6aeOzZL2x3t38/kTyuLgTnF9kySk7+kI4GvAO+NiKe6aSwiRoFRgJGlivh87+f4veq4njd50D/95qIs7T79hj6eThMwmRYGSyWNNT0fLeJpWrskfVabei6W9JvAD4H3RcRDlcf2qqXBo2/qppo5efFT29Grz+x5u8vPvpOhFx/Z83YBxv/HURxYWjrCUZM/mvmtRYth5UvSqrnrJ7XFdqdmk5K/pIU0vhxfiIivtinyMLCy6fkKcs1fssNKAFORVHRPRIzMsbmvA1+MiHFJ7wRukPTPcWxbxRpb+CYP+9QS28C5nf4gZbaPgE8D90XEx2Yotgn418XMiNcA+zwmaqka46OdjwSlSToinoiI8eLpdcCZOLatJoGSjgTdxnZHKWf+ZwNvBb5fzJAA+APgpKLRT9HYVf4CYBfwLPC7CfWaAcln/mW2AaslnULji7EeuLS5QMuF2gtp/Hx2bFstEhN7im5i+76ySlNm+9xGyd0KxUyId5fVZdYqoprkHxETki4HtgDDwPURsUPSNcBYRGwC3iPpQhrTl/cCF0fED0rqdWxbF8RU0kz6cl3G9mVl9foOX8uuqht8I2IzjTP15teubnp8FXBVRc2ZzSj2j3Pgpw9WV18Nse3kb1kFMFnNsI9Z/1i8mKFVp6aVvetH9fZlBk7+ll1FY/5mfWUWs32ycPK37Lyumw0eOfmbdTKLef5mh40QTDn5m3Xm3G+DJsbH2f/gA7m70ZGTv2UVwISzvw0YLV7M8D9enVZ4+5yW6Omak79l5zF/GzwCD/uYzayqm7zM+kngMX+zUs79NnAEqJo7fOvi5G/ZedjHBk/yom3ZOPlbVoGTvw2eqfHn+NWDf5+7Gx05+Vt2nu1jg0aLl7Do1JemFb79rno7MwMnf8sq8Ji/DSbf4WtWwsM+NnAkwhd8zTrzVE8bRD7zN+vAF3xtUDn5m5WoKvlLWgd8gsZuR9dFxB/PUO5i4Cbg1RExVlHzZgdNjo/zywd2VVZfHbFdmvwlXQ+8EXgsIk5r8/45wP8Epret+WpEXFNWrxkUm7lUUI+kYWAjcB6NvXm3SdoUETtbyh0FXAHcAfwHSa/FsW0VG1qyhBe89OVphW+9vePbXcZ2eR8TynwGWFdS5jsRcXpx+MthszIV5UeCtcCuiHggIvYDNwIXtSn3H4GPAs8B38CxbTWJ4kavsiNBN7FdqjT5R8StNDYENqvc9Jh/2QEslTTWdGxoqWo58FDT893FawdJ+g1gZUR8s3jpLhzbVpMYUtJBPbFdqqox/9dKugd4BPhgROxoV6j4UBsATjqiopbtsJc45r8nIka6bUPSEPAx4LJZ/unsY/s4B7fNaievLLFdRfK/Ezg5Ip6RdAHwNaDtQtYRMQqMAowslSf4GVDZTV4PAyubnq8oXpt2FHAacIsaX8oTgE3AOzvU2V1sr1rq2J7nJsef45m/r2x5h65iW9KFnS76zjn5R8RTTY83S/oLSUsjYs9c67bBV+FmLtuA1ZJOofHFWA9cerCdiH3A0unnkm4BPgjMGKeObevW0JIlvPClL0sr/L9vKSvRVWyXzfaZ8y1okk5Q8c+NpLVFnU/MtV6bPxLH/DuKiAngcmALcB/w5YjYIekaSRd20y/HtnWv2Mwl5ShRR2xD2lTPLwLn0LgosRv4ELCw6NSngDcD75I0AfwKWB8R/tlrSapc2yciNgObW167eoay5zi2rVYVLu8w29hOqbM0+UfEJSXvXwtcm9KYWTu57vB1bFtt5J28zEp5eQcbSENO/mYz8to+Nogmx8d5ald1yzvUwcnf8gqY9Ci6DZihJUs48mWJyzt882/r7cwMnPwtK5/522ASDHk9f7OOnPxtEHlJZ7MOvI2jDSwnf7POfOZvA0cQHvYxm1lV6/mb9ReP+ZuV8pm/DRzh5G9Wxhu42yBK3KglGyd/y8pTPW0wqdK1ferg5G/Z+cTfBpKXdzDrzGf+Nmgmxsd5srrNXGrh5G9ZebaPDaIFS5ZwzMsSN3PJpL8HpWxeqGIzFwBJ6yTdL2mXpCvbvP9vJX1f0t2SbpO0pqKPYHYoDaUdKVXVENtO/pbV9AXfuSZ/ScPARuB8YA1wSZsvwF9FxCsi4nTgT2hsem1WPdEY8085yqqqKbad/C27SDgSrAV2RcQDEbEfuBG46HntNO3JCxyRXrXZbFW3jSM1xbbH/C27xGGdpZKaN6QejYjRpufLgYeanu8GzmqtRNK7gfcDi4BzZ9tXs2TpN3llie2UPXyvB94IPBYRp7V5X8AngAuAZ4HLIuLOsnrNYFbz/PdExMic24vYCGyUdCmwWdLROLatYhPj4+xNn+1TR2z/IfC2TuVTzvw/Q2Mf08/O8P75wOriOAv4JG3+VTKbSUWbuTwMrGx6vqJ4bSY3AqPAb+LYtooNL1nCsdXN9ukmtj9ZVmnp75KIuBXY26HIRcBno+F24BhJJ5bVawbVXfAFtgGrJZ0iaRGwHtjUXEDS6qanbwB+gGPb6lLdbJ9uYvtHZZVWMebfbjxqOfBoa0FJG4ANACcdUUHLNhCqOPGPiAlJlwNbgGHg+ojYIekaYCwiNgGXS3o9cAB4kpKfxXQb2y9aAD/6zlw/0qwd+NnP+dUjj/S83WP3H2BqYqLn7QJMPfsrJr57W5a2O5qe7VOBmmK7txd8i4sYowAjyxRa2MvWG67s/XfjoPvP2J+l3c++KEuzyaq6wzciNgObW167uunxFa1/I2lVRW3/Q2yfsjT0srOrqHZW9i7aA298V8/bfeSxi1lx7Kt63i7A5LF/zgtPzhTgf3NLhzerXdunm9guU0Xyn+14lNlBfb6wm2PbujI5Ps7eBx/M3Y2Oqkj+0z85bqRxMWxfRBzys9hsJn28vINj27qyYMkSjnvp6vKCGaVM9fwicA6Nuai7gQ8BCwEi4lM0fopcAOyiMR3ud+vqrA2enGf+jm2r1eG+pHNEXFLyfgDvrqxHNu/kus3WsW218k5eZp318Zi/WXeUtm5PTk7+llWfX/A1697hPuxjVjevrmaDZmL/fvb+5Me5u9GRk79lFUCe24PM6rNg8WKOO/XU3N3oyMnfsguf+tsg8rCPWWce87eBIyHP9jGb2Sw2azE7rChto5ZsnPwtO5/520Dymb9ZZ07+Nmgm9u/niYceKi+YkZO/ZRX09do+Zl1ZsHgxS1/yktzd6MjJ37LzmL8NGkl9P+bf34NSNvAq3MkLSesk3S9pl6Qr27z/fkk7Jd0r6f9IOrmaT2HWxtBQ2pGgjth28rfsIuEoI2kY2Ehj3901wCWS1rQUuwsYiYhXAjcBf1LJBzBrQ0NDSUdpPTXFtod9LLuKLviuBXZFxAMAxRr8FwE7pwtExNam8rcDb6mmabPnmzhwgL0PV7bvTy2x7eRvWVV4wbfdfrtndSj/duBb1TRt9nwLFi1i6apVVVVXS2w7+Vt2iWf+SyWNNT0fLfbNnTVJbwFGgN/q5u/NSkmNI02W2Hbyt+wSZ/vsiYiRDu8n7bcr6fXAvwd+KyLG03tpNjuzWN4hS2w7+VtWFa7nvw1YLekUGl+M9cClzQUknQH8JbAuIh6rplmzQ4lZJf8ytcR2Uu8SphldJulxSXcXxztS6jWDamb7RMQEcDmwBbgP+HJE7JB0jaQLi2J/ChwJ/HURp5sc21aL6Z28Uo4S3cZ2Wb0pG7hPTzM6j8aFhm2SNkXEzpaiX4qIy0s/iVmTKnfyiojNNDZdb37t6qbHr29+r4jtH+LYtopNHDjAL37288rqm21sp0gZ9imdZmQ2FxmXd3BsWy0WLFrEcStXlhfMKGXYp900o+Vtyl1c3F12k6T+/tTWV6oY9umSY9tqoyElHblUdUXi68Cq4u6ym4Eb2hWStEHSmKSxx5+rqGU7rFW5vENNZh/bTzu45zsB0lDSkUtKy6XTjCLiiaapRdcBZ7arKCJGI2IkIkaWLemmuzaIMp751xPbRzm4571iJ68qlneoS0rLB6cZSVpEY5rR864kSzqx6emFNK5ImyXJeObv2Lba9HvyL73gGxETkqanGQ0D109PMwLGImIT8J5iytEEsBe4rMY+2wDJuY2jY9vqMjkxwb49e3J3o6Okm7wSphldBVxVbddsvsi5mYtj2+qwYOFCXnziieUFM/IdvpZVAFPezcUGTTHm38+c/C07534bRE7+ZiW8gbsNoqE+38bRyd+y85m/DZqpyUme3rcvdzc6cvK3rCrczMWsbwwvWMAxy5bl7kZHTv6W3VTKr2P/PLDDiS/4mnWWc56/WZ2c/M1KeKqnDRpJDDn5m3Xm2T42kPp8tk9//9NkAy9lUbfUHwYJu3L9pqQ7JU1IenM1n8DsUFNTU4w/+2zSkaKO2PaZv2U3UUEdiTvO/ZTG2jwfrKBJsxkNDw9z5DHHVFJXXbHt5G9ZVbiNY+muXBHx4+I9jzRZ7Sq84FtLbDv5W3aJ0bpU0ljT89GIGG163m5XrrPm3Dmzbszugm+W2Hbyt6xmcea/JyJGau2MWUXErM78s8S2k79lV9EYTOmuXGY9U+1Uz1pi28nfsqpweYeDu3LR+GKsBy6tpmqz2YkIJiaqmMoA1BTbTv6WVQAHqqgnYVcuSa8G/gY4FniTpI9ExK9X0LzZ8wwNDbHkhS+spK66YtvJ37KramG3hF25ttH4yWxWuyrv8K0jtp38Laugmnn+Zv3kcFjeIal3CXeXLZb0peL9OyStqrqjNpimk3/ZURfHttVlaGgo6cjWv7ICTXeXnQ+sAS6RtKal2NuBJyPiVODjwEer7qgNpukx/7KjDo5tq8v0mf9hnfxpurssIvYD03eXNbsIuKF4fBPwOqnPVzWyvjAFPJ1w1MSxbbUZHh5OOnJJGfNPubvsYJniyvQ+4DhgT3MhSRuADcXTcf0lf9dNpyuwlJa+DXK7y/flaxt4Wcn7W2j0rUwdfa8vti/7TI7YXspHvpHj/3Gj7Xn0nSrMGNvbt2/fIiklriFT/3t6wbe4ZXkUQNJYrjs2c7U9Xz9zp/cjYl2v+lKnfojt+RpfOT/zTO8dDnGdMuyTcnfZwTKSFgBHA09U0UGzGjm2bd5KSf4H7y6TtIjG3WWbWspsAt5WPH4z8O2I8P5M1u8c2zZvlQ77pNxdBnwa+JwjGb5vAAACOUlEQVSkXcBeGl+iMqPlRWqTq21/5j4ygLHt+Jo/bc+ZfBJjZjb/9PctaGZmVgsnfzOzeShL8i+7pb7Gdq+X9Jikns7BlrRS0lZJOyXtkHRFj9pdIul7ku4p2v1IL9ptan9Y0l2SvtHLdnNxXPcmrou2Hdtz1PPkn3hLfV0+A+SYfzsBfCAi1gCvAd7do888DpwbEa8CTgfWSXpND9qddgVwXw/by8Zx3dO4Bsf2nOU480+5pb4WEXErjRkbPRURj0bEncXjp2kEzfIetBsR8UzxdGFx9OQKv6QVwBuA63rRXh9wXPcorov2HNtzlCP5t7ulvicB0w+KVSHPAO7oUXvDku4GHgNujoietAv8GfD7VLZLY99zXPcwros2Hdtz4Au+PSTpSOArwHsj4qletBkRkxFxOo27V9dKOq3uNiW9EXgsIrbX3ZbllyOuwbE9VzmS/7zcaFvSQhpfkC9ExFd73X5E/ALYSm/Ghs8GLpT0YxrDH+dK+nwP2s3JcZ0hrsGx3a0cyT/llvqBUiwB/Gngvoj4WA/bXSbpmOLxC4DzgB/U3W5EXBURKyJiFY3/v9+OiLfU3W5mjuvetu3YnqOeJ/+ImACmb6m/D/hyROzoRduSvgh8F3iZpN2S3t6LdmmcLbyVxlnC3cVxQQ/aPRHYKuleGsnp5og4bKem9TPHdU/jGhzbc+blHczM5iFf8DUzm4ec/M3M5iEnfzOzecjJ38xsHnLyNzObh5z8zczmISd/M7N56P8DwBENUOzGgAAAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "# Choose colormap\n", "cmap = colorcet.cm.fire\n", @@ -332,7 +586,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 16, "metadata": { "scrolled": false }, @@ -343,10 +597,10 @@ "010000" ], "text/plain": [ - "" + "" ] }, - "execution_count": 22, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -362,7 +616,7 @@ }, { "cell_type": "code", - "execution_count": 39, + "execution_count": 17, "metadata": {}, "outputs": [ { @@ -371,10 +625,10 @@ "010000" ], "text/plain": [ - "" + "" ] }, - "execution_count": 39, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -395,23 +649,21 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 18, "metadata": { "scrolled": false }, "outputs": [ { - "data": { - "text/html": [ - "0.01.0" - ], - "text/plain": [ - "" - ] - }, - "execution_count": 32, - "metadata": {}, - "output_type": "execute_result" + "ename": "NameError", + "evalue": "name 'cm' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mcm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mStepColormap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'#64abb0'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;31mNameError\u001b[0m: name 'cm' is not defined" + ] } ], "source": [ diff --git a/pcolormesh.ipynb b/pcolormesh.ipynb new file mode 100644 index 0000000..d443901 --- /dev/null +++ b/pcolormesh.ipynb @@ -0,0 +1,126 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEYCAYAAACdnstHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztvXncJVV17/1dPQMNNEjL1NCNYVDEZjIOQRNETQAVEzWKU8RXwkuiN/gmjjEfMSZcNa9643Sv1yARJ4SgcpE4gFEBUZFmnhQRaLuZh26ggW766V73j6rqrrOfqtp7V+0aznn29/PpTz/nnDpV++yq2r9aa+29lqgqkUgkEonUYVbfDYhEIpHI+BJFJBKJRCK1iSISiUQikdpEEYlEIpFIbaKIRCKRSKQ2UUQikUgkUpsoIpERRORDIvLVvttRhYicICI/bWG/60TkaRWf3yEiL3HcVyttNI6xTERUROa0eZxIpIp48UUiKaq6MPtbRL4ErFbVf+ivRZHI8ImWSCQI8Wk4EpmZRBEZY1L3yvtF5CYRWSMi/y4iC3Kfv1JErhGRR0TktyJydPr+HiJyvog8JCK3ishfVhzjOBG5UUTWishPROQZxvHfKyLXAY+JyJz0vXeLyHUi8piIfFFEdhWR74nIoyLyQxHZKbeP54nIz9L9XysiR+Y+O0FEbku/d7uIvNFo28fT3327iBxT0v63ish3cq9/IyL/kXu9SkQOSf9WEdlXRE4C3gi8J3VxfSe3y0PS3/awiJyd7+8qROTpInJR2ue/FpHXpu8/V0TuEZHZuW3/LO1TRGSWiLwvPX8Pisg5IrJzyTEq+ysSaQVVjf/G9B9wB3ADsBewM3AZ8M/pZ88BHgZeSvKwsCfw9PSzS4D/CSwADgHuB45KP/sQ8NX07/2Bx9J9zAXeA9wKzMsd/5r0+Nvk3vsFsGt6zPuAq4BD0+P9CDg13XZP4EHg2LSNL01fLwa2Ax4BDki33R14Zvr3CcBG4C+B2cBfAXcBUtBHTwPWpvvfA1hJ4qbKPlsDzEpfK7Bv+veXsr40+vuX6X52Bm4GTi45NycAP03/3g5YBbyVxIV8KPAAcGD6+W+Bl+a++x/A+9K/T0n7cwkwH/jfwFnpZ8vSNs+p6q/4L/5r81+0RMafz6rqKlV9CDgNeH36/tuAM1T1IlXdrKp3quqvRGQv4Ajgvaq6XlWvAU4H/qJg368D/jPdx0bg48A2wB/ktvl0evwncu99RlXvVdU7gUuBy1X1alVdD3ybZBAFeBPwXVX9btrGi4AVJKICsBk4SES2UdW7VfXG3DFWquq/qeom4EySQXNX8weo6m3AoyRi+YfAD4C7ROTpwB8Bl6rq5qoONvi0qt6V9vd30v3aeDlwh6r+u6pOqerVwDeBP08/P4v0vInI9unvPyv97GTgA6q6WlU3kIj8a0rch1X9FYm0QhSR8WdV7u+VJE/JkFgHvy3Yfg/gIVV91PjeniXbrsxepIPtKmPbVeaXgHtzfz9R8DoLYC8F/jx1Za0VkbXAC4DdVfUxEhE7GbhbRP4zHfgz7sm16/H0z4UUczFwJImIXAz8hERA/ih97cM9ub8frzhmnqXAc43f+UZgt/TzrwOvEpH5wKuAq1R1Ze67385972ZgE4ZgOvRXJNIKUUTGn71yf+9N4taBZHD/vYLt7wJ2Tp9489+7s2TbpdkLEZH0ePltm6SBXgV8RVUX5f5tp6ofBVDVH6jqS0msjF8B/1bzOJmIvDD9+2LsIhIyvfUq4GLjdy5U1b8CUNWbSMT6GOANJKKS/+4xxncXpFbeaIPD9Vck4kwUkfHn7SKyJA22fgA4O33/i8BbReTFaXB2TxF5uqquAn4GfEREFojIchLXV9HakHOAl6X7mAv8HbAh/X4Ivgq8QkT+RERmp+05Mv09u6YTA7ZLj7mOxF1Th4uBF5HEbVaTuNiOBp4CXF3ynXtJYiYhuADYX0TeLCJz03+/n5+kQCIcp5BYS/+Re//zwGkishRARBaLyCvNAwTur0jEmSgi48/XgQuB20jcV/8MoKq/JAnk/g+SAPvFbLUqXk8SlL2LJEZxqqr+0Nyxqv6aJG7xGZJA8CuAV6jqkyEangraK4G/JwnurwLeTXJdzgL+Nm3jQyRWw1/VPM4tJIPqpenrR0j667I0plLEF4EDUzfSeXWOmzv+o8AfA8eT/J57gI+RBMozziL5jT9S1Qdy738KOB+4UEQeJQmyP7fgMMH6KxLxQVRjUapxRUTuAE4sEoBIJBLpgmiJRCKRSKQ2TiKSLiC7XpKFayvablQkEolE/BCRM0TkPhG5oeTzN6YLZa9PF/geHOS4Lu6s1G3ybMNXG4lEIpGBICJ/SBL7+7KqHlTw+R8AN6vqmjTDw4dUtSi+5kXMdxSJRCITgKpeIiLLKj7Pz6rMsiA0xlVElGR2iAL/W1W/YG6Q5hs6CWC7bbY9/IBl+xbvaa7FgzZn9PP1mx4Zeb1h02Ojn298fPTzJ580XpdNvqnH/HmzjdfzRl4vmLvt6Oezt6vc39xZTqmXnNi4ef3Ia1tfmZh9Z2L+VhPbbzd/6+xNxuW3cXRG6sa1o+2dWrdh9PPHR1+v3+Q3aWzB7Orfs81Tdxh5PXub0e1nLTS+H/ja9cH3Opx2LsT4LVOW2cHGufK9rzcZE/x8r92h3ee/veHuB1R1cajjb7fXAt20frSPNzyw8UYg31FfKBqLHXkb8L2a3x3BVUReoKp3ishTgYtE5Feqekl+g/THfAHg8AMP1p997fuFO5q9W/WgOrXT6Mn67eOjE49uXfPLkde/umt0mv+tq0YXUN++cm3l8XzZZ+mikdf77rXXyOun73Ho6Oc7Padyf7tu84zKz32494mbR17b+srE7DsT87ea2H67+VsXrR295zbdMzpw3HPelSOv11x628jru667feT1rx/6XWX7TA5YtHfl58/8yxePvN7hWUtHXm9/xOjr0NeuD77XoXkuFspowoI5a6oFzTxXvvf1OmOtpO+1O7T7/E/3/aeVBGTT+s0s+9NdRt779el3r1fVZzfdt4i8iEREXtB0X+AoItnqWFW9T0S+TZLc75Lqb0WKqBIN80Y2MW88GzYBM29M80Yxb1Tztbm9TaTM324OLObFaA7aJqaIHLDzqCjYRMX83Py+KVom2/7e6E1utn/XRc0eEHxEpYkAAUlGtBwLd6oWFZtodI056IcWlUkhXVx8OkkWhAdD7NMqIukK2Fmq+mj69x8DHw5xcBfMgcd8QplJmCJjiorZV+bTXWhsomIyrT3GwLWIUcvEHKQfuX70YW+P5fuMvDZFxZdponNd9famyG1rfG4OxObvtV3L+f5sKhI2gZ+GRVT6xvbAE5mOiOwNfAt4c7oANwgulsiuJAngsu2/rqrFvqoCfM3c0EzyE4pNVEx8LRNfmt7I055+jc9tlsk0DBHwdXeZmCK10/XVWVG2N6598/f5nA/fQdNX4K0PZ6aoWKzmprR9rU4iInIWSY64XURkNXAqSQkHVPXzwAdJUv38z3Q8nwrhHrOKSJpKO8h84kgkEom0g6q+3vL5icCJoY/bzhTfubMG5zMdKiGf6Kbty3h69HVvhXYZTIvBGE+bpiVlPrmb7i0TWwzDl7ZjJr+300uc2+Ibv/LdnxkoNrG5IkPTt9t6kjwWbTN260TMgWeSzVrbDBkTX9dg1+6taQOfpwtl0W5GzMTYfKcXWpLuNnRvhY6ZNAnEty0qJua1Ynsgsc06bDteN+7Mmy8sedro/fzrntpiY+xExMT2BDWTA26hLZOm+J4L22wuMaaZtj2by4ZvzMQ3EF+FTVR88X2A8BUZXyb5YXHcGXsRmSR8LQ/b90NbJiZtC7R1NpfFMjFpezaXr7vLpDIQ7+k+Cj3o2txf1plmxrXVt7vKl+jeKieKyBhhLvAyMeNQpqjYnnSHbpmY2GZz8dvRVG9Dc3eZVMVQmlgp0L2omPiKRtOFsTMRETmapP7MbOD0rEJo7vOlwBnAYpKaM29Ki7Q1IqaCj0QikTFHRGYDnyMpsXwg8HoROdDY7OMkyRmXk6z1+0iIY0dLpCHWVCGeLqIm+KaisMVMfGn6dNg0ZuI7m8tKy+tMTKoC8eaNav7We6m2Ik1LoW3LZNIZoHvrOcCt6ZIMROQbJFVDb8ptcyBJ9UuAHwONKnZmtCMic2Zt8cf7+ul9U3tEyjFFxTzZbS/0NDFvPHMhqC+hYya2GEbotCom5or8PE1ndpnupLZFJTTRfQUkiwjz9ZzyCRj3JCkvnbGa6WWUrwVeReLy+jNgexF5StP0J61bIl0PVEPDvHl9Ei7aYiBN8Y2Z+GLe+LbsAebnbcdMfKcId51WpYqmM7uGJipRJEaZP292wUPWygcarjB/F/BZETmBJPfhnUDj9MfRneVJ06dp8+l54bb1Fxs+/tvqGmGmO8fX3dU1oUXF5u4aeloVH2wpVroWla5FYYDupa65E8jP616SvrcFVb2LxBJBRBYCr1bVxh0XRWSCMUXGJio2v3vXMZOmojItHmVJKth1DMXH3WW62kzBC30uh2apmETRmMYVwH4isg+JeBwPvCG/gYjsAjykqpuB95PM1GrM4EVk0le25mNAZhbb0NhEpWt8k2P6WoGhV8SbhI6hmOQtE3ONS1X8BOwxlKai0jZRJPxQ1SkReQfwA5Ipvmeo6o0i8mFghaqeT5Kc8SNpccFLgLeHOPbgRSQSiUQidlT1u8B3jfc+mPv7XODc0MeNIjIgbEWaQgfap1kmxudDc2+ZhH5abTsQb5ud5bMvX9o+t01zeUXLY3xpRUQ26ZNb3DShq/VF3LG5PEy/elNRsa1VsBG69svQAvFNhaAK3xQrtnOL5+zrocVMxp358+YV5D+7tpe22HAWkXRF5ArgTlV9uev3okh0h000bNt7z04yaFovPrRlMvRAfBNRMb9rxkxMzHNta7vvQlSfKo0wflN67bMwg5ZYHyt8LJFTgJuBHVpqy0QwpBXsvnQdrPXNPOtrmQwtEG8jpKiYmK63rl2ZNstkaBVImy6E7QNb7qx0m9cCHwIUuFZV32Bu44uTiIjIEuBlwGlsXTbfC+OW/dOkyeLDvvEdeGwuEV8XSNOBJvjAtMh8ORxRMQnt7mrqyhx3y2RorqZc7qyXkqxWv0JEzlfVm3Lb7EcytfcIVV0jIk8NcWzXBIz/CrwH2Fy2gYicJCIrRGTFgw88FKJtkUgkEnFjS+4sVX0SyHJn5flL4HOqugZAVe8LcWCrJSIiLwfuU9UrReTIsu3SHC5fADjksGdpiMaNAyFXsDdZvV4H8+nUdHk0jZkMLcFj2/v3nd3VhNBBepsr01wRb+JrUYcOvK++bTSFj1kV0Ibtvm1a5KsDXHJn7Q8gIpeRuLw+pKrfb3pgF3fWEcBxInIssADYQUS+qqpvKvvCxs3rtwyOvhfXpC8urMKchOC7+NAc5M2BwebS8BUVk6Z+dV9XZWj3VtdpVnxExTw3tkC67/5s05VDr4gPHQ/zFY2hs2DutgU1Wi6oSsDowhxgP5JFh0uAS0TkWU1Tn1hFRFXfT+JHI7VE3lUlIJFw2NaNmIOQLZeWL22LihlTaFpZsW1RMZk2EHrO7gopKqHxPbd4Jlod9ynBBQN8F4etSsBozZ1FYp1crqobgdtF5BYSUbmiSaNaX2xoWhbm09pMtjzGDd+nV1PUpiUJHLi7yyZC04LDnrO7xklUbDTNCD1us7cGiDV3Fkn9kNcD/57m0dofaHzheImIqv4E+EmTA4YWjaE/wfgMbG3P1DIH/aYDj+37bScJ9KXpwBN84PIQlabxlL5Fxoa1vO6Yzd7qGsfcWT8A/lhEbiJJAf/uprVEIKY9iUQikYnAIXeWkizRCLpMI4rIgLBW6ms5y29T2l6L0Hcg3kbjp+UKy8TX1WXGNGyuR18mLUVOpD6DF5FxW1wYugRsHt9Au4lvWhTbNFLfGUIzLRAfUlRCTx/2vRZC03eKnKEzf/Z2g85ukacVEdmw6bEtg/+4dETEH1v+ptCzu/ou32ujVVFpGJQ3LYemedImjTYf/iadwVki42Z5NMVnYkDoxYihA+1NRcXG0GrC2xiyqJiEnh7eNdG91R+ti4gpCjPdMjEHDp+VsOZixKYDh80SMEXAd5W0r6jY6LveSd/rUCrxFBWTkKvpJ5FxcJ/ZEjCKyMkk1Qw3AeuAk/K5terSuSXStqUxDid7phJFxW//+QcM76nsNlExrLrZxhoe3wJo5vfbpu0EjuPm3nJJwAh8XVU/n25/HPBJ4Oimx3ZNwBiJRCKR4WJNwKiqj+RebkeSDr4xg4uJ2Bj64kKTNp9oQru3TGz5mpomAZzplolPzMTm9my6Wt6k74HB5vYet3EgEFW5s1wSMCIibydZJzIPOCpEo/q+ViI5fG+McapF4sJMFxWffZuYK759RcXEJjJtM/QJNm27zefOWlB0f1flznJCVT8HfE5E3gD8A/CWJvuDKCK94/O0GRrf8q5NLZNfP/S7ys8P2Hnvyv1NmqiE3LeJr6iYA5atrPW0vGcWxr1M9hjESFwSMOb5BvC/Qhx48CIyQ83WQkIn+PMNnoaeEmxiisyki0qbmPeNt6Vi4CsyNnxz6MVxwIo1AaOI7Keqv0lfvgz4DQFoRUTWb3x8y0m3JVYziRfLVmwDQVPMGTVtTwk2RcFmmUy6qJg0ERnb1HFfUTFjEm1n2x66+8pkaOtQHBMwvkNEXgJsBNYQwJUFY2CJRCKRSMSOQwLGU9o4busiYnsCamp5xHUhW7ElcLStFbBhS5XRdrrxSbNMbNX9hmSZmDRdNNzU8oj3/XBwqbG+ALgEmJ9uf66qnlr3gNFdVY7vanbfbABNpwT7psYI7d4ymWmikif0IGq7L5vGVHyZ6ePEbJnnPXmhL1wskQ3AUaq6TkTmAj8Vke+p6i9ablvEwDdGYqsqacM3ZmLiO5sriko5vjXJQ+MrMk33ZxItj+HiUmNdSfKsAMxN/wVZ6VgH34tpaAGwNglumVjcXU2nCLfNTBKVvgndnnifjw9OMZE0L8uVwL7A51T18oJtTgJOAli0a+BpK5FCmvq1bdM2254ibBLaMjGxfv+60ZfjJComQxMZG9HSaI5DAsb5wJeBw4EHgdep6h1Nj+skIqq6CThERBYB3xaRg1T1BmObLwBfAFjy9EW9WSqTRJOMvy74urtCF8XqW1RMpu2vb1HJYxEY84HB5l5qW2TaFoVoeYzimIDxbcAaVd1XRI4HPga8rumxvRIwqupa4McEyPwYiUQikWBYEzCmr89M/z4XeLGISNMDu8zOWgxsVNW1IrINidJ9rOmBXYlmrjvBF5SZU4SN2SJtl+ttOxBvo2/LJE9T15fNUhma+2volsb0tCeByw1PbS6KSTZNwLhlm3Rx4sPAU4BGFclc3Fm7A2em5tIs4BxVvaDJQdtk6Bdfl4ybqJh0PbvLhk1UmlLVP6HjKUNzfw2dgeTKapyAsQ1cZmddB4TNt1HBpFkeA7n4nAgtKibjNmXYhu/sryY0iqfA4C2XoT38jdN9m+KSgDHbZrWIzAF2JAmwN6L3tCeTJhpDxnedSWNR8aymF0XFHR9XGNgtF/NcmufaV1TarjzYN9MnuVzbSztyWBMwAueT5Mv6OfAa4EfpEo5G9C4ikXLaXlDWdJVyFJVqmopKfvaaLX5k66uu3WG+ohK61LAvNsuj68WdvjgmYPwi8BURuRV4iERoGhNFJBKJRCYAhwSM64E/D33czkUktFkbfant0XpgPlomle1pgplNwFwYarNMzIWntjVEfc/+Wn3b6LWz5GnzSraMhCZaIgNi6CZz16Ji0nS2lw2bqJiEFhmTfHuaTic28Q3UNxUVk9Axk65FY3o8MfCE1Y2bvTNC9EUUkYY0tTyGLhxVBH+67Hhdim3FvEnXCSOrju0rKmYafzMjc9uiYmK7dvqOkUTcGTsRGfeLaZxFw5fGJVotwV7fVPa+NHWHmYS0XIYuKjasMwMnbDZX34jIzsDZwDLgDuC1qrrG2GYp8G2S9YBzgc+o6udt+x47EYm0R9Mb1yaQoet++9aUD+3+suErMiGP1beo+Ob6shEtk8a8D/gvVf2oiLwvff1eY5u7geer6gYRWQjckObfuqtqx165syKRSCQyluTzZp0J/Km5gao+qaob0pfzcdSHaIm0TBP31bilorBZMrYSrTYaz/6y7L9pTMXX/TUkzN9mK4VsYs6sC53mPrq3GrOrqt6d/n0PsGvRRiKyF/CfJGU/3m2zQiCKiDdtTuH1rQ4XWmTaXmXctshMY2Czv2x0KTKm4JmC6Lu40ZyO3Tah3Vvm9/ueqr95w1RROeqqBIyIyA+B3Qp294H8C1VVESlcqa6qq4DlIrIHcJ6InKuq91a1NYpIYHwtD1/haLLvoVkuJrb6KX0H6rc3BspxmYLZBb7rUO7FsCItTHoaFUcqEzCq6kvKPhORe0Vkd1W9W0R2B+6rOpCq3iUiNwAvJEkbX8rgRWTSAmhtiobvsW2i0veN2/lA0bHl4jPluOvSwjamBd49SyWHnhLcNhMgWlnerI+m//8fcwMRWQI8qKpPiMhOwAuA/2Hb8eBFZCZjq5Fuw/Zk3rWlYnsg8HUh+Fou3jS0XEyaiIxtTcvQRKZrJt29FYCPAueIyNtIip+8FkBEng2crKonAs8APpG6ugT4uKpeb9uxS1GqvUjq8u4KKIkf7lN1f0mknKaiYdtfU1Fp2zLxvfHNG3vSLZc8TRdORpoxbpaJqj4IvLjg/RXAienfFwHLffftYolMAX+nqleJyPbAlSJykVG7NxKJRCIzEJeiVHeTLEJBVR8VkZtJyizWEpG+/exdY4uBhLY+fI7lvbjP3J/nuezb5TA09xdmGhdjmmzI2WG22VYmvlN8u6brcaRr99amJ570Pmd94RUTEZFlJFUOLy/47CTgJIBFuzacFD5GtJnGxDf4aDJtHYVB08JDQ8M35tK3+8vEJ8bS9fTjacf3DKT7Ynu4Gvq1OJNwFpF0Gfw3gXeq6iPm5+l85S8ALHn6osbVsiaVqpujqWj47q/vanZdp7JoKjLBLRcTjxiL7cY1p9y2PeiHXifS1Epum0mbNdoEJxERkbkkAvI1Vf1Wu02KuGKraW5LimeKzEwTFROby2JIlotvEN8XU4R8RcLMuBzpF5cEjOl2ewOnk9RiV+BYVb2jat8us7OEpKzizar6SZcGL5i77ZYBpunTWt8DS9u+zyprwSYSNnxFJorKKG1bLo0SUPrGXyzYCoLZsImGb5Zfk+i+aoxLAkZIZuKepqoXpd6nzbYduyTYOgJ4M3CUiFyT/jvWo/GRSCQS6RdrAkYRORCYk071RVXXqerjth27zM76KcnCk4gDfa5I98W0VJoWGgqN+WTft2ViMjj3Vx7PuS2+losNm+Vhm/QRKaQyd5YFlwSM+wNrReRbwD7AD4H3qeqmqh33vmJ9pk35rcLXfWW6IHwxXRDj5v6yuRonzf0VEvPcNXU32bCJhm8gfdLHial1G4qyEFTmzgqQgHEOSa6sQ4HfkcRQTiAJZ5TSuoiM2zTR0IRcB9JUNHz3ZxOZoYmKySRZLq3GV3rA1r6ZNk6EIEACxtXANap6W/qd84Dn0beIROrjKxq+WWVtwVPb8a0ukJ6XC42b5VLVvqZWi6/IhF4E21Q0bA8MfT8gjAHWBIzAFcAiEVmsqvcDRwErCrYbIYpIx4SMMzRNRd5UdEyRGbrlYuIrMiZ9Dlw2UTEJLTJNGXdLY/q1MvjV5dYEjKq6SUTeBfxXOiv3SuDfbDuOItIzPnGQpqJRUORmBNuCNNvagaFbLk1FxmRImV19RcXEV2Sa7s/G0CyPIZ3rOrgkYExfeydhjDXWI5FIJFKbViyR+bO32+KKGPd8TOOEzdJo+n3TUrFZRjPdUumTppaJSej71LfvZlrMY+PjGzotl9yEwbmzfKd1TtLFFXr2VWh8RcqWJHDcRKYpQxadthm3325zX00X5Wvba8zAGZyIRNzxHdRDp5a2pQu3WjaW/Q9NZGxWtYmvZTNJNBWNSXo4nHRaF5HQNSwi7rRdj6BpjYqmImNiu5jHzZKZSa7evkXD3/IYLzwSMH4MeFn68p9U9Wzbvnu3RGKMpJyms7F8aVqn21ai1SY63iLjGaMZusj4Wjo+xPsqLNPPzQW9tMMDawJGEXkZcBhwCDAf+ImIfK+o9EeeODsrEolEJh9rAkbgQOASVZ1S1ceA64CjbTvu3RKJhMPXvdTU8mi6P9Nyadv9ZnOP2W6GRSweeW2tmdGhe8zmJm5q1fhaMl3nxFt922i8bMnTxrueyfpNT/Lrh35nvt12AsZrgVNF5BMkt8uLcCiD3oqIzJ21YMtq5KbZOvtO0Nh1beU2CS0aTRma6HQegzHpMSbTVIT6rhs07qLhSKsJGFX1QhH5feBnwP3Az4HKDL7QgyUSOtA+yVN+zdlJXcdITNqet77H8n0qPx830RHP6oCN0rF7CpAtxYyNtpNnRvwJkIARVT0NOC39zteBW2zHdalseAbwcuA+VT3Itn1kcuh6sVPT45kiNDTLK7To5Gnb6uk7zf8kPRz2hDUBo4jMBhap6oMispwk/cmFth27WCJfAj5LUjbRG1uSPZNJm61ltbRyN7M5EPguPjRnN7UdYxgaviI07qKTx9elYBOdrpNlRlrHmoARmAtcmuRe5BHgTao6ZduxS2XDS0RkWe2mRyKRSKRXXBIwqup6khlaXohqUYErY6NERC6ocmeJyEnASenLg4AbfBvTEbsAzZJMtUtsX32G3DaI7WvCkNsGcICqbh9qZyLyfZLfnOcBVbVOue2aYCJibL+iahZBnwy5bRDb14Qhtw1i+5ow5LbB8NvXJnGxYSQSiURqE0UkEolEIrWxioiInEWy6OQAEVmdRvdtuK6i7IMhtw1i+5ow5LZBbF8Thtw2GH77WsMpJhKJRCKRSBHRnRWJRCKR2kQRiUQikUhtaouIiJwhIveJSOF6EEn4tIjcKiLXichh9ZsZvG1HisjDInJN+u+DXbUtPf5eIvJjEblJRG4UkVMKtuml/xzb1lv/icgCEfmliFybtu8fC7aZLyJnp33zuPUBAAAgAElEQVR3eZeLZR3bd4KI3J/rvxO7al96/NkicrWITCuC0WffObav7767Q0SuT4+9ouDz3sa93lDVWv+APyQpYHJDyefHAt8DBHgecHndY7XQtiNJ1r100p6C4+8OHJb+vT1JkrMDh9B/jm3rrf/S/liY/j0XuBx4nrHNXwOfT/8+Hjh7YO07Afhsj9ff3wJfLzqHffadY/v67rs7gF0qPu9t3OvrX21LRFUvAR6q2OSVwJc14RfAojR7ZOs4tK1XVPVuVb0q/ftR4GbATGbUS/85tq030v5Yl76cm/4zZ4fkC/CcC7xY0oRAA2lfb4jIEpLyp6eXbNJb34FT+4ZOb+NeX7QZE9kTyKfmXM2ABiPg+anL4Xsi8sy+GpG6Cw4leWLN03v/VbQNeuy/1N1xDUk664tUtbTvNEkg9zDwlAG1D+DVqbvjXBHpsoD3vwLvATaXfN5r32FvH/TXd5A8EFwoIldKkurJpPf7tmtmamD9KmCpqh4MfAY4r49GiMhC4JvAO9VSx7hrLG3rtf9UdZOqHgIsAZ4jIoMqUeDQvu8Ay1R1OXARW5/8W0VEspIOV3ZxPF8c29dL3+V4gaoeBhwDvF1E/rDj4w+ONkXkTiD/lLAkfa93VPWRzOWgqt8F5oqImeysVURkLskg/TVV/VbBJr31n61tQ+i/9NhrgR8zvQ70lr4TkTnAjsCD3bauvH2q+qCqbkhfng4c3lGTjgCOE5E7gG8AR4nIV41t+uw7a/t67Lvs+Hem/98HfBswq3kNdtxrizZF5HzgL9LZCs8DHtatNX57RUR2y/y8IvIckn7obJBJj/1F4GZV/WTJZr30n0vb+uw/EVksIovSv7cBXgr8ytgsK8AD8BrgR6raSVzCpX2Gj/w4krhT66jq+1V1iaouIwma/0hV32Rs1lvfubSvr75Lj72diGyf/Q38MdOzlQ923GuL2uVxJUmHciRJ8fjVwKkkQURU9fPAd0lmKtwKPA68tWljA7btNcBficgU8ARwfFc3SsoRwJuB61PfOcDfA3vn2thX/7m0rc/+2x04U5IqbLOAc1T1AhH5MLBCVc8nEcGviMitJBMsju+oba7t+xsROQ6YStt3Qoftm8aA+q6QAfXdrsC30+enOcDXVfX7InIy9H7f9kZMexKJRCKR2szUwHokEolEAhBFJBKJRCK1iSISiUQikdpEEYlEIpFIbaKIRCKRSKQ2UUQikUgkUpsoIpFIJBKpTRSRSCQSidQmikgkEolEahNFJBKJRCK1iSISiUQikdpEEYlEIpFIbaKIRAaLiBwgIteIyKMi8jcN9nOjiBxZ8flPROREx30dmWaGbhURURHZt+3jRCJNqZ0KPjKzEZEPAfsW1KMIyXuAH6dVAmujqlvK93bU7khkxhAtkcjgSCvqASwFbuyzLZFIpJooIjMAEdlLRL4lIveLyIMi8tn0/Vki8g8islJE7hORL4vIjulny1KXyltE5Hci8oCIfCD97GiSQlWvE5F1InJt+v4eInK+iDwkIreKyF/m2vAlEfnn3OsRt5CI3CEi7xWR64DHRORHwIuAz6bH2N/4TS8Sketzry8SkStyry8VkT/N7fslZe1OWSoil6WuswvFsdxv+pu/mfbt7ZnbLX3/CRHZObftoWk/zk1f/z8icrOIrBGRH4jI0pJjHCsiN6Vtu1NE3uXStkikC6KITDhphb0LgJXAMmBPkvrVkFSFO4FksH4asBD4rLGLFwAHAC8GPigiz1DV7wP/HThbVReq6sHptt8AVgN7kFQ//O8icpRHc18PvAxYpKpHAZcC70iPcYux7S+A/URkl3RQXg7sISLbp2Vpn51+fwsV7QZ4A0kVuqcC8wDrQC0is4DvANeS9OuLgXeKyJ+o6l3Az4FXG8c4V1U3isgrSQTtVcDitK1nlRzqi8D/q6rbAwcBP7K1LRLpiigik89zSAb1d6vqY6q6XlV/mn72RuCTqnqbqq4D3g8cn3MnAfyjqj6hqteSDJYHU4CI7EVSWve96TGuAU4H/sKjrZ9W1VWq+oRtw3SbK4A/BA5P23ZZ2obnAb9RVZ+67/+uqrek+z0HcInD/D6wWFU/rKpPquptwL+xtaTs10mEMatdf3z6HsDJwEdU9WZVnSIRt0NKrJGNwIEisoOqrlHVqzx+VyTSKlFEJp+9gJXpQGWyB4mFkrGSZLLFrrn37sn9/TiJtVLEHsBDqvqosb89Pdq6ymNbgIuBI0mE5GLgJ8Afpf8u9tyX6+/Ms5TE+lmb/SOxLrL++ybwfBHZPW3jZrZaR0uBT+W+9xAgFPfXq0nqdq8UkYtF5Pmevy0SaY04O2vyWQXsLSJzCoTkLpLBLGNvYAq4F1hi2a8W7GtnEdk+JyR7A3emfz8GbJvbfjeHfdq4GPgE8Dvgo8AaEktgA/A5x3Y3YRVwu6ruV3gg1TUiciHwOuAZwDdUVXPfPU1Vv2Y7iKpeAbwyddu9g8RS2ivED4hEmhItkcnnl8DdwEdFZDsRWSAiR6SfnQX8fyKyj4gsZGu8oMhqMbkXWJbGBVDVVcDPgI+kx1gOvA34arr9NcCxIrKziOwGvDPAb/sZSbzmOcAvVfVGElF8LnCJS7sb8kvg0XRCwDYiMltEDhKR389t83USl95r2OrKAvg88H4ReSaAiOwoIn9uHkBE5onIG0VkR1XdCDxCYtFEIoMgisiEo6qbgFcA+5I8sa8meTIGOAP4CsmAezuwHvhvjrv+j/T/B0Uk89G/niR4fxfwbeBUVf1h+tlXSOIWdwAXAmfX+kE5VPUx4CrgRlV9Mn375yTuu/s82l33+JuAl5PET24HHiCJA+2Y2+x8YD/gnjSulH3328DHgG+IyCPADcAxJYd6M3BHut3JJLGsSGQQyFbrOhKJRCIRP6IlEolEIpHaRBGJRCKRCUBEzkgXDd9Q8vkbReQ6EbleRH4mIoXT9X2JIhKJRCKTwZeAoys+vx34I1V9FvBPwBdCHDRO8Y1EIpEJQFUvEZFlFZ//LPfyF9in8TvRiojsstNTdOkeRvvmzHZqwSbdOPJ6Stdv+fvJTY+PfLZ+4+jC5g0bN0zb/YYnXWar2pk/b3pXzZ87f9p7C+ZuM/J63uxtp20zRxYEaVNGvo/ymP2VYfabSVE/llHUB1DdD/nfPztJI5WQnaqpTVve2vxk7u/1W6+NTY9tbePUuifJs2Fq9LWN+XPmjbyes3D09eztkt84a8HWts6al7ues2s7d4lk13HT67cuda7Nouty5PxkNLmlCkacqnseRvttqPf8b66/6wFVXRzkwMCOyxbo1PrRmdyP37vxRpIZlBlfUNW61sTbgO/V/O4IrYjI0j2W8POzLtzyetZTpw+kGVOLRpuwTu8aeX3/+tGUSSsfuXrL37fcc93IZ7feefu0/d+++n57gx3YZ8n062PfPfeZ9t7+uy0feb10h0OnbbN4wf7T3suzUPaY9p7ZL0WYfZWR77MMs+9MivqyiqK+gOr+yPdD/jfPWZsMApvvyw0eqx/e8vejN21d2P7wiq0L7tdcOr3Nv1njvgh+v51G1+/t9MLR37Tjs7euy9z+wGTbBUu2zubNrvP8NZ0/b/nzU3Udg3//l+FyjcL067ToGi26LvNk560I8z7P43PPw/Du+/13W84xyz6wctqGDZhav5lnvHHXkfeu/OTq9ar67Kb7FpEXkYjIC5ruC1qa4nv4Mw/WuiIC5TceTL6I2G7UjCpRKRKTIiGBajGpM5BNopjAqKAUiQlsFZT89Z5d3y5iAu0IStk5gfqCAu7XahVF13EIEYFu7/1Tnnf2lSEG+IztdpunBSJiPUbqzrpAVQ8q+Xw5yRquYwqSmtaiHRE5+DC9/AduqYuaiAh0JyRFFxIUP5Xk6eKGLBMUVzEJLSQZviI7aWIC9a0TCC8ovmIC4SzpPDar2uW67fIBEuxCMg4iIiJ7k2SA/gsjPtKIKCKOtCkiIZ7ooD0xaeNJeKaKCYSxTiCspZjHR1DALiquNHHHtm2NQLWQDEFEROQskoSku5Ck9zkVmAugqp8XkdNJknlmN8dUEPfY0EUE/G6wvl1afYpInjpuAmhPSGD8xQTqxU2g2tUFzQUF/M5REzGBckHJ4yIuZcKR4eOGbdsagfIxYAgi0he9ikhZsC2KyFaqgpUZPkFLqC8moYK94NdvkygmYLdOoL6gZBSdMxcBMakSFHATFV9843hdWCNQPA585tW/jCISkjZFBPyD69CeS8v2dG0TkTIBcRGPIur2KXQT5DXx6b/QYgLTBaWpmEC5qwuqZ3WBm3UC9QQlJDZRyfARlzLRyFPHUg4tIhn58SC0iCzae76+8F2jpWUuOOX2KCJ5XAc78BMRGHZcxEVE6gpInhACnTFuYgJbBaWumEDYuAnUt07AT1CgW1FxFZSm1HXltSUisHVMiCISmC5FBMbLpdWViOTxjTtBP1ZJRlV/lll2NjGBrYKSFxPoNggP4cQEhicoeZqKi2u7+xKQPKEH+DoiIiJHA58CZgOnq+pHjc+XkpR/WExSSfNNqrq6aVt7SXtS5cOfdO5ff0uj2Sz5QdGkbD1ONqjm+32h7DEyIGVtygalpTscOjIgZQNC/sbOBvvQYpLtL9t/dsz9d1s+0qalOxy6dRDNLbZeuCgRlDlrp7b0yeb7Hh8ZpDMevWnVlkH94RUrRwb8TFDywmATlOzz/HfywrTTC/eZZgUVkV87PidXEDL7bZAISv5aun/9LSMim/VVfjDvUlDaPlZbDzHjiIjMJqnm+VKSmkFXiMj5qnpTbrOPA19W1TNF5CjgIyS1apoduxVL5PDD9LLLLpv2RO0iHkNffFTHpeUbFymzRKoEpIgyUWmSJQCaWSXmOSjrzzxDi5tAd9YJtGehQH9WShPqXmtt0rclIiLPBz6kqn+Svn4/gKp+JLfNjcDRqrpKRAR4WFV3aNrWVk2CSbQ4bl99v9PAl2flI1c3nr3iKyD575hiYlom2WCbDURFVglMf7INYZW4iEp+n1XWSd4yWbxg/y2/Z6HsseW35p/qs6f99asf3jJw5y0TaN86yfadCdeOz146ImjbH7jXFrFbsGTHkXOaf9ioY6FAf1aKC00sjS4FZCDsCeQvxNUkZaLzXAu8isTl9WfA9iLyFFV9sMmBBzXKu+SHcmHfPfeZKFO3joCUfT8vKG2JSVXf77NkceUNnv+sSlBsrq4RUrUYEZOcq6ssHabp6oKtg37eOnEVFF9XV15Q8mJitjkvKFOL5oxYYEMQlC7vxUkRj/nz5hRc/7fvIiIrcm/4JmB8F/BZETmBpCT2ncCmym840Ko7yxcXVxbUd690mQYB3F1aNneWTUTyg0tGkf8/T5Gry8fN1dTF5XsuQs+Mq+vqAjd3F9gtlCI3F1S7uqA8GA/l7i7o3uXVlXgMQThCu7N23Xc7fd2/jGYuqZoB5uLOMrZfCPxKVRungx+MJRLKChkqIVxaRRQJSNH75mBTZJ34WCZNXVyZKLgOANl2ppiYrq788UO4ukyqrBMot1BguqgUWSb5fZiuLpju7sqTd3eZv2eIFkodhiAYA+UKYD8R2YfEwjgeeEN+AxHZBXhIVTcD7yeZqdWYwYjIJHDrnbePPBnfcs91TtMc1+ldtVKflAmIbdsyQfERk1AuLpt7y6TK3dXE1QWjcQVb7MTEJihQ7vb6zZpVhVaJ6erKHycTslJ3V0X8pCtBaepWjoLhjqpOicg7gB+QTPE9Q1VvFJEPAytU9XySvFofERElcWe9PcSxB+HOaprW3MeV0qdLC9zcK+Dm0vIRkTKK3F6mq6vKLTKUmXNNZ835LGAEP3cXlLu8fDFdXVDf3QV+M7ygm9o04ygefbuz+qQVS2STbnR6uq6TEnpSqWuNNCU/8yfDxzLxDb6Xubh83VsmNndXW9aJSdHsLqi2UHwwXV35Y/hYJ+AfkAdqrUPxsUjGUUBmOq1YIocc9iz94WXnbXltDo6u8Q/X1NB9WSLQXYC9jiVS5jsvc8lk+Dy9wjAtE+jeOoF+LBRf6wTK155AmIB83WzQ4yoioS2RvZ+xs77rS3888l7oTMGhcBaRdEXkCuBOVX151bamiNSliTk9VBEB9zQoNpdWmYiUiUcRVYLSlpj0cb7acj26uLugG0FpKibQzN0VKu/aOAiJeT2FdjVNqoj8LfBsYIcuRKRpvfCuRAT6s0ZCiEieMkGxxU3GxTKB8BkH2hIUqCcqRUIC08UE/ASl6ywHQxMS2wLjmSwiTjEREVkCvAw4DfjbVlvE5MdCqqb7VsVGZj1125HBacGSHYME1zPMldIZRTO7imb+DDlmklE2u6tqqnCesqnC4BY/MWd4lcVQYLoguIhKUczE3G9R7CRrG1SfYxiNn7SV5cB31l5b+Gan6BNbAsZ0m9cCHwIUuFZV32Bu431cF0tERM4lSda1PfCuIktERE4CTgJYstceh1/960tqN6pKRIZoiUB3M7VcrJG6lkgZRRZKKDdI35ZJRh/nz2ahQLmVAm6iUmaZQL+urrbP890rN7L70rm1vpvhKiBtVDb0tUTScMMt5BIwAq/PJ2AUkf2Ac4CjVHWNiDxVVe9r2lariIjIy4FjVfWvReRISkQkTxN3lq+AwPiICDR3k0BYIXl4xcpCV0cRM1VMYPiCUhcXNxe4B+LbEJM+LBLfxKADEBGXBIz/AtyiqqeHaie4ubOOAI4TkWNJLPMdROSrqvqmsi9M6XrvlOfj7sIqSsxoLj40cXVrTS2aMzL4mG4tX8zBqGhwKhpcitxdIaaQuqyAb9vNlVHH3eU8XTi3Oh7cpgxnmG4vCCMq5vTg7FhQ7NI0FzEChee4ys1lLl60LUwNfY6bUqe8cCCqcme5JGDcH0BELiNxeX1IVb/ftFFeU3xdLZGDDt1P/+NHn97y2iYmLgIy1HrLedoM2kI4t1adwafMWvG1TtoOwo+jpQl+5zbkec3jeo77tD67EJIqK6RMQEJbIvsv31M//Z3RBeXHLPtAlSXyGpI07yemr98MPFdV35Hb5gJgI/BaYAnJqvVnqeraJm3tJO3JuFsZroRIE5+34Mwge94iaTvIblIUmIVq68Q3QFuWesPHOsn3f8gBp2gxo2uKejAKaIGzhWIW0ypKvZKdgyIR8BGWMtfmozetKrVKsjbCqGVSN2WOLfA+NKsEsuvy7L6bcSeQvzCWpO/lWQ1crqobgdtF5BZgP5L4SW1aWWxoWiJNKXUVMCxLJGNIa0faemrNCPH02vZ6ky6tE/A719As/U3TeFhG3bjYpFolrpYnJOe3ykqoQw1LZA5JYP3FJOJxBfAGVb0xt83RJMH2t6TJGK8GDmlaTySKSAuEDs7CsIUko+laBBi/QHxGG4ICfucd2p+t57qWqI2HhL7XfkH5+exbRADSuPW/sjUB42n5BIxpNcNPAEeT1BE5TVW/0bStgxeRkAIC3ZnBoX3oMB5CAuNhnUC3gw+EFRRolqRzXBak9mVtjqOI9MVgRaRKPDLGTUSg2yfTUEJSth6hag1CRpdTSIeUUw3crwEILyjQXubnMvrOudbH2iEYhjurT1oRkawD6hZhaiIgMAwRgckREtf0GzZRaTP1RqjFbTB+ggLNqmO2VVIAZoaQtCEiRQ/iB+587MwTkQxXMXERj4xxEBGo59aCdoUE2lkdbeKzWhpmhqD4Th8tKmrma6WAu6WSx0VcQpdhHpqQuIp/6Cm+UUQKTLE85tQ+X+oKCAxDRKAbl0YdqyS0kGSEFJQ+0m/AeAoKuIsK2IXFhyIBgWbJHPuIkbjcw1FEAmMTkSY0raI2pBQKQxUSCJtZ1qRJPifoTlCgvwqZeeoKCrhZKVAuKhmu4lImHBmmgMDw8m6VUeVVGIKI2BIwisjJJCVxNwHrgJPyubXqMlYiYhMQGJYrK8NHRKB7IYFurZI845IgENq/tnzzNWXUFZSMusLiS5GAQJg6JX0LSd8i4piAcQdVfST9+zjgr1X16KZt7WTFegiaCkiflK1kt+XWyqha1Q5bb8J8ri1gZHU7bBWTbPDNi0k2QOfFJBvEi8RkpxfuE0RIzH0UlX3N2mKujnfN2wWjadqrVsdD0t9FJV/z56osv1NGncHLJeOBWe43376Molxe01bMw0hmhDzmdVSEi8BUfb/ouONAWY68AfAc4FZVvQ1ARL4BvBLYIiKZgKRsR5IOvjFjYYm4CAgM05WV0Yb/G9pZS9CXVWLiWmCpbwsF3AYS1+vPN3VO1YOI63UE1TnuyqyVupQJSOhqpl2tag9dlKooE/ribfddCTyQe2tLAkaX3Fnp+28nqQk1jyQl/G+atnXwIhJKQKD/fDtdCQl0LyZtCUnGEAUF6pV/bRNfQYHq2ZM+mbjNa87H0vCpZAr9rwHKs8+SxV2JSKMEjMb2bwD+RFXf0rStgxURV/EA95t2yCIC7jESqB9EbTte4iomv1kzfZ/77VRe7z1PXUEB96A8dGulhMbmJi0TFLBPyfcRFhtt1w+Cbu77C065vW8RsdYTMbafBaxR1eo52g60IiJZQZWqC7UMH/HIGAcrJKNOmmlo1yqBfiyTIiHJ00RU+q4pDsMRFZe4m+1erbtwuA7jmOpoACLikoBxv8x9JSKvAE4N0eZWRcSk7EKtIxwZ42KFZNS1RsDPHTFJYpLhIiptCAq0KyrQjbC4FlNyffgLKSy2NWNDX1zct4iAUwLGTwEvIakpsgZ4R15k6tKpiIRm3AQkoyshgfqLz5pkiC0SlJBikmETlS7iKGAXFPAXFWhPWHwr89XxKIRiXNaFDUFE+sKlxvoCkgpY80mmBJ+rqqdWfacLEfG5wcZNRKCekEB4MYHu4ya+YpJRR1T6sFKgnqhAOGFpWuK1LWEZ14k0UUSqNkhy0G+nqutEZC7wU+AUVf1F2XfaFpFxFhBwn8LZRnA0pJhAu66uumKSUSUqfQkKuIsK+Lt6bfdGj/XBgzBU70NoETn88MP0sssuG3lvm222HaSIWBcbaqIy69KXc9N/4X1gjgxkYU8n2BYj3nLPdaVCkg1KpphkA1jVYkUYHQzLFi3C9IWL+UG3qGxrXlDyA3mRoGQiUFdM8t8zBaVokaMpduYCxyIWGK/zpW1hdJEjTF/oCNMXO8L0ksB5qoRl3EXCxPd+H+JD46TjFBNJl9RfCewLfE5V31uwzUnASQA77bbt4R867xVBG1pHPIZ6QYVcTJZRd3aNTyqMJjUrQsVOmlonGV1YKRDOUoF6s5bGhaYPh33f60OwRBxyZ80HvgwcDjwIvE5V72jaVq/AuogsAr4N/DdVvaFsu9DurLoXWN8XVhltiEhGk6maXQjKkNxdGV3FUqA7UckYkriE9CIM7d7uW0Qcc2f9NbBcVU8WkeOBP1PV1zVtq/fsLBH5IPC4qn68bJtQItLkohvaRZbHV0SgnRk1XQpKm8H4UGKS4WulgP9CR3CrtQHjKyyhXc9DvqcHICLWxYYi8oN0m5+n60ruARZrwym61piIiCwGNqrqWhHZhkTpPtbkoDbG3bStoo6AgHuyxoxsoKgSk/yAUxY7geL4CRQn6jPjJ5AISn4ADR0/yQ/6IQTFN5ZitrUolmImjITp8RQzaSQkfVuUTqRI5ItiKzBdWNpYr2USQkCGfB/nSe7p3mO1ewL5i2418NyybVR1SkQeBp7CaD4ub1yy+O4OnJmaS7OAc1T1giYHrSIKSFjyA8NMExRoLirm96tEpUhQoL6owPRAfZGoQHnAHootziKLxTeIX8W+e+4zcW7oPK3fy1OFFv4uIrIi93pLAsY+cZmddR3Qes6DSX9y6UNATIYuKJCIik1QYOtA7TrDC4ZhpYCHqBjp7WG6qIC/tQLdCEsdIYn3cCUPVLjM7gTyF+SS9L2ibVan7qwdSQLsjeh9xXoov+lMuPjamr4ZIs3F0GZ5QXeB+YxQK+jBLaYC7nEV8Eu/ntEkHUmeoa7vcMHl/g2dxffwgw/Ty39w8ch7c3ffoWnurLcDz8oF1l+lqq9t2tbeRGSSZ2pkhH5y6WoNwBBFpYtqjCFFJaSgQDNRAfeAPdQXlVBpXIZyP/tklghd2dBXRMApd9YC4CsknqWHgOOzIlZN6EVEooD40fcCsqaiMo5WCnQ3hRjaExVobq2ESteeMZR8VzZ8Mm4PQUT6onMRmXQBGbrrKhRNphBHUalnpcCwRcVXUNoSk7tXbgRg96Vza30f/Es2RBEJTJmITHr8I4SADF08yggtKmXlWENmG4YwZX/bShiZMRRRaaOkLQwjlbtJ2b1cdn8GF5FnHqw/P+vCkffmH7zbzBaRKCDljKtw2AiZiqULKwXC1ZIPndY+o46gQHeiErpGekaX932dMtbHLPtAFJGQtCUiQxSQuuIxqcJRxZBFBdoXljxNygPnCWWlQLPZXyHT3/dVKz3DxwrJrukoIoGZKSLSRfqSrql7nur+rjbT2UP31grUF5amtG2luAiKi3UyDmLiKiTjIiIisjNwNrAMuAN4raquMbZZSpIbcRZJtvbPqOrnrfseFxEZdwEZunhktJFqP1QCyb4D9VBfWKA7cakSE2hHUOpYJ0MWE19rZAxE5F+Ah1T1oyLyPmAnMxu7iMwj0YQNIrIQuAH4A9WSYFj2vXEQkaEJCPiJyLgICHRbryVEze+uXGDgLyzQr7jYxASau7zqWCc2MRmykIyxiPwaOFJV7xaR3YGfqOoBFds/BbgaeF4UkRYIVZmwiKZlR0Mk0euz8FfT2ildWivQnrDkaSIyLkICUUyK8LFGOpqdtZLRZInOubNEZK2qLkr/FmBN9trYbi/gP0lqR71bVT9n3ffQRWQmCEhb9aqhvqgMqYJkl8LSNLYC9YQF/MUlNL5uLghff35chSS0iBy630F68ae+NfLeji87wLZi/YfAbgUffQA4My8aIrJGVXeq2NcewHnAK1T13qq2RhHxxCcVgo02xTxxf7kAAA4uSURBVKOIOoIyJDHJ00RYQlgrEEZYwC4u0J3AVAkJhIuZNBWToQnJEESkCl93VvqdM4Dvquq5ldsNfbHhkEQklIB0LR5F+AjKUIXEZKjCAvXEBdwEJqOp0NgEJGMoQgLF12ZbY0aVmIyBiPz/wIO5wPrOqvoeY5sl6TZPiMhOwOXAq1X1+qp9u9QTiThiG8R8xaMqwWEZLpXtzLbYBCX7XX0vGLWJeFH7zHNi/tasH4r6bekOhxaud1i8YP/CqawLZY/CWUpz1k4VxgyKUuFnFNVZMTEFxlUEmlDWnvWrH572Wzbf9/jI756zdmqkfxbKHiP9uHjB/iP9vXSHQ0fOy/67LXd6+NlnyeJWhOT21fcXXoNj8pD1UeAcEXkbsBJ4LYCIPBs4WVVPBJ4BfEJEFBDg4zYBAQdLJA20fBnYFVCSYM6nqr4TMnfWOFkiVSLStFxtE1yFBcKn9jbp0t1QRd1zFcpiAX+rBeyWSx4fK8ZGlZhlhFpbMmTXVoZ5zYVOBR/aEmkTF0tkCvg7Vb1KRLYHrhSRi/IF4GcCbQpIW+JRtv8qUXG1UOpaJ208JZbtr+qcVVktZb97/92We1ksQHG5QvytFqi2XDLyIuMy8IfA1qa2KSt+1ZZFAluvuQEUquodl8qGdwN3p38/KiI3k9TqrSUiTcpmDpXQT7U2qirWueAqKi6Ckv/truc1u/HatjJ9xcXmDivqgzJhAX93GMDCReWxljJxga3Wi4/INMFFOHyyBfvg6tbqgjIXV1M2r98Y1JJsE68zKiLLSAqaXF7w2UnASQA77VZ+sc8kmgpImYvEd1ubwOTbEkJQhiYmJkXHCy0s4BdngXriAu4CA91YDT6FsVwwYyNF9GGNZAzJ5d4Hzmc1XQb/TeCdqvqI+Xm66OULkMREqvY1idaIDzYB8REPF3wsFxcrxSYopmVmO9f5AbyvG9LHagktLNBAXErcYhk2gTGpiseYuO63aa0Sn5hepHucRERE5pIIyNdU9Vu27ScN39TQ4D+VFMKLRxWuwmKzUkK7vIYgKHlcrZY6wgJ2qwWKz0t2/qrEBexP/2Zw30dwqrAdN3Qa+Ug1LgkY0+32Bk4H9iKZSHWsqt5RtW+riKRL5L8I3Kyqn/Rseykz0RrxndljUjXTx5WyQaeoHeZNbbNSTOG0WSnjJigZZlvqusIgvNUC1ZYLjF5HIWIUrlS1qakFMtPGkhq8D/iv3DqR9wHvLdjuy8BpqnpR6n3abNuxyxV0BPBm4HoRuSZ97+9V9btlX1gwdxuH3brRhU8zND7rQWwCEkI4XPZXdIO3LSo+bi9zoB7SNREyxgLVwgLV4gLlrkqb9WLS5NpzOUao2u1DKWY1cF4JHJn+fSbwEwwREZEDgTmqehGAqq5z2bHL7Kyfkiw88cJlBsW4WyM++bF8Z2LVuYGL1h24PmkWHc8cCIYsKjCsAaMLYYEw4pJhnl9XsXHFZQZh6Drt48qmxzYUZSDYRURW5F47J2AEdk1n2gLcQ7Luz2R/YK2IfAvYB/gh8D5V3VS141Zt2SFNxRsidRaoQfkitbrblgmN2Y6+RQXG11oBdzcYuAkL2K0WsItLhs808VDYXFZV44dNQIZ2/gPwQIMEjFtQVU1XpZvMAV5IMgP3dyQxlBNIwhmltO4QtQmJizUyTi6tIldWkRUScoVzU8r2a4qLzVqxBeuL+sFMa2HSRFjKBumhXEv5drisNQgpLFA9gLexANY1xmF78HSxPoZyjrtEVV9S9pmI3Csiu+cSMN5XsNlq4BpVvS39znnA8+hbRFwYd7dWKHyT+hXhOkXTZRaOi3ts3IQFhikudVdAm7+zzMXq4g7L0/W0WhePRRSPRpwPvIUkh9ZbgP9TsM0VwCIRWayq9wNHASsKthuhExEJ4dYaJ2vERtM6Fnl85vX7fq9IaMZZWKCeuEB3g1PTFdCuopLhKy4haCOD9KSMDS1iTcCoqptE5F3Af6Wzcq8E/s22484skRBurXGkiVugSkDqiocPZccwxcXFHWabFVbm3svEpawfM3EpG/RcxQXqC0yeJoNZG+kz8r/JZyJIX7HMcUvQuvWcDXvsUtUHgRcXvL8CODH3+iLA6wmiU3dWDLQX41qiNcNVQJrkSapKj9G2uDS1WsAuLhBOYPIMOSFf9hvqlG1ui3GsMTTkc9wHncdEqoTEZo1Mgkuryap0F/EIlWDPtp8ikalqn1lXoohMXOpaLeBuuUA4gYHxmmp665239yYkQyoh4EuX4jG17knWXDoe11QrIjJvdnXQtomQzATqzMIKJR4hjucjMCHEBZq7xcBdYMBfZPLMpOu7yW8dgnCAaxmIX3bTmAHSmiXiknmzDpNgjeRxXVRYt7xqnq6KFPlYMS7WS5WougqMzQK8f/0t1viVj8iAWy6xSWbSBWSczqNH7qyPAS9LX/6Tqp5t23er7qwqIZnpbq2uaKsmgct+q0qpltGHwEAzF1keX6HJmLRYYRSQwWHNnSUiLwMOAw4B5gM/EZHvFWVtz9N6TKQtIZlJ1LVC+i5qU0doXK0YW3xo1lO3tboFpxbNsVqCLlYM+AkN2FfvuzJU8al7/0YBaQ1r7izgQOASVZ0CpkTkOuBo4JyqHXcSWK8rJFVEa6QaXwEpyNPTmB2fvdS6ja2dQxcZ8Bca8Jv67VrOOAQhRaluCeWm3L1yIwC7L50bfN9FApKcA6vXJwRt5866FjhVRD4BbAu8CIcKtp3NzqoTI5kJ1sg6vSt4pl5XAWlDOOruv0xwfK2ZLkUGwgtNRl3ByVMnJtkk1lOGj5gM4eFwCFN4N0w9yW/WTLv2W82dpaoXisjvAz8D7gd+DlQmX4SOp/iWCcnQrZGyVcQuUyVXPnL1yCBw//pbOi0+VUTb4lEH1zYViY2P0LhMRFiwZEfnFf0hhQbcYjRF2AqKueIiPK4LO01cxaTpfd2GBQJVVsgwCJA7C1U9DTgt/c7XAWtWzlZEZI4sKP3MV0jGzRq55Z7rgl9Ys566bScr1F1pY/76Ti908zPXFZuhCw34iU1GSNEBu/D4uNZss9NcMgT0YZU0eWgcMNbcWSIyG1ikqg+KyHKSlesX2nbcmiWyeMH+3qml6wjJEMzf0EwtmtNaxt66tL3wyXf/NtFxEZuhCA34iQ24zTwzsWUCKMKW3j+jTulklyqXfYqJC208NLaENXcWMBe4NEmbxSPAm9IgeyWiWpRWvhmHHPYs/eFl5wHlTzyhi8+0fZGV+UldTVzz5jNvYnMwKBpQigakssHMNiD6uLTGZeVsGa5WThkuEwTyVK2hKaIqxYyJbw103/K3roWobA+IVRbLONQI8Zmdtf9uyzlm2QeurIpXeB9/1i566oKXj7z31ifODHqMULQeEymzSCbFrTXmJu6MwEcEiwTH16pxndjga9WAf6p/36JkLmtqwD+nWdn6mbKV/31bJj7ZlIc6zborWrFEROR+EpNpaOwCPNB3IyqI7WvG0NsHw29jbF89lqpqsGldIvJ9kt+a5wFVPTrUMULRiogMFRFZMURzMCO2rxlDbx8Mv42xfRFfZvXdgEgkEomML1FEIpFIJFKbmSYirikC+iK2rxlDbx8Mv42xfREvZlRMJBKJRCJhmWmWSCQSiUQCEkUkEolEIrWZOBERkTNE5D4RuaHk8yNF5GERuSb998GO27eXiPxYRG4SkRtF5JSCbUREPi0it4rIdSJy2MDa11sfisgCEfmliFybtu8fC7aZLyJnp/13uYgsG1j7ThCR+3P9d2JX7cu1YbaIXC0iFxR81lv/Ge2oamPvfRhJ6DSLb0d8Cfgs8OWKbS5V1ZdXfN4mU8DfqepVIrI9cKWIXKSq+bz9xwD7pf+eC/yv9P+htA/668MNwFGquk5E5gI/Tauv/SK3zduANaq6r4gcD3wMeN2A2gdwtqq+o6M2FXEKcDOwQ8FnffZfnqo2Qv99GGECLRFVvQR4qO92lKGqd6vqVenfj5LcJHsam70S+LIm/AJYlKZvHkr7eiPtk3Xpy7npP3N2yCtJqrcBnAu8WNKscgNpX6+IyBKSOtqnl2zSW/9lOLQxMhAmTkQceX7qbvieiDyzr0akboJDgcuNj/YE8gmYVtPDQF7RPuixD1M3xzUkNREuUtXS/kuzkD4MPGVA7QN4deqqPFdE/DI2NudfgfcAm0s+77X/UmxthH77MJIyE0XkKpI8NwcDnwHO66MRIrIQ+CbwTlV9pI82VGFpX699qKqbVPUQYAnwHBE5qMvj23Bo33eAZaq6HLiIrU/9rSMiLwfuU9UruzqmL45t7K0PI6PMOBFR1Ucyd4OqfheYKyJmorNWSX3l3wS+pqrfKtjkTiD/ZLUkfa8TbO0bQh+mx14L/Bgwk9Jt6T8RmQPsCDzYbevK26eqD6rqhvTl6cDhHTbrCOA4EbkD+AZwlIh81dim7/6ztrHnPozkmHEiIiK7Zf5dEXkOSR90doOkx/4icLOqfrJks/OBv0hnaT0PeFhV7x5K+/rsQxFZLCKL0r+3AV4K/MrYLKviBvAa4Efa0apal/YZ8a3jSOJOnaCq71fVJaq6DDiepG/eZGzWW/+5trHPPoyMMnGzs0TkLOBIYBcRWQ2cShLcRFU/T3JT/JWITAFPAMd3eYOQPGW9Gbg+9ZsD/D2wd66N3wWOBW4FHgfeOrD29dmHuwNnSlLKcxZwjqpeICIfBlao6vkkIvgVEbmVZJLF8R21zbV9fyMix5HMhHsIOKHD9hUyoP4rZeh9OFOJaU8ikUgkUpsZ586KRCKRSDiiiEQikUikNlFEIpFIJFKbKCKRSCQSqU0UkUgkEonUJopIJBKJRGoTRSQSiUQitfm/Qj3CFQE6FvQAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import BoundaryNorm\n", + "from matplotlib.ticker import MaxNLocator\n", + "import numpy as np\n", + "\n", + "\n", + "# make these smaller to increase the resolution\n", + "dx, dy = 0.05, 0.05\n", + "\n", + "# generate 2 2d grids for the x & y bounds\n", + "y, x = np.mgrid[slice(1, 5 + dy, dy),\n", + " slice(1, 5 + dx, dx)]\n", + "\n", + "z = np.sin(x)**10 + np.cos(10 + y*x) * np.cos(x)\n", + "\n", + "# x and y are bounds, so z should be the value *inside* those bounds.\n", + "# Therefore, remove the last value from the z array.\n", + "z = z[:-1, :-1]\n", + "levels = MaxNLocator(nbins=15).tick_values(z.min(), z.max())\n", + "\n", + "\n", + "# pick the desired colormap, sensible levels, and define a normalization\n", + "# instance which takes data values and translates those into levels.\n", + "cmap = plt.get_cmap('PiYG')\n", + "norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)\n", + "\n", + "fig, (ax0, ax1) = plt.subplots(nrows=2)\n", + "\n", + "im = ax0.pcolormesh(x, y, z, cmap=cmap, norm=norm)\n", + "fig.colorbar(im, ax=ax0)\n", + "ax0.set_title('pcolormesh with levels')\n", + "\n", + "\n", + "# contours are *point* based plots, so convert our bound into point\n", + "# centers\n", + "cf = ax1.contourf(x[:-1, :-1] + dx/2.,\n", + " y[:-1, :-1] + dy/2., z, levels=levels,\n", + " cmap=cmap)\n", + "fig.colorbar(cf, ax=ax1)\n", + "ax1.set_title('contourf with levels')\n", + "\n", + "# adjust spacing between subplots so `ax1` title and `ax0` tick labels\n", + "# don't overlap\n", + "fig.tight_layout()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1. , 1. , 1. , ..., 1. , 1. , 1. ],\n", + " [1.05, 1.05, 1.05, ..., 1.05, 1.05, 1.05],\n", + " [1.1 , 1.1 , 1.1 , ..., 1.1 , 1.1 , 1.1 ],\n", + " ...,\n", + " [4.9 , 4.9 , 4.9 , ..., 4.9 , 4.9 , 4.9 ],\n", + " [4.95, 4.95, 4.95, ..., 4.95, 4.95, 4.95],\n", + " [5. , 5. , 5. , ..., 5. , 5. , 5. ]])" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}