{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Deep Learning: Downscale climate models\n", "\n", "> * 🔬 Data Science\n", "* 🤖 Machine Learning\n", "* 🥠 Predictive Modeling\n", "* 📈 Graphing and Visualization\n", "* 🖥️ Requires ArcPy\n", "* 👟 Ready To Run!\n", "\n", "In this notebook, an essential component for making climate projections is addressed with spatial and non-spatial machine learning methods, namely climate downscaling. Downscaling is calibrating a global climate model with spatially-limited local temperature measurements.\n", "\n", "For climate models, calibration can be done in various ways. One way of performing calibration is the so-called statistical downscaling that forms a statistical relationship between global climate models and temperature observations. For the inherently spatial problem of downscaling, we evaluate regression methods based on their method of handling spatial information.\n", "\n", "1. Non-spatial (generic) machine learning\n", "2. Spatial machine learning\n", "3. Non-spatial machine learning with geoenriched predictors\n", " \n", "## What is Statistical Climate Downscaling\n", "### Overview\n", "Scientists use Global Climate Models (GCMs) to model mesoscale to large scale (100s of kilometers) atmospheric dynamics of our planet. GCMs solve analytical equations, such as Navier-Stokes equations, to model global air flow and energy transfer in the atmosphere. In addition to global dynamics, weather is impacted by local effects such as distance to a water body, proximity to mountains, etc. These local effects are not modeled in GCMs. Thus, GCMs need to be calibrated with measurements collected at a finer scale. Fine-scale measurements of weather such as local temperature comes from weather stations that collect data continuously. In this analysis, we will be using weather measurements from a subset of weather stations in U.S. and CGCM1 climate model. Statistical downscaling is summarized below:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Predicting Climate Over Contiguous U.S.\n", "Predicting future temperature using a GCM requires forming a statistically-valid relationship between coarse-scale climate model and fine-scale weather measurements using observed data. We will use 19 variables from CGCM1:\n", "\n", "- Mean sea level pressure\n", "- Airflow strength (simulated at 3 pressure levels)\n", "- Surface zonal velocity (simulated at 3 pressure levels)\n", "- Meridional velocity (simulated at 3 pressure levels)\n", "- Surface vorticity (simulated at 3 pressure levels)\n", "- Geopotential height (simulated 2 pressure levels)\n", "- Zonal velocity (simulated at 2 pressure levels)\n", "- Near surface relative humidity\n", "- Near surface specific humidity\n", "\n", "We will establish a relationship between simulated climate variables and measured average temperature to calibrate GCM. This process is known as downscaling and it can be done a number of ways (Wilby and Wigley, 1997). In this blog we will explore statistical downscaling. For the sake of demonstration, we will downscale GCM variables for a time snapshot: March 12, 2012 in the contiguous United States. \n", "Downscaling methods we will explore are summarized below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Approach\n", "#### (1) Introduce Climate Downscaling\n", "#### (2) Explain Regression Problem\n", "#### (3) Methods we used: SVM/Ridge Regression method/GWR/RF" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's get started. First, we will import libraries that contain spatial and non-spatial regression models used throughout the study. In this notebook, we will use ArcPy and ArcGIS API for Python to perform spatial machine learning functions and data management, and scikit-learn will be used for non-spatial machine learning methods." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", "import shutil\n", "import zipfile\n", "\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt\n", "import matplotlib.cm as cm\n", "import pandas as pd\n", "import numpy as np\n", "from sklearn.linear_model import Ridge, RidgeCV\n", "from sklearn.svm import SVR\n", "from IPython.display import display\n", "from IPython.display import IFrame\n", "\n", "from arcgis.gis import GIS\n", "import arcpy" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "IFrame(src='https://arcgis.com/apps/Minimalist/index.html'\\\n", " '?appid=295c7fed9bbb413f8d8394b0effe1f71',\n", " width=\"100%\", height=\"800px\") " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even though it is a must, we recommend making sure ArcPy is properly imported by running the GetInstallInfo function and querying the version of the ArcPy installed on ArcGIS Notebook Server." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'10.7.0'" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arcpy.GetInstallInfo()['Version']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Connect to GIS. GIS enables interactive widgets to be embedded in analytical work and internally uses the ArcGIS API for JavaScript." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "agol_gis = GIS(set_active=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Automatically authenticate to your WebGIS using your credentials" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "gis = GIS(\"home\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 0. Environment setting " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 0.1 Unzip uploaded GDB " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from shutil import unpack_archive" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's copy the sample geodatabase used from this sample from `/arcgis/samplesdata/` to `/arcgis/home`" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "GDB succesfully copied to /arcgis/home/climate_input.gdb/\n" ] } ], "source": [ "arcgis_dir = os.path.abspath(os.path.join(os.sep, 'arcgis', ''))\n", "home_dir = os.path.join(arcgis_dir, 'home', '')\n", "\n", "def copy_sample_gdb_to_home(gdb_zip_name):\n", " \"\"\"Given the full filename (with extensions) of a gdb zip file in \n", " /arcgis/samplesdata/, will copy and unzip that gdb to /arcgis/home/\n", " Will return the full path to the unzipped gdb in home\"\"\"\n", "\n", " # Get the full paths of all the source and destination files to copy\n", " gdb_dir_name = gdb_zip_name.split(\".zip\")[0]\n", " gdb_zip_path_src = os.path.join(arcgis_dir, 'samplesdata', \n", " gdb_zip_name)\n", " gdb_dir_path_src = os.path.join(arcgis_dir, 'samplesdata', \n", " gdb_dir_name)\n", " gdb_zip_path_dst = os.path.join(home_dir, gdb_zip_name)\n", " gdb_dir_path_dst = os.path.join(home_dir, gdb_dir_name)\n", "\n", " # If the gdb has been copied/unzipped to home dir before, delete it\n", " if os.path.exists(gdb_zip_path_dst):\n", " os.remove(gdb_zip_path_dst)\n", " if os.path.exists(gdb_dir_path_dst):\n", " shutil.rmtree(gdb_dir_path_dst)\n", "\n", " # Copy the zip file to home, unzip it\n", " shutil.copy(gdb_zip_path_src, gdb_zip_path_dst)\n", " zip_ref = zipfile.ZipFile(gdb_zip_path_dst, 'r')\n", " zip_ref.extractall(home_dir)\n", " zip_ref.close()\n", " \n", " # Return the output full path to /arcgis/home/unzipped_gdb/\n", " return os.path.join(gdb_dir_path_dst, '') # Adds trailing slash\n", "\n", "gdb_path = copy_sample_gdb_to_home('climate_input.gdb.zip')\n", "print(f\"GDB succesfully copied to {gdb_path}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 0.2 Set up working environment" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "arcpy.env.workspace = gdb_path\n", "outputSpace = home_dir\n", "climate_output_gdb_path = os.path.join(home_dir, \"climate_output.gdb\")\n", "if os.path.exists(climate_output_gdb_path):\n", " shutil.rmtree(climate_output_gdb_path)\n", "arcpy.CreateFileGDB_management(outputSpace, \"climate_output.gdb\")\n", "outGDB = os.path.join(outputSpace, \"climate_output.gdb\")\n", "\n", "arcpy.env.overwriteOutput = True ## overwrite the output results" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 0.3 List the feature class in my working environment" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "datasets = arcpy.ListDatasets(feature_type='feature')\n", "datasets = [''] + datasets if datasets is not None else []" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Data Wrangling \n", "### 1.1 Filling Missing Station Values\n", "#### There are a lot of them and we don't have much data to start with, so let's do FIMV.\n", "#### There are some average surface temperature values are missing, and we need to fill the missing values." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "station_fc = \"station\"" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "station_sdf = pd.DataFrame.spatial.from_featureclass(station_fc)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "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", "
T_MONTHLY_MAXT_MONTHLY_MINT_MONTHLY_MEANT_MONTHLY_AVG
0-9999.011.417.8-9999.0
1-9999.0-9999.0-9999.016.7
222.69.316.0-9999.0
3-9999.0-9999.018.818.8
421.69.515.615.8
\n", "
" ], "text/plain": [ " T_MONTHLY_MAX T_MONTHLY_MIN T_MONTHLY_MEAN T_MONTHLY_AVG\n", "0 -9999.0 11.4 17.8 -9999.0\n", "1 -9999.0 -9999.0 -9999.0 16.7\n", "2 22.6 9.3 16.0 -9999.0\n", "3 -9999.0 -9999.0 18.8 18.8\n", "4 21.6 9.5 15.6 15.8" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# show the temperature info from my station data\n", "station_sdf[[\"T_MONTHLY_MAX\", \"T_MONTHLY_MIN\", \"T_MONTHLY_MEAN\",\n", " \"T_MONTHLY_AVG\"]].head()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# As you can see, there are many -9999.0 in SUR_TEMP_MONTHLY_AVG. \n", "# -9999 represents missing values, and we need to use fill mising value \n", "# tool to fill those fields\n", "# You can use the number of spatial neighbors for interpolating the \n", "# missing values; here we'll use five numbers of neighbors\n", "\n", "station_fmv = os.path.join(outGDB, \"station_filled_missing_value\")\n", "\n", "arcpy.stpm.FillMissingValues(station_fc, station_fmv,\n", " \"T_MONTHLY_MAX;T_MONTHLY_MIN;T_MONTHLY_MEAN;T_MONTHLY_AVG\",\n", " \"AVERAGE\", \"K_NEAREST_NEIGHBORS\", None, None, \"TIME\", 6, \"WBANNO\",\n", " None, None, None, None, -9999, None)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "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", "
T_MONTHLY_MAXT_MONTHLY_MINT_MONTHLY_MEANT_MONTHLY_AVG
024.5411.4017.8018.04
123.069.8216.4816.70
222.609.3016.0016.48
324.8012.0418.8018.80
421.609.5015.6015.80
\n", "
" ], "text/plain": [ " T_MONTHLY_MAX T_MONTHLY_MIN T_MONTHLY_MEAN T_MONTHLY_AVG\n", "0 24.54 11.40 17.80 18.04\n", "1 23.06 9.82 16.48 16.70\n", "2 22.60 9.30 16.00 16.48\n", "3 24.80 12.04 18.80 18.80\n", "4 21.60 9.50 15.60 15.80" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "station_fmv_sdf = pd.DataFrame.spatial.from_featureclass(station_fmv)\n", "\n", "station_fmv_sdf[[\"T_MONTHLY_MAX\", \"T_MONTHLY_MIN\",\n", " \"T_MONTHLY_MEAN\", \"T_MONTHLY_AVG\"]].head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 Plot Local Station Data" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# create a function which is able to plot the spatially enabled data frame\n", "def visualizeFeatureClassBySEDF(myMap, fc, renderer_type,\n", " method, class_count, col,\n", " colors, alpha, line_width):\n", " sdf = pd.DataFrame.spatial.from_featureclass(fc)\n", " # change the spatial reference of the SEDF same as the map\n", " sdf.spatial.project(spatial_reference = {'wkid': 3857}) \n", " sdf.spatial.plot(map_widget= myMap,renderer_type=renderer_type,\n", " method=method, class_count=class_count,\n", " col=col,colors=colors,alpha=alpha,\n", " line_width=line_width)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": true }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "879c6a16b3d1425aab14335c83616150", "version_major": 2, "version_minor": 0 }, "text/plain": [ "MapView(layout=Layout(height='400px', width='100%'), zoom=4.0)" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sdf_map = agol_gis.map('United States',zoomlevel=4)\n", "sdf_map" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "visualizeFeatureClassBySEDF(sdf_map, station_fmv, renderer_type='c',\n", " method='esriClassifyNaturalBreaks',\n", " class_count=7,col=\"T_MONTHLY_AVG\",\n", " colors='Reds' ,alpha=0.5, line_width=0.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.3 Resolve Scale and Discretization Differences\n", "Global climate modeling is simulated at locations around the world. One of the global climate variables is visualized." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# First, let's plot the GCM data\n", "\n", "def visualizeFeatureClassBySEDF_unique(myMap, fc, symbol_type,\n", " cmap,line_width,marker_size, \n", " alpha):\n", " sdf = pd.DataFrame.spatial.from_featureclass(fc)\n", " # change the spatial reference of the SEDF same as the map\n", " sdf.spatial.project(spatial_reference = {'wkid': 3857})\n", " sdf.spatial.plot(map_widget= myMap,symbol_type=symbol_type,\n", " cmap=cmap,line_width=line_width, \n", " marker_size=marker_size)" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [], "source": [ "# We add the GCM data to the map above; as you can see, the locations of GCM data points do not coincide with the locations of stations\n", "visualizeFeatureClassBySEDF_unique(sdf_map, \"airFlow500_201203\", \n", " symbol_type='simple',\n", " cmap=[90,180,172,225],\n", " line_width= 0.5,\n", " marker_size=7,\n", " alpha=\"0.5\")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Calling arcpy.ga.EmpiricalBayesianKriging(\"specHum2m_201203\", \"shum\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"airT2m_201203\", \"air\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"geoHgt850_201203\", \"hgt\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"geoHgt500_201203\", \"hgt\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"relHumSurf_201203\", \"rhum\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"specHum850_201203\", \"shum\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"specHum500_201203\", \"shum\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"vort850_201203\", \"omega\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"vort500_201203\", \"omega\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"merid850_201203\", \"vwnd\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"merid500_201203\", \"vwnd\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"zonal850_201203\", \"uwnd\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"zonal500_201203\", \"uwnd\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"airFlow850_201203\", \"wspd\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"airFlow500_201203\", \"wspd\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"airFlowSurf_201203\", \"wspd\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"vortSurf_201203\", \"wspd\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"meridSurf_201203\", \"vwnd\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"seaLevP_201203\", \"slp\", ...\n", "Calling arcpy.ga.EmpiricalBayesianKriging(\"zonalSurf_201203\", \"uwnd\", ...\n" ] } ], "source": [ "# Since the locations of GCM data do not coincide with the\n", "# locations of the real station, first, we need to interpolate GCM data\n", "# and get the raster surfaces of GCM \n", "\n", "arcpy.env.extent = arcpy.Extent(-124.76308372794, 24.5230996116049, -66.9498936574818, 49.3843680121063) # setting the extent as US polygon - this is important, as otherwise some parts of the US won't get the interpolated values\n", "# run through a loop to create all GCM raster surfaces using EBK\n", "gcm_lst = [(\"specHum2m_201203\",\"shum\"),(\"airT2m_201203\",\"air\"),\n", " (\"geoHgt850_201203\",\"hgt\"),(\"geoHgt500_201203\",\"hgt\"),\n", " (\"relHumSurf_201203\",\"rhum\"),(\"specHum850_201203\",\"shum\"),\n", " (\"specHum500_201203\",\"shum\"),(\"vort850_201203\",\"omega\"),\n", " (\"vort500_201203\",\"omega\"),(\"merid850_201203\",\"vwnd\"),\n", " (\"merid500_201203\",\"vwnd\"),(\"zonal850_201203\",\"uwnd\"),\n", " (\"zonal500_201203\",\"uwnd\"),(\"airFlow850_201203\",\"wspd\"),\n", " (\"airFlow500_201203\",\"wspd\"),(\"airFlowSurf_201203\",\"wspd\"),\n", " (\"vortSurf_201203\",\"wspd\"),(\"meridSurf_201203\",\"vwnd\"),\n", " (\"seaLevP_201203\",\"slp\"),(\"zonalSurf_201203\",\"uwnd\")]\n", "\n", "for gcm,value in gcm_lst:\n", " arcpy.env.workspace = home_dir + \"/climate_input.gdb\"\n", " file_name = str(gcm)+\"_ebk_surface\"\n", " output = os.path.join(outGDB, file_name)\n", " print(f'Calling arcpy.ga.EmpiricalBayesianKriging('\\\n", " f'\"{gcm}\", \"{value}\", ...')\n", " arcpy.ga.EmpiricalBayesianKriging(gcm, value, None, output,\n", " 9.94450736020055E-02, \"NONE\", 100, 1, 100,\n", " \"NBRTYPE=StandardCircular RADIUS=15 ANGLE=0 \"\\\n", " \"NBR_MAX=15 NBR_MIN=10 SECTOR_TYPE=ONE_SECTOR\",\n", " \"PREDICTION\", 0.5, \"EXCEED\", None, \"POWER\")" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "# We can list all raster files in the output gdb to make sure our \n", "# EBK outputs are in there\n", "\n", "# Reset the current workspace\n", "arcpy.env.workspace = home_dir + \"climate_output.gdb\" " ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# extract the values from the raster surface to the points of stations. \n", "# The \"extract multi value to points\" method can help us to do that\n", "arcpy.env.workspace = home_dir + \"/climate_input.gdb\"\n", "\n", "# since we don't want to change the original input station data, \n", "# we make a copy to climate_output.gdb\n", "arcpy.CopyFeatures_management(station_fmv, \n", " home_dir + \"/climate_output.gdb/station_GCM\") \n", "\n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\"\n", "arcpy.sa.ExtractMultiValuesToPoints(\"station_GCM\", \n", " [[\"specHum2m_201203_ebk_surface\", \"specHum2m\"], \n", " [\"airT2m_201203_ebk_surface\", \"airT2m\"], \n", " [\"geoHgt850_201203_ebk_surface\", \"geoHgt850\"],\n", " [\"geoHgt500_201203_ebk_surface\",\"geoHgt500\"],\n", " [\"relHumSurf_201203_ebk_surface\",\"relHumSurf\"],\n", " [\"specHum500_201203_ebk_surface\",\"specHum500\"],\n", " [\"specHum850_201203_ebk_surface\",\"specHum850\"],\n", " [\"vort850_201203_ebk_surface\",\"vort850\"],\n", " [\"vort500_201203_ebk_surface\",\"vort500\"],\n", " [\"merid850_201203_ebk_surface\",\"merid850\"],\n", " [\"merid500_201203_ebk_surface\",\"merid500\"], \n", " [\"zonal850_201203_ebk_surface\",\"zonal850\"],\n", " [\"zonal500_201203_ebk_surface\",\"zonal500\"],\n", " [\"airFlow850_201203_ebk_surface\",\"airFlow850\"],\n", " [\"airFlow500_201203_ebk_surface\",\"airFlow500\"],\n", " [\"airFlowSurf_201203_ebk_surface\",\"airFlowSurf\"],\n", " [\"vortSurf_201203_ebk_surface\",\"vortSurf\"],\n", " [\"meridSurf_201203_ebk_surface\",\"meridSurf\"],\n", " [\"seaLevP_201203_ebk_surface\",\"seaLevP\"],\n", " [\"zonalSurf_201203_ebk_surface\",\"zonalSurf\"]],\n", " \"NONE\")" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "# now the station data gets the values from 20 GCM raster surfaces, \n", "# we can check them by listing the files in our data\n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\"\n", "fields = arcpy.ListFields(\"station_GCM\")\n", "station_predictors = [\"airT2m\",\"geoHgt850\",\"geoHgt850\",\n", " \"geoHgt500\",\"relHumSurf\",\"specHum500\",\n", " \"specHum850\",\"vort850\",\"vort500\",\"merid850\",\n", " \"vort500\",\"merid850\", \"merid500\",\"zonal850\",\n", " \"zonal500\",\"airFlow850\",\"airFlow500\",\"airFlowSurf\",\n", " \"vortSurf\",\"meridSurf\",\"seaLevP\",\"zonalSurf\"]\n", "station_predictand = [\"T_MONTHLY_AVG\"]" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "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", "
airT2mgeoHgt850geoHgt850geoHgt500relHumSurfspecHum500specHum850vort850vort500merid850...zonal850zonal500airFlow850airFlow500airFlowSurfvortSurfmeridSurfseaLevPzonalSurfT_MONTHLY_AVG
0289.849213164.443024164.443024164.44302484.12443510.73386810.733868-0.000187-0.0001871.560632...-0.228055-0.2280554.3007534.3007533.4807433.4807431.6411731019.5000610.16813418.04
1288.776917159.583237159.583237159.58323781.95317110.15522110.1552210.0097560.0097561.874495...0.0811500.0811504.5692214.5692214.8042244.8042242.3062591018.8183590.03188616.70
2288.779755161.611938161.611938161.61193882.47177910.26257010.2625700.0070050.0070051.700716...0.0112860.0112864.3801004.3801004.4022244.4022242.0130571019.1153560.13222516.48
3292.476349165.336807165.336807165.33680783.30391711.52295411.522954-0.002438-0.0024382.187343...-0.947112-0.9471124.9866184.9866183.7407923.7407921.5715671019.343567-0.97135018.80
4288.211273163.101761163.101761163.10176182.62401610.28277410.2827740.0056270.0056271.584904...0.0275540.0275544.1575634.1575634.3063434.3063431.7470681019.3632810.21642115.80
\n", "

5 rows × 23 columns

\n", "
" ], "text/plain": [ " airT2m geoHgt850 geoHgt850 geoHgt500 relHumSurf specHum500 \\\n", "0 289.849213 164.443024 164.443024 164.443024 84.124435 10.733868 \n", "1 288.776917 159.583237 159.583237 159.583237 81.953171 10.155221 \n", "2 288.779755 161.611938 161.611938 161.611938 82.471779 10.262570 \n", "3 292.476349 165.336807 165.336807 165.336807 83.303917 11.522954 \n", "4 288.211273 163.101761 163.101761 163.101761 82.624016 10.282774 \n", "\n", " specHum850 vort850 vort500 merid850 ... zonal850 \\\n", "0 10.733868 -0.000187 -0.000187 1.560632 ... -0.228055 \n", "1 10.155221 0.009756 0.009756 1.874495 ... 0.081150 \n", "2 10.262570 0.007005 0.007005 1.700716 ... 0.011286 \n", "3 11.522954 -0.002438 -0.002438 2.187343 ... -0.947112 \n", "4 10.282774 0.005627 0.005627 1.584904 ... 0.027554 \n", "\n", " zonal500 airFlow850 airFlow500 airFlowSurf vortSurf meridSurf \\\n", "0 -0.228055 4.300753 4.300753 3.480743 3.480743 1.641173 \n", "1 0.081150 4.569221 4.569221 4.804224 4.804224 2.306259 \n", "2 0.011286 4.380100 4.380100 4.402224 4.402224 2.013057 \n", "3 -0.947112 4.986618 4.986618 3.740792 3.740792 1.571567 \n", "4 0.027554 4.157563 4.157563 4.306343 4.306343 1.747068 \n", "\n", " seaLevP zonalSurf T_MONTHLY_AVG \n", "0 1019.500061 0.168134 18.04 \n", "1 1018.818359 0.031886 16.70 \n", "2 1019.115356 0.132225 16.48 \n", "3 1019.343567 -0.971350 18.80 \n", "4 1019.363281 0.216421 15.80 \n", "\n", "[5 rows x 23 columns]" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# here are the variables we will use later\n", "station_GCM = os.path.join(outGDB, \"station_GCM\")\n", "station_GCM_sdf = pd.DataFrame.spatial.from_featureclass(station_GCM)\n", "station_GCM_sdf[station_predictors + station_predictand].head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.4 Create Prediction Location\n", "\n", "Now, we have almost prepared all input data for our analysis, but we also need to prepare prediction locations. \"US_Prediction_Grid\" is the predicted grid that we created in advance, and we need to use the \"Extract multi value to points\" method again to extract values from raster surface data." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "arcpy.env.workspace = home_dir + \"/climate_input.gdb\"\n", "station = \"US_Prediction_Grid\"\n", "\n", "# since we don't want to change the original input US_Prediction_Grid \n", "# data, we make a copy to climate_output.gdb\n", "arcpy.CopyFeatures_management(station,\n", " home_dir + \"/climate_output.gdb/US_Prediction_Grid_GCM\") \n", "\n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\"\n", "arcpy.sa.ExtractMultiValuesToPoints(\"US_Prediction_Grid_GCM\",\n", " [[\"specHum2m_201203_ebk_surface\", \"specHum2m\"], \n", " [\"airT2m_201203_ebk_surface\", \"airT2m\"], \n", " [\"geoHgt850_201203_ebk_surface\", \"geoHgt850\"],\n", " [\"geoHgt500_201203_ebk_surface\",\"geoHgt500\"],\n", " [\"relHumSurf_201203_ebk_surface\",\"relHumSurf\"],\n", " [\"specHum500_201203_ebk_surface\",\"specHum500\"],\n", " [\"specHum850_201203_ebk_surface\",\"specHum850\"],\n", " [\"vort850_201203_ebk_surface\",\"vort850\"],\n", " [\"vort500_201203_ebk_surface\",\"vort500\"],\n", " [\"merid850_201203_ebk_surface\",\"merid850\"],\n", " [\"merid500_201203_ebk_surface\",\"merid500\"], \n", " [\"zonal850_201203_ebk_surface\",\"zonal850\"],\n", " [\"zonal500_201203_ebk_surface\",\"zonal500\"],\n", " [\"airFlow850_201203_ebk_surface\",\"airFlow850\"],\n", " [\"airFlow500_201203_ebk_surface\",\"airFlow500\"],\n", " [\"airFlowSurf_201203_ebk_surface\",\"airFlowSurf\"],\n", " [\"vortSurf_201203_ebk_surface\",\"vortSurf\"],\n", " [\"meridSurf_201203_ebk_surface\",\"meridSurf\"],\n", " [\"seaLevP_201203_ebk_surface\",\"seaLevP\"],\n", " [\"zonalSurf_201203_ebk_surface\",\"zonalSurf\"]], \n", " \"NONE\")" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "# now the US_Prediction_Grid data gets the values from 20 GCM raster \n", "# surfaces. We can check them by listing the files in our data\n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\" \n", "fields = arcpy.ListFields(\"US_Prediction_Grid_GCM\")" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "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", "
OBJECTIDFID_USdownFID_US_polairFlow850airFlowSurairT2m_201geoHgt500_geoHgt850_merid500_2merid850_2...zonal850zonal500airFlow850_1airFlow500_1airFlowSurfvortSurfmeridSurfseaLevPzonalSurfSHAPE
01656677483.785594.60107286.886993117.762001117.7620010.8948560.894856...1.9824751.9824754.0889004.0889004.6385794.6385790.9745921013.3960571.815351{\"x\": -105.51129143475185, \"y\": 31.06161320108...
12656687483.810094.62228287.040985117.556000117.5560000.9277910.927791...1.9721881.9721884.1078874.1078874.6657584.6657581.0193391013.3709111.796262{\"x\": -105.39566505429877, \"y\": 31.06161320108...
23656697483.861464.66634287.350006117.151001117.1510010.9971580.997158...1.9606181.9606184.1273014.1273014.6929184.6929181.0654411013.3469241.774962{\"x\": -105.2800386738457, \"y\": 31.061613201088...
34656707483.888354.68925287.489014116.954002116.9540021.0335901.033590...1.9315061.9315064.1683434.1683434.7564234.7564231.1779401013.2924191.715283{\"x\": -105.16441229429194, \"y\": 31.06161320108...
45656717483.916074.71277287.639008116.760002116.7600021.0711801.071180...1.9161671.9161674.1893044.1893044.7838134.7838131.2280211013.2719121.687406{\"x\": -105.0487859138388, \"y\": 31.061613201088...
\n", "

5 rows × 43 columns

\n", "
" ], "text/plain": [ " OBJECTID FID_USdown FID_US_pol airFlow850 airFlowSur airT2m_201 \\\n", "0 1 65667 748 3.78559 4.60107 286.886993 \n", "1 2 65668 748 3.81009 4.62228 287.040985 \n", "2 3 65669 748 3.86146 4.66634 287.350006 \n", "3 4 65670 748 3.88835 4.68925 287.489014 \n", "4 5 65671 748 3.91607 4.71277 287.639008 \n", "\n", " geoHgt500_ geoHgt850_ merid500_2 merid850_2 \\\n", "0 117.762001 117.762001 0.894856 0.894856 \n", "1 117.556000 117.556000 0.927791 0.927791 \n", "2 117.151001 117.151001 0.997158 0.997158 \n", "3 116.954002 116.954002 1.033590 1.033590 \n", "4 116.760002 116.760002 1.071180 1.071180 \n", "\n", " ... zonal850 zonal500 \\\n", "0 ... 1.982475 1.982475 \n", "1 ... 1.972188 1.972188 \n", "2 ... 1.960618 1.960618 \n", "3 ... 1.931506 1.931506 \n", "4 ... 1.916167 1.916167 \n", "\n", " airFlow850_1 airFlow500_1 airFlowSurf vortSurf meridSurf seaLevP \\\n", "0 4.088900 4.088900 4.638579 4.638579 0.974592 1013.396057 \n", "1 4.107887 4.107887 4.665758 4.665758 1.019339 1013.370911 \n", "2 4.127301 4.127301 4.692918 4.692918 1.065441 1013.346924 \n", "3 4.168343 4.168343 4.756423 4.756423 1.177940 1013.292419 \n", "4 4.189304 4.189304 4.783813 4.783813 1.228021 1013.271912 \n", "\n", " zonalSurf SHAPE \n", "0 1.815351 {\"x\": -105.51129143475185, \"y\": 31.06161320108... \n", "1 1.796262 {\"x\": -105.39566505429877, \"y\": 31.06161320108... \n", "2 1.774962 {\"x\": -105.2800386738457, \"y\": 31.061613201088... \n", "3 1.715283 {\"x\": -105.16441229429194, \"y\": 31.06161320108... \n", "4 1.687406 {\"x\": -105.0487859138388, \"y\": 31.061613201088... \n", "\n", "[5 rows x 43 columns]" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Check the extracting data\n", "US_Prediction_Grid_GCM = os.path.join(outGDB, \"US_Prediction_Grid_GCM\") \n", "output_GCM_sdf = pd.DataFrame.spatial.from_featureclass(\n", " US_Prediction_Grid_GCM)\n", "output_GCM_sdf.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Explore the Relationship between variables " ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "station_GCM = os.path.join(outGDB, \"station_GCM\")\n", "station_GCM_sdf = pd.DataFrame.spatial.from_featureclass(station_GCM)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "gcm_variables = station_GCM_sdf[station_predictors]" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "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", "
airT2mgeoHgt850geoHgt850geoHgt500relHumSurfspecHum500specHum850vort850vort500merid850...merid500zonal850zonal500airFlow850airFlow500airFlowSurfvortSurfmeridSurfseaLevPzonalSurf
0289.849213164.443024164.443024164.44302484.12443510.73386810.733868-0.000187-0.0001871.560632...1.560632-0.228055-0.2280554.3007534.3007533.4807433.4807431.6411731019.5000610.168134
1288.776917159.583237159.583237159.58323781.95317110.15522110.1552210.0097560.0097561.874495...1.8744950.0811500.0811504.5692214.5692214.8042244.8042242.3062591018.8183590.031886
2288.779755161.611938161.611938161.61193882.47177910.26257010.2625700.0070050.0070051.700716...1.7007160.0112860.0112864.3801004.3801004.4022244.4022242.0130571019.1153560.132225
3292.476349165.336807165.336807165.33680783.30391711.52295411.522954-0.002438-0.0024382.187343...2.187343-0.947112-0.9471124.9866184.9866183.7407923.7407921.5715671019.343567-0.971350
4288.211273163.101761163.101761163.10176182.62401610.28277410.2827740.0056270.0056271.584904...1.5849040.0275540.0275544.1575634.1575634.3063434.3063431.7470681019.3632810.216421
\n", "

5 rows × 22 columns

\n", "
" ], "text/plain": [ " airT2m geoHgt850 geoHgt850 geoHgt500 relHumSurf specHum500 \\\n", "0 289.849213 164.443024 164.443024 164.443024 84.124435 10.733868 \n", "1 288.776917 159.583237 159.583237 159.583237 81.953171 10.155221 \n", "2 288.779755 161.611938 161.611938 161.611938 82.471779 10.262570 \n", "3 292.476349 165.336807 165.336807 165.336807 83.303917 11.522954 \n", "4 288.211273 163.101761 163.101761 163.101761 82.624016 10.282774 \n", "\n", " specHum850 vort850 vort500 merid850 ... merid500 zonal850 \\\n", "0 10.733868 -0.000187 -0.000187 1.560632 ... 1.560632 -0.228055 \n", "1 10.155221 0.009756 0.009756 1.874495 ... 1.874495 0.081150 \n", "2 10.262570 0.007005 0.007005 1.700716 ... 1.700716 0.011286 \n", "3 11.522954 -0.002438 -0.002438 2.187343 ... 2.187343 -0.947112 \n", "4 10.282774 0.005627 0.005627 1.584904 ... 1.584904 0.027554 \n", "\n", " zonal500 airFlow850 airFlow500 airFlowSurf vortSurf meridSurf \\\n", "0 -0.228055 4.300753 4.300753 3.480743 3.480743 1.641173 \n", "1 0.081150 4.569221 4.569221 4.804224 4.804224 2.306259 \n", "2 0.011286 4.380100 4.380100 4.402224 4.402224 2.013057 \n", "3 -0.947112 4.986618 4.986618 3.740792 3.740792 1.571567 \n", "4 0.027554 4.157563 4.157563 4.306343 4.306343 1.747068 \n", "\n", " seaLevP zonalSurf \n", "0 1019.500061 0.168134 \n", "1 1018.818359 0.031886 \n", "2 1019.115356 0.132225 \n", "3 1019.343567 -0.971350 \n", "4 1019.363281 0.216421 \n", "\n", "[5 rows x 22 columns]" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gcm_variables.head()" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "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", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
airT2mgeoHgt850geoHgt850geoHgt500relHumSurfspecHum500specHum850vort850vort500merid850...merid500zonal850zonal500airFlow850airFlow500airFlowSurfvortSurfmeridSurfseaLevPzonalSurf
airT2m1.0000000.5236390.5236390.5236390.0246020.4536660.453666-0.158646-0.1586460.115328...0.115328-0.511303-0.5113030.3472360.3472360.0156080.015608-0.0230690.447385-0.424286
geoHgt8500.5236391.0000001.0000001.0000000.2898580.4939860.493986-0.291146-0.291146-0.217218...-0.217218-0.486174-0.486174-0.025860-0.025860-0.486201-0.486201-0.2643050.988291-0.476785
geoHgt8500.5236391.0000001.0000001.0000000.2898580.4939860.493986-0.291146-0.291146-0.217218...-0.217218-0.486174-0.486174-0.025860-0.025860-0.486201-0.486201-0.2643050.988291-0.476785
geoHgt5000.5236391.0000001.0000001.0000000.2898580.4939860.493986-0.291146-0.291146-0.217218...-0.217218-0.486174-0.486174-0.025860-0.025860-0.486201-0.486201-0.2643050.988291-0.476785
relHumSurf0.0246020.2898580.2898580.2898581.0000000.6566110.656611-0.077643-0.0776430.327936...0.327936-0.486520-0.4865200.3783290.3783290.0796830.0796830.3550620.345988-0.514809
specHum5000.4536660.4939860.4939860.4939860.6566111.0000001.000000-0.002722-0.0027220.377871...0.377871-0.661489-0.6614890.2371640.237164-0.040895-0.0408950.3333650.511869-0.591803
specHum8500.4536660.4939860.4939860.4939860.6566111.0000001.000000-0.002722-0.0027220.377871...0.377871-0.661489-0.6614890.2371640.237164-0.040895-0.0408950.3333650.511869-0.591803
vort850-0.158646-0.291146-0.291146-0.291146-0.077643-0.002722-0.0027221.0000001.000000-0.190525...-0.1905250.4375140.4375140.0910730.0910730.1931340.193134-0.185647-0.2703640.541750
vort500-0.158646-0.291146-0.291146-0.291146-0.077643-0.002722-0.0027221.0000001.000000-0.190525...-0.1905250.4375140.4375140.0910730.0910730.1931340.193134-0.185647-0.2703640.541750
merid8500.115328-0.217218-0.217218-0.2172180.3279360.3778710.377871-0.190525-0.1905251.000000...1.000000-0.377805-0.3778050.5108180.5108180.6157990.6157990.929763-0.231771-0.304147
vort500-0.158646-0.291146-0.291146-0.291146-0.077643-0.002722-0.0027221.0000001.000000-0.190525...-0.1905250.4375140.4375140.0910730.0910730.1931340.193134-0.185647-0.2703640.541750
merid8500.115328-0.217218-0.217218-0.2172180.3279360.3778710.377871-0.190525-0.1905251.000000...1.000000-0.377805-0.3778050.5108180.5108180.6157990.6157990.929763-0.231771-0.304147
merid5000.115328-0.217218-0.217218-0.2172180.3279360.3778710.377871-0.190525-0.1905251.000000...1.000000-0.377805-0.3778050.5108180.5108180.6157990.6157990.929763-0.231771-0.304147
zonal850-0.511303-0.486174-0.486174-0.486174-0.486520-0.661489-0.6614890.4375140.437514-0.377805...-0.3778051.0000001.000000-0.254708-0.2547080.0549830.054983-0.283296-0.4594070.879255
zonal500-0.511303-0.486174-0.486174-0.486174-0.486520-0.661489-0.6614890.4375140.437514-0.377805...-0.3778051.0000001.000000-0.254708-0.2547080.0549830.054983-0.283296-0.4594070.879255
airFlow8500.347236-0.025860-0.025860-0.0258600.3783290.2371640.2371640.0910730.0910730.510818...0.510818-0.254708-0.2547081.0000001.0000000.6890660.6890660.403269-0.046564-0.163162
airFlow5000.347236-0.025860-0.025860-0.0258600.3783290.2371640.2371640.0910730.0910730.510818...0.510818-0.254708-0.2547081.0000001.0000000.6890660.6890660.403269-0.046564-0.163162
airFlowSurf0.015608-0.486201-0.486201-0.4862010.079683-0.040895-0.0408950.1931340.1931340.615799...0.6157990.0549830.0549830.6890660.6890661.0000001.0000000.631892-0.5120120.184889
vortSurf0.015608-0.486201-0.486201-0.4862010.079683-0.040895-0.0408950.1931340.1931340.615799...0.6157990.0549830.0549830.6890660.6890661.0000001.0000000.631892-0.5120120.184889
meridSurf-0.023069-0.264305-0.264305-0.2643050.3550620.3333650.333365-0.185647-0.1856470.929763...0.929763-0.283296-0.2832960.4032690.4032690.6318920.6318921.000000-0.265171-0.282589
seaLevP0.4473850.9882910.9882910.9882910.3459880.5118690.511869-0.270364-0.270364-0.231771...-0.231771-0.459407-0.459407-0.046564-0.046564-0.512012-0.512012-0.2651711.000000-0.463801
zonalSurf-0.424286-0.476785-0.476785-0.476785-0.514809-0.591803-0.5918030.5417500.541750-0.304147...-0.3041470.8792550.879255-0.163162-0.1631620.1848890.184889-0.282589-0.4638011.000000
\n", "

22 rows × 22 columns

\n", "
" ], "text/plain": [ " airT2m geoHgt850 geoHgt850 geoHgt500 relHumSurf \\\n", "airT2m 1.000000 0.523639 0.523639 0.523639 0.024602 \n", "geoHgt850 0.523639 1.000000 1.000000 1.000000 0.289858 \n", "geoHgt850 0.523639 1.000000 1.000000 1.000000 0.289858 \n", "geoHgt500 0.523639 1.000000 1.000000 1.000000 0.289858 \n", "relHumSurf 0.024602 0.289858 0.289858 0.289858 1.000000 \n", "specHum500 0.453666 0.493986 0.493986 0.493986 0.656611 \n", "specHum850 0.453666 0.493986 0.493986 0.493986 0.656611 \n", "vort850 -0.158646 -0.291146 -0.291146 -0.291146 -0.077643 \n", "vort500 -0.158646 -0.291146 -0.291146 -0.291146 -0.077643 \n", "merid850 0.115328 -0.217218 -0.217218 -0.217218 0.327936 \n", "vort500 -0.158646 -0.291146 -0.291146 -0.291146 -0.077643 \n", "merid850 0.115328 -0.217218 -0.217218 -0.217218 0.327936 \n", "merid500 0.115328 -0.217218 -0.217218 -0.217218 0.327936 \n", "zonal850 -0.511303 -0.486174 -0.486174 -0.486174 -0.486520 \n", "zonal500 -0.511303 -0.486174 -0.486174 -0.486174 -0.486520 \n", "airFlow850 0.347236 -0.025860 -0.025860 -0.025860 0.378329 \n", "airFlow500 0.347236 -0.025860 -0.025860 -0.025860 0.378329 \n", "airFlowSurf 0.015608 -0.486201 -0.486201 -0.486201 0.079683 \n", "vortSurf 0.015608 -0.486201 -0.486201 -0.486201 0.079683 \n", "meridSurf -0.023069 -0.264305 -0.264305 -0.264305 0.355062 \n", "seaLevP 0.447385 0.988291 0.988291 0.988291 0.345988 \n", "zonalSurf -0.424286 -0.476785 -0.476785 -0.476785 -0.514809 \n", "\n", " specHum500 specHum850 vort850 vort500 merid850 ... \\\n", "airT2m 0.453666 0.453666 -0.158646 -0.158646 0.115328 ... \n", "geoHgt850 0.493986 0.493986 -0.291146 -0.291146 -0.217218 ... \n", "geoHgt850 0.493986 0.493986 -0.291146 -0.291146 -0.217218 ... \n", "geoHgt500 0.493986 0.493986 -0.291146 -0.291146 -0.217218 ... \n", "relHumSurf 0.656611 0.656611 -0.077643 -0.077643 0.327936 ... \n", "specHum500 1.000000 1.000000 -0.002722 -0.002722 0.377871 ... \n", "specHum850 1.000000 1.000000 -0.002722 -0.002722 0.377871 ... \n", "vort850 -0.002722 -0.002722 1.000000 1.000000 -0.190525 ... \n", "vort500 -0.002722 -0.002722 1.000000 1.000000 -0.190525 ... \n", "merid850 0.377871 0.377871 -0.190525 -0.190525 1.000000 ... \n", "vort500 -0.002722 -0.002722 1.000000 1.000000 -0.190525 ... \n", "merid850 0.377871 0.377871 -0.190525 -0.190525 1.000000 ... \n", "merid500 0.377871 0.377871 -0.190525 -0.190525 1.000000 ... \n", "zonal850 -0.661489 -0.661489 0.437514 0.437514 -0.377805 ... \n", "zonal500 -0.661489 -0.661489 0.437514 0.437514 -0.377805 ... \n", "airFlow850 0.237164 0.237164 0.091073 0.091073 0.510818 ... \n", "airFlow500 0.237164 0.237164 0.091073 0.091073 0.510818 ... \n", "airFlowSurf -0.040895 -0.040895 0.193134 0.193134 0.615799 ... \n", "vortSurf -0.040895 -0.040895 0.193134 0.193134 0.615799 ... \n", "meridSurf 0.333365 0.333365 -0.185647 -0.185647 0.929763 ... \n", "seaLevP 0.511869 0.511869 -0.270364 -0.270364 -0.231771 ... \n", "zonalSurf -0.591803 -0.591803 0.541750 0.541750 -0.304147 ... \n", "\n", " merid500 zonal850 zonal500 airFlow850 airFlow500 \\\n", "airT2m 0.115328 -0.511303 -0.511303 0.347236 0.347236 \n", "geoHgt850 -0.217218 -0.486174 -0.486174 -0.025860 -0.025860 \n", "geoHgt850 -0.217218 -0.486174 -0.486174 -0.025860 -0.025860 \n", "geoHgt500 -0.217218 -0.486174 -0.486174 -0.025860 -0.025860 \n", "relHumSurf 0.327936 -0.486520 -0.486520 0.378329 0.378329 \n", "specHum500 0.377871 -0.661489 -0.661489 0.237164 0.237164 \n", "specHum850 0.377871 -0.661489 -0.661489 0.237164 0.237164 \n", "vort850 -0.190525 0.437514 0.437514 0.091073 0.091073 \n", "vort500 -0.190525 0.437514 0.437514 0.091073 0.091073 \n", "merid850 1.000000 -0.377805 -0.377805 0.510818 0.510818 \n", "vort500 -0.190525 0.437514 0.437514 0.091073 0.091073 \n", "merid850 1.000000 -0.377805 -0.377805 0.510818 0.510818 \n", "merid500 1.000000 -0.377805 -0.377805 0.510818 0.510818 \n", "zonal850 -0.377805 1.000000 1.000000 -0.254708 -0.254708 \n", "zonal500 -0.377805 1.000000 1.000000 -0.254708 -0.254708 \n", "airFlow850 0.510818 -0.254708 -0.254708 1.000000 1.000000 \n", "airFlow500 0.510818 -0.254708 -0.254708 1.000000 1.000000 \n", "airFlowSurf 0.615799 0.054983 0.054983 0.689066 0.689066 \n", "vortSurf 0.615799 0.054983 0.054983 0.689066 0.689066 \n", "meridSurf 0.929763 -0.283296 -0.283296 0.403269 0.403269 \n", "seaLevP -0.231771 -0.459407 -0.459407 -0.046564 -0.046564 \n", "zonalSurf -0.304147 0.879255 0.879255 -0.163162 -0.163162 \n", "\n", " airFlowSurf vortSurf meridSurf seaLevP zonalSurf \n", "airT2m 0.015608 0.015608 -0.023069 0.447385 -0.424286 \n", "geoHgt850 -0.486201 -0.486201 -0.264305 0.988291 -0.476785 \n", "geoHgt850 -0.486201 -0.486201 -0.264305 0.988291 -0.476785 \n", "geoHgt500 -0.486201 -0.486201 -0.264305 0.988291 -0.476785 \n", "relHumSurf 0.079683 0.079683 0.355062 0.345988 -0.514809 \n", "specHum500 -0.040895 -0.040895 0.333365 0.511869 -0.591803 \n", "specHum850 -0.040895 -0.040895 0.333365 0.511869 -0.591803 \n", "vort850 0.193134 0.193134 -0.185647 -0.270364 0.541750 \n", "vort500 0.193134 0.193134 -0.185647 -0.270364 0.541750 \n", "merid850 0.615799 0.615799 0.929763 -0.231771 -0.304147 \n", "vort500 0.193134 0.193134 -0.185647 -0.270364 0.541750 \n", "merid850 0.615799 0.615799 0.929763 -0.231771 -0.304147 \n", "merid500 0.615799 0.615799 0.929763 -0.231771 -0.304147 \n", "zonal850 0.054983 0.054983 -0.283296 -0.459407 0.879255 \n", "zonal500 0.054983 0.054983 -0.283296 -0.459407 0.879255 \n", "airFlow850 0.689066 0.689066 0.403269 -0.046564 -0.163162 \n", "airFlow500 0.689066 0.689066 0.403269 -0.046564 -0.163162 \n", "airFlowSurf 1.000000 1.000000 0.631892 -0.512012 0.184889 \n", "vortSurf 1.000000 1.000000 0.631892 -0.512012 0.184889 \n", "meridSurf 0.631892 0.631892 1.000000 -0.265171 -0.282589 \n", "seaLevP -0.512012 -0.512012 -0.265171 1.000000 -0.463801 \n", "zonalSurf 0.184889 0.184889 -0.282589 -0.463801 1.000000 \n", "\n", "[22 rows x 22 columns]" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Let's create a correlation martix to check the relationship\n", "cor_table=gcm_variables.corr()\n", "cor_table" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1QAAAMFCAYAAABzs7cPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3Xe8XFW9/vHPExJAEnpvEqSKSlFEsIGoiCiIoj8FVIqAXguWi9gQEfUq1+61glKliREBRRAVBAGRiFRRQAldIQRCCBAp398fax2yM5kTDiez9jqZed6v13kxs6c8ex+yz8zaa63vUkRgZmZmZmZmT9+42jtgZmZmZma2qHKDyszMzMzMbJTcoDIzMzMzMxslN6jMzMzMzMxGyQ0qMzMzMzOzUXKDyszMzMzMbJTcoDIzMzMzMxslN6jMzMzMzMxGaXztHTCzRZ+kicCGwNLALOCGiJhdd6/MzMzMynMPlVkfkHSApEskzZT0eP7vJZL2L5y7nKSTgRnAhcBJwO+BeyWdKGm5kvlmZmZmtblBZbaIk3QE8EHgh8D2wEbAK/L9D0r6YsH4o4EngGdHxNIRsVZELANskrcfXTDbzMzMrDpFRO19MLOFIOkeYNOIuKvLY2sAV0fESoWyZwGrRsRDXR6bCPwrIpYukW1mZmY2FriHymzRp4V8fGHcC2wxzGObk4YCmpmZmfUtF6UwW/T9CPidpK8CVwEzgWWAzYCPAEcVzP4k8CtJZ3bJ3hl4T8FsMzMzs+o85M+sD0h6N/BO4DnAJOBB4Drg+Ij4QeHs5wB7dMk+MSL+WjLbzMzMrDY3qMzMzMzMzEbJQ/7M+oCk8aQKf89h7lpQ1wLnR8RjLe7HusBOpHlb50bEjW1lm5mZmdXgHiqzRZykzYAzSI2Yq5k7j2lTIIA3RMTVhbKvj4hn59vbAmcBF+fcl+Xs35XINjMzMxsL3KAyW8RJupw0V+r/ujz2fmDviNiyUPasobLoki4CjoqI4/P9PYH3RcSLS2SbmZmZjQUum27WA5JObdzep+X4TYDvD/PYkcCzC2Y3r8hsBJzcuH8KsHHBbDOzRU7lzwszK8ANKrPeeI2kofWevtly9vXAfw3z2Lvz46VMkLSPpH1JjavFG4+NBxYrmG1mtiiq+XlhZgW4KEWfk7QisBvzFiu4DpgSEff2cfaGnbkRcUPByIuASyXdACwp6fhuT4qIdxbI3g/4uaSPMv8cqseBXQtkDrmMVK4d4K+k3rLL8/1tgb8XzDazPiNpefLf7Yi4r/b+FFLz88LMCvAcqj4m6ZXAT4FrmH/R1ecBu0XE+f2ULemZwKk55x+N3PXyfrwtIm4tkLsk8GZgHeBQ4H+6PS8iPtvr7Jw/AdiO+deCuiAiHi2ROYJ9WhaYEBHTa+Sb2aIh//36LLAPsAqpwE4A/waOAQ4r9XdM0g7A3sx/4e+YiDivUGbVzwsz6z03qPqYpL8Ch0TEz7o89kbgf4YqtPVLtqTfAn8mfQA/1Ng+kfTB9cKI2L7XuY2cccDvgVdGxH9K5Yx1klaIiBm198PMxj5JPyRd9Poc816A2xz4FHBzROxXIPfDwMHAUcx/4W9/4IiIKDYkz58XZv3DDao+Jmk2sEJEzOny2BLAfRGxVD9lS3ow58734ZRzZ0TExF7nduQ8WfmuDZImR8S0xv23kq5+Cvh5RPy4YPYawP+RhvqdTfpCdA6wFXAH8MaImFoq38wWfZLuAyZHxMwujy1PalAtVyD3TmD7iPhbl8c2Jq3jt3qvcztyWv28MLMyXJSiv10GfD73zjwp3/9cfrzfsm8DXj/MYzsBPR/u18WFkrZuIWfIk2tMSXoP8HVgKvAn4EuS3lcw+3vAbOC/gdWA3wK/AJbPj32tYLaZ9YdHgOEaLqvlx0uYCNw5zGP/AopccOzQ9ueFmRXgHqo+JmkdUhnrLYB/Mnc4w7OAKyk0n6hmdp67NQW4lvmHjjyHNHer6EKzkr4L7E5abPc2GqXFI+LQAnnNtaCuAQ6IiEvz/a2AYyNik17n5vefDqwREf+RtBRwPzAxIh7N8yL+HRErlMg2s/7QGHr3I+Yfevcu4H8j4hsFco8mDTU8nHkL+mwGHAL8MyLe1evcjn1o9fPCzMpwg2oASNqAjmIFEXFjv2bn6oJv6swFTm+jQIKkY4Z7LCJ6vuaIpAciYpl8ezqwcjRO7ObjBbL/BTw3IqZLWp00zG/5iJgpaRngxohYtUS2mfUPSa8hVQzt/Lt9fEScWyhzcVIxjHcAa5AaMwLuAo4HPlN6blPbnxdmVoYbVGaLOEn/IS2iC7Az8OyI+Fd+bDngpohYqVD294EXAecCLyUNqVwCOJb05WhORLy9RLaZWa/kv5WTgAcj4v7a+9NvJJ0aEW/Nt/eJiGEbkmaLIjeo+pikl0fEhfn2OOCjNIoVAF+MiMf7MLv1Mrgd+c8a7rGI+GeBvM90bDp1aJK1pF1Iwyv36HVufv/FgQ8Bk0nDdW4Avk9aA+ty4CBX+zOzpyP3br+M9Hlxcen1qLqsW3hti6M4Wv28qEXS/aTRC1Fy1IRZLW5Q9bGOoWCHkMZpH54fPgT4acF1kapk1y6Dm/fhCeYOHRkSABGxWMlsGwwDsvjpPAbxmAeFpHMiYsd8e1NSj/f9pL+bKwI7RsRfCuRWWbewYx8G4vNC0lnAyqSLbm9j7qiKefTLYsaSvhwRH823ty89d9vqc4Oqj3UUK/gbqSDDdfn+xsAvImL9fsoeC2Vwu+SuBnwGuCgiTiqctRSwPmnoypMi4pKSubWzB0HNxU9rqX3MtXu7B0XHBbhzSAuSfynf/yiwQ0S8ukBu1XULh9mn1j4v2jRoixlLmhkRy+bb7pEbAG5Q9bGOD6m7I2KVjseLrX9RK1vSTGDtiHigy2PLAbcM/ZFrU14D64aIWKdgxjuBbwP/AR5uPBQR8cxSubWzB0WtxU9rqnnMY6G3e1B0fF78G1gnIh7J95cAbo+IlQvkVl+3cJj9Kv55UcugLGYs6ULgbuCvwMeAI7o9z5Uc+4cbVH1M0mPAUO/AFqRqbLfkx1YBro6I1fopeyyUwR1mvzYFflviS0Ej41/AO2pcOa+ZPShqLX5aU81jrtnbLekA5vaMNSveHRMRR5XIzLnPJlW86+yROyEiri+Y+xDwOlIP5Imk/+dz8mMTSA2bEhfgrgc+FRE/6/LYG4H/iYhn9zp3BPtV/POippIXc8eK/D3nAFKP3F6kf9edIiL2bXXHrJjxtXfAiupsODTHaD8fOKEPs99DGiJ0HHPL4EIqg3sCaShFUZIuauRCWhzyOcydQ1bKf4ALCmeMxexBMbT46XyNC8ouflpTzWOusuirpCNIi5N/lfl75T4i6VkR8YkCubuTFuM+E7iQeS9GXSLpPRFxaq9zs7uBo/PtR0ifEZfm+y8krWVYwvuBKZI+wjDrFhbKfVLFz4uaLpS0dUT8sfaOlBIRdwOflyRS5dv9ShXisrHBPVTWt2qVwZW0V8em2cBVpatG5dwtgc9GC+ttjZXsQaFKi5/WVPOYa/V2S7oH2DQi7ury2Bqk3v2eL4Mg6Wbg7RFxcZfHXgKcGBGTe507gv1aG1g6Iv5a6P1rr1tY5fOiJg3YYsaSZpP+DT9Re1+sHDeo+pyklUhfCq6OiIfzUIJXkf5g/7Zfszv2YxLp3/qstjJrkLQNqXLSWs3NpGEFRatF1cweJKqw+GlttY5ZlRZ9VVqc+3kLaFBdExErFsh9kLQo+MNdHlsKuDsiJs3/SrOnRwO2mLGkP5B6qOYbPmz9ww2qPibpDcBJpCtejwDvJl3pvYy0xsfnC17drZIt6VMR8YV8e0XSuOUdSF+Gzgf2yF3xPSdpR+CBoap2ktYjffF6Lmn4yj7dviT1MP8m4GRSGeB5vhRFxD9K5dbONiutzd7uPORvF+Yf8rcZ8BHgrIj4eIHc00l/qw9pnrP579jhwFIR8cZe5+aMScBhpEbz2aShhycCryX9DvYu8XdEkqLxJShfGHoTqeF8RkRc1OvMRlbVz4tBJGmFqLAuoqTPA28nLXjf2SN39DAvs0WMG1R9TNK1wMERcbakXYEfA9tFxFRJWwCnRbmy6VWyO6pFHU2aWP3+/PA3gUciYu9e5+a8y4EDI+LSfP/3pAbld4B9c/aeJbJz3n2kilWtn9Q1sweVWl78dCyoccyqsOirpHczfK/cDwplLg98l9SgeIy5DbnxwM+A95X6fUs6gTQn7nTSvKXHSPOmvk+6GDc5InYqkNv8vHgT6QvviaQvvHuSehVO63Vuzmv980LSqRHx1nx7n4gYtqeoJLW8mLGkrUn/hlcFbgd2jQLrmi0g//xhHopouSy/leMGVR/TvOsgjAPmRMSEbo/3S7bmXf/qVuAFEXFPvr8iafjhmr3Oze8/A1glIh7LFX7uJJX/vSMPf7w6ItYokZ3zvwZcGRHHl8oYi9mDQpUWP62p5jFrDCz6WkMe3rchcxtyN0RjjaZCmXcD60XErNywuwdYJiIekvQM4LZC88aanxeXA4dGxK/y/R1IpfG36HVufv/WPy8k3Q8sHxGhimsjqeXFjPOQu4tJxar2AbaIiFf1OscGm6v89bc7JO0QEb8mlaR9RNIWEfEXSZuRKiv1W3bkqjrjSH+s7208NoP0hagN25BKOt+R799Lx4K3BWwFvF/Sp0gLnz4pIl7er9n5C9jSwKw+76V5ceP2/wLfjHkXP/1foOeLn1ZW85iPAS4irZfTbdHXY4HiV5fzBalnk/6eXR/lK4Up/4xr3C5NwNCE/ehyv5Tme69LarAPOQ/4ScHsprY+Ly4CLpV0A7CkpK4XwCLinQWym+8/rnlfjcWMC0VuAmwbEY9LOgS4pVBOV/kc7ipcqKJvuEHV3w4DzsrDsa4HPgicl4cWvCw/3m/Zk0jDRUT6sNwcuCI/tgHpymcplwMHKi1Guh/wq8ZjzwJKV4w6Kv/U0Gq20to0nyVdbVyF/P9baVHQY4DDIuLRQtk7MHeNoOZaPcdEe+twbQHs2rj/LVI1vCIG8ZiBFwGv7Sw8ERGzJR1KukDTc5KOjIgD8u11gF+QvuwHcLuk1xeaT7Qcae7Sm0jLIAz1yC0uaQppyF+p+WO/I5Uv/zmwc77/FUk/IP0tvWRBL14IS0hqlidfFhi6KLMUZRtzNT4v3gK8mbQ2UpB6XquLiH9J+hBwA2nuda+NH7oQERFzcsGZNj3G8P+WXLSpT3jIX5/LVaHWAqZGxBOSXkUawnJ5RFzYb9n5C0jT9IiYnR/bCnhWRJxSKHsT4CxgMnAT8IqIuDM/dhhpHsDeJbIHTf4Ssh7wOeZfP+ZTpKu9+xXIHSrjfRTzFwzYnzRE6Ju9zs3ZtRY/Hbhjzu9fZdHXjnk9PyVVFfwg6XfwVeCZEfGmArk/IxWU+XRzHkue7/JZYGKJ3JyxPPAF0t/O75Augv0U2BT4E7Bv5IXhe5zbOYfomxFxZX5sN+CDpXrYG58X6wI30uLnRe4x+T2p97XnlSpHQwUXM5b0H1IV2iFvJQ3nfVLJXrku30tWBz5OKjDzo1K51i43qAaApLd0m1gr6c0R8dN+zK58zCtGxL0d25YD/lNyLoKkYVdcj8KVhNrOzj2fkyNivgVf85ezmyNiuQK5dwLbR5fyt5I2Bs6PiNV7nZvffxrzXuXcozGh/cXA9yJiswK5A3fM+f1fCUwBrmWYRV8j4ncFcpvzeu4ENo6IB/L9ScBNEbFaiVxg1W5/o/Iwx3+VarwOsoqfF7Nq/f/UAhYzjogvFsj7zFM9JyI+2+vcBZG0LOni8oZt5lo5blANgOEmn0qaEREr9GN2zWNuZC1Dxzj4oSuQhfI6KwmtRurFuTgiXlEqt0a2pLtIV3S7fcl/NulLfokvnTOBtYe+4HY8thxwSxQq9PIU+1Vs8dNBPOZGRuuLvuaGzaakOUy/J104eCw/Ng6YWagnchqwZwy/sO9JEdF5pb3X+3AgcEFEXF0yZ6zk1syW9EvgcxHxxzZzc/bALWbcKf/9ujoilq+9L9YbnkPVIknjSauDb8H8X7QPKJA3VJp0nKR1mXdy8bNIa44UUSu75jE39uFVwJGkceqdVYyKjZfu1nDJPUc9H5Y0BrL/Fzhf0o+Yfxjau4AvFcqdQpobeDhwdUfuIaRhSq2LiNsKvv0gHvNQxr20Py9xImnI8NDfjq2BP+TbzwPu6PaiHvgk8CtJZzL/ObUz8J5CuU0vBP5b0tKkAgW/zz9/KTx5v1ZuzexbSP+/z2D+tZEOLZhLRBxX8v2HI+lK0tDhkyPi9hZzT2D+HrmXk5aTsT7hHqoWSTqF9IH4K+Zf/PTTBfK6lSYd8i/SxP0je51bM7vmMTf24RbS3J5TmP//c+kKXZ37Mo40j6yVXrk2syW9huHX6jl3Qa9diMzFSfNJ3gGswdx/a3eRFuX8TKk5Caq3+OnAHXPOrrbo6wL2aRPSsLzh1rVZ2Pd/DrAH859TJ5bsCeyyH5OBbfPPmwBKDOEdK7k1srvMH3tSROxTKLPqYsZ57uPuwE7An0kFME6Lwov9dhlyOJu0zMhvSuZau9ygapHSGhBrR8SslnN/HxHbtplZO7vyMf8bWKNS46lpKdLq7AdHxLALKS7q2bXk4W6TgAejXPWzZl6VxU879mFgjlkVF30dZJI2IjUqtgNeQlqG4fcR8dF+zK2d3SZVWMx4mP1YmtRo3Z1Udfi3EbFL6Vzrb25QtUjSJcDuUaBakY0dkj5Ouor9pWjxBGv0zjXdAexfqsdmLGTn/GVIH4wizdsquh6VpA2Zt4T4taXH/6vS4qeN/IE6ZlVa9DVfnHgv6Xf9q4g4U9IRpF65K4GPlJq/1bEf65Ku5As4t435Lfli1CzSMNILSOdy8QuQtXJrZjeGx88nGlUee5zZ+mLGC9iXCaR/3weS1qjq+RQYSS8A5kTEtfn+ysA3mNsjd1BEPNjrXKvDDaoW5T9gPwB+zfyLn3ZdYK9HuZ0VdYbMAW4HfhYRZ/VTduVj3oC0QORKdKwlUrK3RvOXZp3dxhevGtmSzomIHfPtTUm/7/tJ/89XBHaMiL8UyH0mqdzuZqQ1XIbmmaxHGoL2toi4tde5OfseUo/M7NxTdC+plPUjkpYEbi/RuBjEY87ZzR6q6aQvgk/k+wLuKzEkS9L/kXorziE1oi4HViCtr7YX8GhEvLVA7vWRy8BL2pZU0vti0jn1MuANUaCqYcc+/BB4ab47NJfowtLzXWrl1sweZnh8AEREkbm+HQ2qNwBfiYgN8mMiFVyZr5hUD/NFWox7D+CNpHlkJwGnlJiTmb+HfHZoaF+er7YGqbd7d1ID8r29zrU63KBqUZ7UfTCpDG9zbk1EobUucu7nSB/Ex5Emn65NmntyEumP6buAL0fE//ZLduVjvop0Jfk05p9D9dte5w2iji+755CqZH0p3/8osENEvLpA7m9JY+8Pi0ZJY6Wy0ocCL4yI7XudmzNOJS08OrT46eLA30kXafYD1i0xbGUQjzlnzwGOyHffC2ww1POZj/32KFChS6lU+uYRcbekNYFbgZUi4r7cqLwhIlYpkNvskbsIOGroQp+kPUkL+76417nD7MuqpEn725KGDk+PiPX7Nbd2ds5fDfgMcFFElFhcF0nnki5+/ZA0fPbmiDgwP7YecF7hi453keYFnkKqWnl9qaycNx1YM9JiwssBdwPPjYgblKr8XRIRa5fcB2uPG1QtUio/vHXpk7hL7mWkydvXN7ZtDBwXES9SWvD2lBJ/yGplVz7mB4Dlonx1qKG8zgpC84lCixbWyu5oUP2bNGzkkXx/CdKX3RILRD4IrBBdijDk3BkRMbHXufn9ay1+OnDHnLOrLPqar+KvGhGP5mGNDwBL5fuLAfdEgUIvHefU3aQvgo/m+8Vyu+zHFqRGxStIPWMPk+YT7dGPubWzO/ZjCVKDvUh5fM1dzHgyqZJla4sZ54wXRcRlpd6/S979wPIREUoFOY6MiGc2Hq+2Fpj1nsumt+vfpKuNbduYNJG76RZgI4CI+FMez9xP2TWP+QzSsIK2Kvjc1HH/Y8y9st6v2eMlvYLU29g5f+sJYMlCubcBrwd+1uWxnSh4fufekc7hIS8pldcwiMe8wEpnETGFVE6+hEuBH0j6CWlY0FWkstrfAf4r3y9hgqR9SOdUkHoDH82Pjafgkg9DlBbsnglcCJwJ/HdEdP6N6Zvc2tldbEQqKFREpEqR66nLYsakuUVFqoU28i9TWqfwzcBqEfG+fKF18SizDth1wFuAnwBvo/GdIPc+z7cwvS263EPVIkn/BbyGtEbO3c3HotAk0Jx7FmnS66Gk+UNrkUoRLxcRr5f0PNKcog36JbvyMf+E9AX0IuafK1ekp6gj/74SQ5HGUrbSIqTNP157xNzKUS8GvhcRmxXIfSXpi/S1zLtWz+akIgK7lZhnovmrKHZVold0QI+5ZvY6wHeBdUlfMi8kDZNaC7gZeFOJL3+SLmDec+rgiLg8P7YD8PmI2KrXuR37MDkippXMGEu5NbM1/zzjpUjn8+ER8cUCedXOqcY+vIV0bk0hfWYsI2lLUgGpVxXIeympRy6Ax4GXRsTf82MfAV4UBeZDWh1uULUoTwLtJqLQJNCcuwLpj8ibSFcZHyNdbf5ARExXKtm6dERM7ZfsysfcuebEkyLis73O65I/o42hOWMtu7EPa5P+3xZZN0fSiqR/V51r9ZwehQpxqHsVxfmU+jsyaMdc+/fdSZJIwy47r+q3QtKywIRS/687sjYg9cytSaoUekpE3NCvubWyJe3VsWk2cFUUquY4Fs4pSdeTKi1fOXTxT6na350lhojnzKWBDUlDKWc1tm8EzBoa8miLPjeoBki+QrQyaSx8K/N7amfXPOZaBr1B1Y80bxXF15GGrHyRNIx1HdJQyykR8b0Ku1dEzWMexN/3WCBpZ1Kxgl+QftfPJPX2vyMizuy33NrZbRoL55Ske0kFXmLos0rSeFKDqtQUABsQblD1sZF0sZdqZNTKrnnMjX0YtuJZoaFRnYU1rgC2gLnlcEsNKa2VLWkSaQjnc4Czge+RvpS8ljQsbe+I+EeBXEXjj6akbUg9NwLOiIiLep3ZZR9uAraMxsK6uXjD1IhYr0DewB1zrWxJtzGyq/jPfKrnLAq5XfbjGtLCr+c3tm0HfDsinttvuTWyc3GEByLiknx/PeB45q6NtE9E3NXr3I59qHI+S/o18OOIOL7RoHo7aemH1/c4a0ycU9YeN6haIun/kSZTXwccE7l6Un7su1FgLYKn6GIXBYca1squecyNfbi5Y9PKpAnet0eZqoLd1hNpauP/c6vZStUFVwNOB3YjDen8J/B94N2kalE7FchtVkJ7E2k9kRNJv4M9gf0i4rRe53bswz3AZs2hInmC81VRZh2qgTvmWtlK6z89pYj4fT/kdtmP+4CVI+KxxrbxpBLiPV/zq3ZujWylRaoPjLlzTn9PGu73HWBf4JGI2LPXuR37UOV8VipA8WvSXMStSQspb0haZqOnQx3Hyjll7XGDqgWSDgLeT6r+9jLSl8+dhq4CNb+w9Di32cUu0pX7TZvPiXIlh6tk1zzmBezTYsAhpPHSX2szu18plXVeLyJm5Sub9wDLRMRDSqWmbyvUuGiu1XM5cGhE/Crf3wE4IiK26HVuxz58hVRd7xvMXWPtQODciPjvAnkDd8xjJXvQSDofOCcijmhsO5j0ebldv+XWyNa8i+uuAtxJWnbiDkkrkRabXaPXuR37UPN8Xoo0pHKdnP2LiHiwZKYNBjeoWiDpH8BrI08ylfRZ0sJ920fELWppLYJBnFszVub05CuOt0fEagUz1gQeirz4aN62PPCMKDzxte3sfIVzckTMVlow8V5gYkQ8ImlJ0u+6dG/NdNIXkyfyfQH3tXBFexxwAKkc7xrAXaSyvEdFxOMF8gbumMdQ9uaki3ArMe8w2kP7NHdjUlW0icz9oj0b2CUKrt9YK7dGdkeD6g3AVyJXu83n88wSF3g79qHaOdVlXyYAv40Ca8p15FQ5p6w9blC1QF0WepX0fuBg4NXA5aX/gOVMN6gqkfRa4Eclr/zl3oN9I+KaxrbnAT+MiBeVyq2RLelUYFng58DOpCGVfwd+AOwHrBsRuxTIncPcdbbeC2ww1IiUNJHUkKtSsr6UQTzmsUDSAcDXSUOUXgv8CtiBNG+t2IKvtXIb+eOBbYDVSb0nlzWHyPdbbtvZks4lleH/IWno7s0RcWB+bD3gvBJD08cqpcWMHyo5FaD2OWXt8MK+7biFNOzsyqENEfFtSQ+RxvAuUWm/rIAuk1GXIi002/N5ch02bDZoACLimnwFtLS2s98DfAHYhTT2/wrgp8AlwJ+AdxXKPYl0BRnSVeV1gKFeuR2Ba7q9aGFJekdEnJBv7zvc8yLi6ALxA3fMlX/fQw4GdoyIi5RKPL8xX5h5W8HMmrkA5LlExQudjJXcCtkfJp3HXwFuJM05HfIO0rpnPTdGzqlaqp5T1g73ULUgz6EiIr7S5bE9gc+VuCIk6fCOTR8FvtzcUKq7uVZ2zWNu7EPnZNTZpDUoHiicexPpj/ZNjW3rA78ufcWxZvYgkHR25CIbec5FNxERw1aYXNTUPOax8PvuGGp5L6lwwROle93bzu1yAaqr6KPKhjWzG/uwYnSsa5aHT8+JiIcL5FU/p7ppqYeqyrls7XIPVTsujYiLuz0QESeSut1LWLvj/ikd20q2pmtl1zzmFFCvas/RwBRJnyJVvFsP+BxpaEdfZkv6S3QpiCBpakRs2S+5jS8iIvW+3RqNqmBtGKRjHgu/b+B2SZMjYhpwA/CGPIftP32W+/ZC7ztWc2tnD/kNaYmLJ0XE/ZKmAn11Pne50NrUxvfgWueytcg9VC1QoSp+NrY8xR9toGzvWJ7o+9+kD6u1gVuBHwFfi/Jrb1XJVpeCLvkD+97CV/Gr5Oac2cDSpf+fdskduGOumS1pb+DfEfGrPDzop6S5ggdG2cVPq+RauwbpfJZ0zFM9JyL2KZi/Nz6n+p4bVC3o9oerwj5sAOwOrAncAZwcPV53Yaxlt53b5Y/27sDJjfsREcOOHbeRk3R8vvlW4NSOhyeT/ra9rF9+RAe1AAAgAElEQVRyO/bhD6S1n/5WMqeRN3DHPFayO/ZjcWDxaLnEcxu5nUOfJH0jIj5UKq92bs3sQT+fx4Ja57KV5SF/7ZCkdWHYxU+JiH8WDN+ZNKzwF6QCGRsBU/Mk0TNL5dbMrpHbeYVL0q4lr3p1I+nVpImuq0TEzpK2JK3P9Ls+y/7HMLcDuBgotdBsrdymC4BzJB1LKrP85FWxQhO6B/GYx0o2SmsFTWreL/l5USl3Qsf9dwJtNGxq5dbMHvTzGUlLM3/58n47p6xl7qFqgaQnSH8whmtQReEJkdeQupbPb2zbDvh2RDy3VG7N7JrH3MhrdcKppA8AHyTNW/pERCwr6TmktT1e3G/ZSgsmv4dUmn1OiYyxlNvIb31C9yAec+1sSTuShs2uxryfHaU/L1rP7RwWr1QJrXgp/lq5YyB7IM9nSZuQLrRuxtzvZJGD++qcsva5QdWC2kP+JN1HqirzWGPbeGB6lF+Qs0p2zWNu5LXdoPoH8MqImDb04Zw/OO+OiBX7MVvS/W39/xwLuTUN4jHXlM+pLwPHRYGqa2MpV9Is0tIiQ182ryAVTCjag1Art3Z2zh+481nSBaTf8+HAzaQhjl8ELomIHxfMrXIuW7vcoGpB7aIU+WrQORFxRGPbwcBOEbFdP2bXyFUqzPDkXeAeYEXm/YAsWaDhbmD1iHh8qDEnaUnSwo2rl8qtmS3pBOAnEXFWqYyxlNvIXxHYCVgtIr4saQ1gXETcXjBz4I65ZrakGcCK0fKHdI3cWqM4ao4eGQMjVwbufM4XWleJiEeHGpRKi5NfGxHrFsytci5bu9ygasEY6KHamLSQ30TSeOW1SWsj7RIR1/djdo3cxgfkk5sa90X5D8ifAn+JiC80GjUHA5tH4dXYa2VLOo20uO+lzD8W/539lpuztwWmAFOBl0TE0nnbQRGxc8HcgTvmmtmSvgxcHy0vdFor19o1iOezpLuA9SLiIaW1E7cnLVJ+R8mL3j6nBoMbVC2QtDZwJ2n9gU0qjVkeD2wNrJH35bKIeLSfs9vOlbTOUz0nIm4pmL86qRG5Eqmy4T+BB4CdI+JfpXJrZkv6zHCPRcRn+y03Z/+F9MXjt43hlUsCt0TEqgVzB+6Ya2ZLugjYilRUZ55zKCJe3oe5i1HhM7JW7hjIHrjzWdJPgLMj4lhJXwJ2BuaQ1sTatWBulXPK2uUGVYsk3QC8MCJmVsieQG5cRMSpuZubiJjdr9k1j7kWSQJeCKxDuur4p5LDDMdK9iBRY/J6ozdwHHBPFJ4rV0vNY66VLWmv4R6LiOP6LTdn3wBsFRH3l8wZK7m1s2sZC3/Dct4ewNLA8SW/F9Q8p6w9Lpverm8AP5H0P8DtzNvFXnLy6fOAM0lXYtYirT2xLbAXaS2KYmplt52bx6M/5dWJ0kMpgHHMLce7GMOPz++bbEmvAN7B3PXGfhztlImvkgv8VdJrIuLcxrZXAdeUDh7EY66VXeuLVuUveN8ATm37M7JibtXsAT2fhxpTq0bBQhRNbjQNBvdQtSjPsemm9NyaPwA/iIgTGt3rE4EbImLNUrk1s9vOXdDwiabCw9A2BX4OLEH6cFwLeAR4Y0RcVSq3Zrak/YD/IZVrvwV4JvAu4NMRcVS/5ebsrUnrq/0S+H/A8aS5ELtExOUFcwfumMdA9j7M+4X3hIjoXEC8n3JrfUZWya2ZPYjns6Tlge8AbwYejYiJknYh9RAeUio3Z1c5p6w9blANAKXKNitERKhRylstlPWulV3zmGuRNBU4GfhaPm4BHwb2jIgX9GN2Hi7zlmajLTfupkTEBv2W28haA3g7aXjlraQry3cUzhy4Y66ZLelTpMVev0r6wrsO6Zz6cUR8od9yrV2DeD5LOoVUhOJw4K/5QuvKpLLpJT8vfE4NADeoBkCeALp/RExtjFfeirTI7Vb9mF3zmHP+q4G3kUq07ixpS2CZksMpJD0ALB8Rjze2LQbcF4XL9tfKlnQvqezuo41tSwB3Fp7fUiU35ywLHAg8H5jEvEOEdiiYO3DHXDNb0s3AdtEoZKNU+ObCiHjKAjiLWq61axDPZ0n3kOZUP9pxoXVmRCxbMNfn1ADwHKrCJJ0TETvm2xcxzDybKFvp5dPALyV9H1hC0ieA/wL2K5hZO7vaMUv6APBB0lCKN+fNDwPfAl5cMPps0rCJ0xvbdiYNqyitVvbFwNckfSxSKdyJ5IUa+zQX4DTSHLXTSf+u2jKIx1wzeyJpLbume4Fn9FNurc/Imp/NY+R7wSCezzNJlWjvGtog6ZnN+4XUOpetRW5QlXd84/YPa+xARPxC0o7A/sAFpLHSb4yIP/drds1jBj4EvDIipkn6WN72N2CjwrmLAadI+jNz1956AXCGpCf/HUaZwhi1st9NGmr4QL7iugLpC8HuPc4ZK7mQKleuGC0te9AwiMdcM/sc4ERJHycNiVoH+AJw7gJftejl1vqMrPnZXP17AYN5Pv8QmJKH4I2TtA1pHtn3C+fWOpetRR7y1zJJq5LWI1iJRhW0KLjgm6TDh3loDqmq0DkR8e9+yq58zHcDq0fE443hhksCN0fE6iUyc261whi1siV9AziJtM7YGqThKrf3MmMs5ebss4GPR8TVbeQ1cgfumGtmS1oG+D9SVdIJwKPAT4APRMGlN2rlNvJb/4ysmVsrexDP5zy394PAAcydu/V94FtR8Mtw7XPK2uEeqhZJ2hU4AbgJeA5wHfBc4A9AyT/aGwJvBP7E3N6DrUgLse4MfFfSbhFxTh9l1zzmi4CPk65ADTkQOL9AVtOFwLSIuFlpod0jgMeAT0bhhX0rZotUXXA26cvBiQWzxkIuwN7A2ZIuA+a5KBARw11I6IVBPOaa2T8n9SDsQ/qiPT0inpD0S+B1fZhb7TOy4mdzzexBPJ+3A86IiG80Pqc2A1alY8HdHqt2Tll7xtXegQHzeWDfiNgCmJ3/ewBQehjaOOBtEfGyiNgjIl5GKlX6eERsDbwX+FKfZdc85g8Bb5Q0DVha0t+BtwAfKZQ35LvAUFGIr5IumARwZOHcatkR8UFSifb3khrNl0n6s6Siv+taudkXcuaqwAaNn/VLhg7iMVfOfjHpb8m3gHtj7iLZL+vTXKj3GVkrt1r2gJ7PtT4ja55T1hIP+WuRpAeGKp5p7tpI44B/RcQqBXNnkkqId63Alm/fHxFL90t2xdzFgAdJ49E3Jc3dug34U+OPaBFD/74kjQfuztn/IQ3lWKlfszv2Y03gGNIctqLrx9TKlTQL2DAiSk+kfqr9GIhjrpWtVDlzTdKV7YnAbhExQ9KsEn+ra+cOZVf6jKySWzu7Yz/6/nyu9TlV85yy9riHql1357HSANPyhMj1SBP6S/oHqcJd03vydkhd0LP7LLtKbm7A3UAqIX5ZRJwWEX8s3ZjKHsj/vrYFrouIB/P2Cf2cLWmSpLfn4RM3kIYa7tWvucA/SWPwWzeIx1wzOyJmkYYo/xGYqrROUPGroLVyqfcZWSu3avYAns/VPqcqnlPWEs+hatdRwEuBKcDXSXNqniB1PZe0H/AzpYpzd5CulDwOvCk/vhGpzHg/Zdc85hOBX0j6JqkARnONjWLrUJEmvV4OLE4aXgDwElKFwdKqZEs6DXgtcAXp6t9eETG9ZGbN3OwE4ExJ/8f88w9KrnM2cMdcOVs5I4BPSLoK+A2wZMHMmrlQ7zOyVm617AE9n2t9RtY8p6wlHvJXkdL6BxMj4voWsiaQSpWuQVpz4dJoqWRpreyKuTcP81BExLMKZ29Imif2j8b9JSLimpK5tbJzg/nkiLi1VMZYys3ZVf59DeIx18yW9P8i4icd27YAdokClTpr5w6zL619Ro6F3DazB/h8rvE5NWbOKSvHDSozMzMzM7NR8hwqMzMzMzOzUXKDyszMzMzMbJTcoKpI0gGDlFsz28c8GNmDllsz28c8GNmDllsz28c8GNk1j9nKcYOqrlonVc2T2cfc/7k1swctt2a2j3kwsgctt2a2j3kwst2g6kNuUJmZmZmZmY2Sq/wtpGW1WKwyyjXhZvI4y45y7b6l115uVK8DmP7gw6w06Rmjfv3CqJW9sLmPLTv6RdRnzLiPFVZYflSvHf/gfaPOnT5rNistPXHUr18YC5P98MRVRp1734x7WX6FFUf9+lq5S8Qjo37tvffdz4rLj/7vwRyNbimURfmYa+VOnzP683H2zHuYuOzKo359rdxVdc+oXnfPzAdZedlJo85dGAub/Z8llxnV62bMmMEKK6ww6tyFsfDZGmXuvaywEH9HJsy+d9Svnf7AbFZaZvTn5L3jVx/V6x6ceQ+TRnlO3fvvaTw4c/roftlWlBf2XUirMIGvL7ZO67nbHvyG1jMH2YzX7Fsld+VLT6uSW9NVL/xg7V1o3fqP/7Va9k2LbVIlt+Yx13LsDVvV3oXWfWDJo2rvQuvu2nD72rvQusc0ugvLC2utS0+skgtwwqqfbD3zC+97YeuZNjIe8mdmZmZmZjZKblCZmZmZmZmNkhtUZmZmZmZmo+QGlZmZmZmZDQxJR0u6W9K1vXg/N6jMzMzMzGyQHAvs2Ks3c4PKzMzMzMwGRkRcCMzo1fu5bLqZmZmZmfXEC8ZNjAfi8ar7cBNzrgOaCx4eGRFHlspzg8rMzMzMzHrigXicb4xvf43Wptc/dsMjEbFlW3ke8mdmZmZmZjZKi1wPlaSzgT0i4v6O7ZcBSwArAM8A7sgP7QF8GVgPeBw4KyI+3t4em5mZmZkNCIEmqO4+PNZu3CLXQxURO3VpTAnYJiI2Bw4FTo2IzfP9acBXImJjYAvgJZJe2/Z+m5mZmZlZfZJOBi4FNpJ0u6R3Lcz7jekeKkk/B9YGlgS+GRFHSpoGbAlMAn4FnA9sA+wK3NL5HhHxUH4OEfEfSVcAa+X3PxZ4GNgYWAfYB9grv99lEbF3uaMzMzMzM7O2RcTuvXy/Md2gAvaNiBmSngFcLmlKx+MbAftExHtH8maSlgN2Br7Z2Lw8sD2wC3AW8BJgv5y3eURcubAHYWZmZmY2CCQxbnzlIX8tG+tD/g6UdBXwR1JP1QYdj98SEX8cyRtJGg+cDHwrIv7ZeOisiAjgGuDfEXFNRDwBXAdMHua9DpA0VdLUmdQtC2lmZmZmZvWM2R4qSdsBryLNjXpI0gWkoX9Ns5/GWx4J3BgR3+jYPif/94nG7aH7XX8/uY79kQAbaMl4GvtgZmZmZta/BJow1vtsemvMNqiAZYH7cmNqY2Dr0b6RpM/n99uvVztnZmZmZmY2lpuP5wDjJV0NfI407O9pk7QW8ClgE+AKSVdKcsPKzMzMzMwW2pjtoYqIOUC38uaT83+nA8/t8rpjgWMb928Hus6Ma1bxi4hpzfdzhT8zMzMzs6dJuCiFmZmZmZmZjYwbVGZmZmZmZqM0Zof8mZmZmZnZIkagCR7yZ2ZmZmZmZiPgHiozMzMzM+sJSS5KYWZmZmZmZiPjBpWZmZmZmdkoecifmZmZmZn1hotSmJmZmZmZ2Ui5h2ohLb32cmx78Btaz/39h85oPXOQvfLXm1XJfWLOnCq5Nd09e6nau9C6h5eo8+8LYNbsCVVyax5zLeMH8BP31o12rL0LrfvldWvX3oXWbbz2Y1Vyb9z/hCq5AA//5JOtZ0a0Hjk6wkUpzMzMzMzMbGTcoDIzMzMzMxulARyAYGZmZmZmJQjQYh7yZ2ZmZmZmZiPgHiozMzMzM+sNwTj3UJmZmZmZmdlIuEFlZmZmZmY2Sh7yZ2ZmZmZmPSI0zkP+zMzMzMzMbATcQ2VmZmZmZr0h0GKD1WczWEdrZmZmZmbWQ25QmZmZmZmZjdKYa1BJOkzSQR3bpkla6Slet7ekNRr3XynpCklXSvqDpPUbz7snb79S0n6N1+wl6cb8s1evj83MzMzMrJ+JtA5VzZ+29dMcqr2Ba4E78/3vAW+IiOslvRc4JD8H4NSIeH/zxZJWAD4DbAkE8GdJZ0bEfS3su5mZmZmZLYJ62qCS9GlgT+A2YDrwZ+B04DvAysBDwP4R8TdJ6wBH5+33APtExK2jzJhGagidKOlhYBtSo2iZ/LJlmdvQGs5rgPMiYkbOOQ/YETh5JMduZmZmZmaDp2cNKklbArsBW+T3vYLU2DkSeE9E3CjpRcB3ge2BbwPHR8RxkvYFvgXsmt/uw5Le3nj7NRaUERE/lfR+4KCImJqfux9wdm5gPQBs3Xi/3SS9HLgB+HBE3AasSWqkDbk9b+t2rAcABwCsvfwy3Z5iZmZmZjZ4hNehWggvBc6IiIcjYhZwFrAk8GLgNElXAj8AVs/P3wY4Kd8+Ib9+yNcjYvOhH+b2LnXLGM6HgZ0iYi3gGOBreftZwOSI2BT4DXBc3t7t/3x0e+OIODIitoyILVea9IwF7IKZmZmZmfWzXg7569YgGQfcnxtFT6Vr42UEGfM/SVoZ2CwiLsubTgXOAYiIextPPQo4It++Hdiu8dhawAUjyTMzMzMzM4A6hSFq6mUP1R+AnSUtKWkS8DrSnKmbJb0FQMlm+fmXAG/Lt/fMrx9NxpBZwNL59n3AspI2zPdfDVyf92H1xmt2GdoOnAvsIGl5ScsDO+RtZmZmZmZmXfWshyoiLpd0JnAVcAswFZhJaix9T9IhwATglPycA4GjJX2UXJRiITIAjgW+3yhKsT8wRdITpAbWvvl5B0raBXgMmEGu/BcRMyR9Drg8P+/woQIVZmZmZmZm3fS6bPpXIuIwSUsBFwJfjYibSdXy5hER00jFKTq3H9Zl2+QFZeTnTAGmNJ53ev7pfK9PAJ/otvMRcTSp8qCZmZmZmT1NEmjAhvz1ukF1pKRNSMUojouIK3r8/m1lmJmZmZmZPaWeNqgiYo9evl+tDDMzMzMzGx2N62WZhrFvsI7WzMzMzMysh9ygMjMzMzMzG6Vez6EyMzMzM7NBJdC4wSpK4R4qMzMzMzOzUXIPlZmZmZmZ9YgYN2Bl091DZWZmZmZmNkpuUJmZmZmZmY2Sh/yZmZmZmVlPyEUpzMzMzMzMbKTcQ2VmZmZmZj2jcYPVZzNYR2tmZmZmZtZDblCZmZmZmZmNkof8mZmZmZlZb7gohZmZmZmZmY2Ue6jMzMzMzKxHxLjF3ENlZmZmZmZmI+AGlZmZmZmZ2Sh5yJ+ZmZmZmfWEXJTCzMzMzMzMRmrMNagkHSbpoI5t0ySt9BSv21vSGo37r5R0haQrJf1B0vqN592Tt18pab/Ga/aSdGP+2avXx2ZmZmZmZv2ln4b87Q1cC9yZ738PeENEXC/pvcAh+TkAp0bE+5svlrQC8BlgSyCAP0s6MyLua2HfzczMzMz6gsaNuT6bonraoJL0aWBP4DZgOvBn4HTgO8DKwEPA/hHxN0nrAEfn7fcA+0TEraPMmEZqCJ0o6WFgG1KjaJn8smWZ29AazmuA8yJiRs45D9gROHkkx25mZmZmZoOnZw0qSVsCuwFb5Pe9gtTYORJ4T0TcKOlFwHeB7YFvA8dHxHGS9gW+Beya3+7Dkt7eePs1FpQRET+V9H7goIiYmp+7H3B2bmA9AGzdeL/dJL0cuAH4cETcBqxJaqQNuT1vMzMzMzOzkXBRioXyUuCMiHg4ImYBZwFLAi8GTpN0JfADYPX8/G2Ak/LtE/Lrh3w9IjYf+mFu71K3jOF8GNgpItYCjgG+lrefBUyOiE2B3wDH5e3d/s9HtzeWdICkqZKmTn/w4QXsgpmZmZmZ9bNeNqi6NUjGAfc3G0cR8exhXt+18TKCjPmfJK0MbBYRl+VNp5IadkTEvRExJ28/CnhBvn07sHbjbdZimGGCEXFkRGwZEVuuNOkZI9klMzMzMzPrQ71sUP0B2FnSkpImAa8jzZm6WdJbAJRslp9/CfC2fHvP/PrRZAyZBSydb98HLCtpw3z/1cD1eR9Wb7xml6HtwLnADpKWl7Q8sEPeZmZmZmZmIyI0ru5P23o2hyoiLpd0JnAVcAswFZhJaix9T9IhwATglPycA4GjJX2UXJRiITIAjgW+3yhKsT8wRdITpAbWvvl5B0raBXgMmEGu/BcRMyR9Drg8P+/woQIVZmZmZmZm3fS6bPpXIuIwSUsBFwJfjYibSdXy5hER00jFKTq3H9Zl2+QFZeTnTAGmNJ53ev7pfK9PAJ/otvMRcTSp8qCZmZmZmY3CoBWl6HWD6khJm5CKURwXEVf0+P3byjAzMzMzM3tKPW1QRcQevXy/WhlmZmZmZmYj0eseKjMzMzMzG1ASaFwv696NfYN1tGZmZmZmZj3kHiozMzMzM+uZcYsNVlEK91CZmZmZmZmNkhtUZmZmZmZmo+Qhf2ZmZmZm1hvSwK1D5R4qMzMzMzOzUXIPlZmZmZmZ9YzLppuZmZmZmdmIuEFlZmZmZmY2Sh7yZ2ZmZmZmPSHhohRmZmZmZmY2Mu6hMjMzMzOznnEPlZmZmZmZmY2IG1RmZmZmZmaj5CF/ZmZmZmbWI/I6VGZmZmZmZjYyblCZmZmZmZmNkof8mZmZmZlZb3gdKjMzMzMzMxupMdegknSYpIM6tk2TtNJTvG5vSWs07h8r6WZJV+afzfN2SfqWpJskXS3p+Y3X7CXpxvyzV6+PzczMzMysv6WiFDV/2tZPQ/72Bq4F7mxs+2hE/LTjea8FNsg/LwK+B7xI0grAZ4AtgQD+LOnMiLiv9I6bmZmZmdmiqacNKkmfBvYEbgOmA38GTge+A6wMPATsHxF/k7QOcHTefg+wT0TcOsqMaaSG0ImSHga2WcBbvAE4PiIC+KOk5SStDmwHnBcRM3LOecCOwMlP53dgZmZmZmaDo2d9YpK2BHYDtgDeRGrgABwJfCAiXgAcBHw3b/82qWGzKXAi8K3G2324MVTvSmCNBWXkXqipwJ4RsXlEPJzf5wt5WN/XJS2Rt61JaowNuT1vG257t2M9QNJUSVOnP/hwt6eYmZmZmQ0mqe5Py3o5yPClwBkR8XBEzALOApYEXgyclhtGPwBWz8/fBjgp3z4hv37I13PDaPOI2Jy5w/i6ZQznE8DGwAuBFYCP5e3dfsuxgO3zb4w4MiK2jIgtV5r0jAXsgpmZmZmZ9bNeDvnr1iAZB9yfG0VPpWvjZQQZ3d8s4q58c46kY0i9Y5B6ntZuPHUtUoPtdtKwv+b2C0aaZ2ZmZmY26OSy6QvlD8DOkpaUNAl4HWnO1M2S3gJPVtjbLD//EuBt+fae+fWjyRgyC1h66E6eF4UkAbuSClYAnAm8M+/L1sDM3Pg6F9hB0vKSlgd2yNvMzMzMzMy66lkPVURcLulM4CrgFtKcppmkxtL3JB0CTABOyc85EDha0kfJRSkWIgPgWOD7jaIUJ0pamdSrdSXwnvy8s4GdgJtIDb598nvPkPQ54PL8vMOHClSYmZmZmZl10+uy6V+JiMMkLQVcCHw1Im4mVcubR0RMA7bvsv2wLtsmLygjP2cKMKXxvPneOz8vgPcN89jRpMqDZmZmZmY2CjXWgqqp1w2qIyVtQipGcVxEXNHj928rw8zMzMzM7Cn1tEEVEXv08v1qZZiZmZmZ2ShILkphZmZmZmZmI+MGlZmZmZmZ2Sj1eg6VmZmZmZkNsEErSjFYR2tmZmZmZtZD7qEyMzMzM7OecVEKMzMzMzMzGxE3qMzMzMzMzEbJQ/4W0mPLrsSM1+zbeu4rf71Z65mD7Lc7HF4ld8uPbFUlt6ZnbTe99i60bv2bz66WfdO6O1XJrXnMtVzF3rV3oXXPvPbM2rvQut03fEntXWjdcv/4U53gKR+ukwtc90i0nvlE+5GjInnIn5mZmZmZmY2QG1RmZmZmZmaj5CF/ZmZmZmbWIwKvQ2VmZmZmZmYj4R4qMzMzMzPrGclFKczMzMzMzGwE3KAyMzMzMzMbJQ/5MzMzMzOz3hDIRSnMzMzMzMxsJNxDZWZmZmZmPSI0zkUpzMzMzMzMbARabVBJukDSlvn2NEkrNR7bTtIvepSzqqRfSLpK0l8lnT2K99hY0pWS/iJpvV7sl5mZmZmZ9ZeeD/lTKjyviHii1+/9NBwOnBcR38z7tOnTebGkxYBdgTMi4jMF9s/MzMzMrP8IcFGKp0/SZEnXS/oucAXwDkmXSrpC0mmSJj3N9ztM0kGN+9fmjMmS/ibph3nbiZJeJeliSTdK2iq/ZHXg9qHXR8TV+X3m6QWT9G1Je+fb0yQdKukPwFuBDwH7STp/dL8VMzMzMzPrd73sodoI2Ac4FPgZ8KqImC3pY8BHSL1Gnc6X9Hi+PQn42why1gfeAhwAXA7sAbwU2AX4JKln6TvAqZLeD/wGOCYi7hzBez8SES8FkLQh8GBEfGUErzMzMzMzMxi4ohS9bFDdEhF/lPR6YBPg4jT6j8WBS4d5zSsiYjqk3iPgoGGe13RzRFyTX3Md8NuICEnXAJMBIuJcSc8CdgReC/xF0nNH8N6njuA5SDqA1KBjzTXWGMlLzMzMzMysD/WyQTU7/1ek+Uu7L8R7Pca8wxGXbNye07j9ROP+EzSOJyJmACcBJ+Vhfi8H/r2A94W5x7BAEXEkcCTAps97bozkNWZmZmZm1n9KzBj7I/ASSesDSFoqD597OqYBz8+vfz6w7tN5saTtJS2Vby8NrAfcCtwCbCJpCUnLAq98mvtlZmZmZmbDEEIaV/WnbT2v8hcR9+RCDydLWiJvPgS44Wm8zRTgnZKuJM2TejqvBXgB8G1JQz1dP4yIywEk/QS4GrgR+MvTfF8zMzMzM7Mn9aRBFRHTgOc27v8OeGGX523XuD2547ELgAvy7YeBHYaJa+bs3W0fIuLLwJeH2deDgYO7bO/cn8OGyTczMzMzs24EDFhRisEqEm9mZmZmZtZDblCZmZmZmZmNUs/nUJmZmZmZ2eDSuMHqsxmsozUzMzMzM+sh91CZmZmZmVnPyEUpzMzMzMzMbCTcoDIzMzMzMxslD/kzMzMzM7PekECD1WczWEdrZmZmZmbWQw6vU0QAACAASURBVO6hWkjjH7yPlS89rfXcJ+bMaT1zkG35ka2q5E792p+q5Nb07Hc/VHsXWqeHZ1fLXmpcnd93zWOu5dbbBu/f9mPPX7v2LrTu/gkr196F1i07fkKV3MdvvqlKLsADSz3WeuYTj0frmTYyblCZmZmZmVnPuMqfmZmZmZmZjYh7qMzMzMzMrHfGDVafzWAdrZmZmZmZWQ+5QWVmZmZmZjZKHvJnZmZmZmY9IQnJRSnMzMzMzMxsBNxDZWZmZmZmveOiFGZmZmZmZjYSblCZmZmZmZmNkof8mZmZmZlZz2ici1KYmZmZmZnZCFRtUEmaJmmlxv3tJP2iR++9naSZkq7MP4c2HttR0t8l3STp443t60q6TNKNkk6VtHgv9sXMzMzMbCBIoHF1f1rW7z1UF0XE5vnncABJiwHfAV4LbALsLmmT/PwjgK9HxAbAfcC7auy0mZmZmZktGp6yQSVpoqRfSrpK0rWS3pp7lo6Q9Kf8s35+7sqSpki6PP+8JG+fJOkYSddIulrSbiPIPUzSQY3710qanH/+JumHeduJkl4l6eLcs7TVU7z1VsBNEfHPiPgP/H/27jzszqq+9//7w2BAQEBAK1ODgooGCBJRVJCpFq0H9BIEpZWobRxrf22p0/nZ6lF7pMMPtRY9WBXrBBpPFS2KMkTQyhAlkDDVVIKEUCsyB8Qmz/f3x14puw87yebJZu/k2e/Xdd3Xs+51r+mO/OH3+q57bc4GjkvnF8iOBOa3dp8FXra+dUqSJEkaX/0cSnEMsKKqfgcgyfZ0Mjn3VNXBSV4DfBh4KfAROhme7yfZEzgf2Bd4D3B3Ve3Xxtixa/yLk6xu5W2BG/pY097ACcA84Erg1cALgGOBd/NQIHRIkquBFcCpVXUtsBtwS9dYy4HnADsBd1XVqq763XpNnmRem5s9dtq+j+VKkiRJY2LMDqXoJ6BaDPxNktOAb1bVpZ1kDl9qz78EnN7KRwPPaM8BHpdku1Z/0prKqrqza/wjqup26Hz3BJzK+t1UVYtbn2uBC6uqkiwGZrY2PwZ+s6ruS/IS4GvAPkCv/4VrHfUPr6w6EzgT4Fl77dazjSRJkqTpb70BVVX9a5KDgJcA/zvJd9Y86m7W/m4GHFJVD3SP0bbTPdLAYxX/fUviVl3lB7vKE133E7R3qqp7ut7hvCRntAMwlgN7dPXfnU4G63ZghyRbtCzVmnpJkiRJfcoIDoYYpX6+odoVuL+qPg/8DfCs9ujErr8/bOXvAG/t6jt7LfXdW/7WZtmauZI8C9irjz7d6/6NFsjRvqvaDPglnS2C+7QT/R5DJ3N2blUVcDFwfBviFODrj2ROSZIkSeOln/BxP+CKJIuA/wl8oNXPSHI58EfAH7e6twFz2sET1wFvbPUfAHZsh0hcDRzRx7xfBR7f5n0T8K99vdFDjgfWzPdR4KTqWEUnuDsfuB74cvu2CuAdwJ8kWUrnm6pPPcI5JUmSJI2Rfrb8nU8n+PgvLfHz91X1vkltb+ehzFV3/X10Mj6T62dOul8ALGjlB4AXrWVZs7r6zO0qL1vzrKo+BnxsLe90HnBej/qf0jkFUJIkSdIjFcbuUIrx2uAoSZIkSQPUzyl/DzM5syRJkiRJELLZeOVsxuttJUmSJGmADKgkSZIkaYqmtOVPkiRJknqKh1JIkiRJkvpgQCVJkiRJU+SWP0mSJEmDEcBT/iRJkiRJ/TBDJUmSJGlA4qEUkiRJkqT+GFBJkiRJ0hS55U+SJEnSwMRDKSRJkiRJ/TBDJUmSJGkwAmS8cjbj9baSJEmSNEAGVJIkSZI0RW75kyRJkjQggc38HSpJkiRJUh/MUEmSJEkaiADxUApJkiRJUj9GGlAlWZZk5677w5N8c0Bjb5/kG0muTnJtktd2PVudZFG7zu2q3yvJ5Ul+kuScJI8ZxFokSZIkTU/TOUP1FuC6qjoAOBz4264A6YGqmt2uY7v6nAacXlX7AHcCrx/qiiVJkqRNWegcSjHKa8jWG1Al2SbJP7dMz5IkJ7bM0mlJrmjX3q3tLkm+muTKdj2/1W+b5DNJFie5Jskr+pj3vUlO7bpfkmRmu25I8g+t7gtJjk7yg5ZZOrh1KWC7JAG2Be4AVq1jvgBHAvNb1WeBl61vnZIkSZLGVz+HUhwDrKiq34HOVjo6mZx7qurgJK8BPgy8FPgInQzP95PsCZwP7Au8B7i7qvZrY+zYNf7FSVa38rbADX2saW/gBGAecCXwauAFwLHAu+kEQh8DzgVWANsBJ1bVROu/VZKFdAKsD1XV14CdgLuqak3QtRzYrY+1SJIkSQIgMGaHUvQTUC0G/ibJacA3q+rSTjKHL7XnXwJOb+WjgWe05wCPS7Jdqz9pTWVV3dk1/hFVdTt0vqECTmX9bqqqxa3PtcCFVVVJFgMzW5vfBhbRyTo9Bfhukkur6h5gz6pakeTJwEWt3z095qlekyeZRyeYY4+dtu9juZIkSZKmo/WGj1X1r8BBdAKr/53kz9c86m7WNd4hXd8n7VZV99LZTdkzOFmHVZPWt1VX+cGu8kTX/QQPBYmvBf5vdSwFbgKe3t5pRfv7U2ABcCBwO7BDkjX9d6eT3XqYqjqzquZU1Zydt9vmEb6WJEmSpOmin2+odgXur6rPA38DPKs9OrHr7w9b+TvAW7v6zl5LffeWv7VZtmauJM8C9uqjT7efAUe1/k8Engb8NMmOSWa0+p2B59M5vKKAi4HjW/9TgK8/wjklSZKk8ZaM9hqyfjY47gdckWQR8D+BD7T6GUkuB/4I+ONW9zZgTjt44jrgja3+A8CO7RCJq4Ej+pj3q8Dj27xvAv61rzd6yPuB57XtfBcC72hbC/cFFrZ1XEznG6rrWp93AH+SZCmdb6o+9QjnlCRJkjRG1vsNVVWdT+dwif/SvpH6+6p636S2t/NQ5qq7/j46GZ/J9TMn3S+gswWPqnoAeNFaljWrq8/crvKyNc/atr6H9a+qf6ETJD5M2wJ4cK9nkiRJkvqw2XgdSjFebytJkiRJA9TPKX8PMzmzJEmSJEnjaEoBlSRJkiQ9TMbvd6jG620lSZIkaYAMqCRJkiRpitzyJ0mSJGlwNhv+b0GNkhkqSZIkSZoiM1SSJEmSBsdDKSRJkiRJ/TBDtYEe2OYJXP3sPxr6vP+x8rFDn3OcPfnw20cy775vuH8k847S9U978aiXMHS3Xv3jkc19z9OeNZJ5R/nOo/KX93101EsYurt3fNGolzB0u9513aiXMHQL9zhxJPPe84QZI5kX4M8ufefQ5/zm6luHPqf6Y0AlSZIkaXDioRSSJEmSpD6YoZIkSZI0GAlsNl45m/F6W0mSJEkaIAMqSZIkSZoit/xJkiRJGhwPpZAkSZIk9cMMlSRJkqTByXjlbMbrbSVJkiRpgAyoJEmSJGmK3PInSZIkaTD8HapNU5Idkrx5Ut1fJbk2yfVJPpp0jhtJsiDJjUkWtesJrX5GknOSLE1yeZKZw38TSZIkSZuSTT6gSrI5sAPw5q665wHPB/YHZgHPBl7Y1e3kqprdrv9oda8H7qyqvYHTgdOGsX5JkiRpWklGew3ZRrPlL8lpwM1VdUa7fy9wL/AbwIuBAj5QVeckORz4C+A2YDZwDfCUJIuA7wL/F9gKeAwQYEvg5+tZwnHAe1t5PvCxJKmqGtArSpIkSZpmNpqACjgb+DBwRrt/JZ0s0THAAcDOwJVJLmnPDwZmVdVNbXverKqavWawJBfTCbgCfKyqru+a6zNJVgNfpROkFbAbcAtAVa1KcjewE3D7o/CukiRJkqaBjWbLX1VdBTwhya5JDgDupJN9+lJVra6qnwPfo7N9D+CKqrqp11hJ9gb2BXanEygdmeSw9vjkqtoPOLRdv7emW69lrWX8eUkWJll45x2/fMTvKkmSJE1b2Wy015BtNAFVMx84HjiRTsZqXZsgV67j2cuBy6rqvqq6D/gW8FyAqrq1/b0X+CKdTBfAcmAPgCRbANsDd/QavKrOrKo5VTVnx8fv1OerSZIkSZpuNraA6mzgJDpB1XzgEuDEJJsn2QU4DLiiR797ge267n8GvDDJFkm2pHMgxfXtfmeAVv9SYEnrcy5wSisfD1zk91OSJEmS1mVj+oaKqro2yXbArVV1W5J/Ag4Brqaz/e7tVfXvSZ4+qd8vk/wgyRI62ah3AkcCi1u/b1fVN5JsA5zfgqnNgQuAT7ZhPgV8LslSOpmpkx71F5YkSZKmldGctDdKG1VABdC+b1pTLuDP2tXdZgGwYFLdqycN9YYeY68EDlrLvL8CTpjKmiVJkiSNp40uoJIkSZK0iQqw2cb2VdGja7zeVpIkSZIGyIBKkiRJkqbILX+SJEmSBqKAGrNDKcxQSZIkSdIUmaGSJEmSNCCBjFfOZrzeVpIkSZIGyIBKkiRJkqbILX+SJEmSBsctf5IkSZKkfpihkiRJkjQwHpsuSZIkSeqLAZUkSZIkTZFb/iRJkiQNRvwdKkmSJElSn8xQSZIkSRocD6WQJEmSJPXDgEqSJEmSpsgtf5IkSZIGZ7PxytmM19tKkiRJ0gCZoZIkSZI0IKE8lEKSJEmS1I9pEVAl2SHJmyfVrU6yqF3ndtXvleTyJD9Jck6Sx7T6Ge1+aXs+c7hvIUmSJGlTs8kHVEk2B3YA3jzp0QNVNbtdx3bVnwacXlX7AHcCr2/1rwfurKq9gdNbO0mSJEn9CpDNRnsN2UYTUCU5rTvLlOS9Sf40yV8nWZJkcZIT27PDk1yc5IvAYuBDwFNaNuqv1zFHgCOB+a3qs8DLWvm4dk97flRrL0mSJEk9bUyHUpwNfBg4o92/kk6W6BjgAGBn4Mokl7TnBwOzquqmtj1vVlXN7hpvqyQLgVXAh6rqa8BOwF1Vtaq1WQ7s1sq7AbcAVNWqJHe39rdPXmiSecA8gCc+aY8NfG1JkiRJm6qNJqCqqquSPCHJrsAudLbjzQa+VFWrgZ8n+R7wbOAe4IqqumkdQ+5ZVSuSPBm4KMni1u9hU7e/vbJR1aOOqjoTOBNg31kH9mwjSZIkjaMawba7UdrY3nY+cDxwIp2M1bq23K1c10BVtaL9/SmwADiQTrZphyRrAsndgRWtvBzYA6A93x64YyovIUmSJGk8bGwB1dnASXSCqvnAJcCJSTZPsgtwGHBFj373AtutuUmyY5IZrbwz8Hzguqoq4OI2PsApwNdb+dx2T3t+UWsvSZIkqS+BjPgaso1myx9AVV2bZDvg1qq6Lck/AYcAV9PZfvf2qvr3JE+f1O+XSX6QZAnwLeCfgP+TZIJO0PihqrquNX8HcHaSDwBXAZ9q9Z8CPpdkKZ3M1EmP7ttKkiRJ2tRtVAEVQFXt11Uu4M/a1d1mAZ1tfN11r5401H700LYAHtyj/lfACVNZsyRJkqTxtNEFVJIkSZI2XR5KIUmSJEnqixkqSZIkSYMzgoMhRskMlSRJkiRNkQGVJEmSJE2RW/4kSZIkDUYCHkohSZIkSeqHGSpJkiRJA1FAeSiFJEmSJKkfZqg20Iz6FXuvvm7o8z4w44ChzznO9r7pvJHMmwdWjmTeUbr16h+PeglDd88BzxrZ3I8b0b/3KN95VOZfeMOolzB0L1/w7lEvYege/I/bR72EoVt99ttHMu+Mcy4fybwAFz3vtKHPee8nFgx9TvXHgEqSJEnS4HgohSRJkiSpH2aoJEmSJA1M4aEUkiRJkqQ+GFBJkiRJ0hS55U+SJEnSgITyUApJkiRJUj/MUEmSJEkaHDNUkiRJkqR+GFBJkiRJ0hS55U+SJEnSYAQq/g7VRi/Jrknmr+XZgiRzWvlVSRYnuSbJt5Ps3Orfm+TWJIva9ZKu/u9KsjTJjUl+ezhvJEmSJGlTtMkFVEm2qKoVVXX8+toBHwGOqKr9gWuAt3Y1Ob2qZrfrvNbnGcBJwDOBY4Azkmz+qLyIJEmSpE3e0AKqJDOT3JDkH5IsSfKFJEcn+UGSnyQ5OMk2ST6d5MokVyU5rvWdm+QrSb4BfKeNtaQ92zrJ2S0LdQ6w9Zop27VNkgCPA1asZ5nHAWdX1YNVdROwFDj4UfjnkCRJkqadar9DNcpr2Ib9DdXewAnAPOBK4NXAC4BjgXcD1wEXVdXrkuwAXJHkgtb3EGD/qrojycyuMd8E3F9V+yfZH/gxQFX9Z5I3AYuBlcBPgLd09XtrktcAC4E/rao7gd2Ay7raLG91kiRJkvQwww7hbqqqxVU1AVwLXFhVRSfomQm8CHhnkkXAAmArYM/W97tVdUePMQ8DPg9QVdfQ2dpHki3pBFsHAru2+ne1Ph8HngLMBm4D/rbV9/qCriZXJJmXZGGShb+8866+X16SJEma9pLRXkM27IDqwa7yRNf9BJ1sWYBXdH3btGdVXd/arFzHuA8LeugES1TVv7Wg7cvA81rdz6tqdQvsPslD2/qWA3t0jbE7PbYJVtWZVTWnqubstOMO61iWJEmSpOlsYzuU4nzgD9s3TyQ5sI8+lwAnt/azgP1b/a3AM5Ls0u5/C7i+tXtSV/+XA0ta+VzgpCQzkuwF7ANcMfXXkSRJkjSdbWy/Q/V+4MPANS2oWga8dD19Pg58Jsk1wCJaAFRVK5K8D7gkyX8CNwNzW5+/SjKbTmZrGfCG1ufaJF+m8y3XKuAtVbV6YG8nSZIkTXOjOBhilIYWUFXVMmBW1/3ctTx7Q4++ZwFn9WpfVQ/QOeq815yfAD7Ro/731rHODwIfXNtzSZIkSVpjY8tQSZIkSdpkhep5ztv0NV75OEmSJEkaIAMqSZIkSZoit/xJkiRJGphxO5RivN5WkiRJkgbIDJUkSZKkwQgQD6WQJEmSJPXBgEqSJEmSpsgtf5IkSZIGJNSY5WwMqDbQg9mKpZs/Y+jz3rtyy6HPOc6W7vWSkcz72M3uH8m8o3TP05416iUM3eOu/vHI5r7ngNH8e4/ynUflqEUfGvUShu7uw48f9RKG7rErbx/1EoZu8+NPH8m8Kx/cfCTzAhz5L+8Y+pzbrbx16HOqPwZUkiRJkgaigPJQCkmSJElSPwyoJEmSJGmK3PInSZIkaWAq45WzGa+3lSRJkqQBMqCSJEmSpClyy58kSZKkgSk85U+SJEmS1AczVJIkSZIGJB5KIUmSJEnqjwGVJEmSJE2RW/4kSZIkDUzFQyk2OUl2SPLmSXWrkyxq17ld9XsluTzJT5Kck+QxrX5Gu1/ans8c7ltIkiRJ2tRs8gFVks2BHYA3T3r0QFXNbtexXfWnAadX1T7AncDrW/3rgTuram/g9NZOkiRJUp+KzrHpo7yGbaMJqJKc1p1lSvLeJH+a5K+TLEmyOMmJ7dnhSS5O8kVgMfAh4CktG/XX65gjwJHA/Fb1WeBlrXxcu6c9P6q1lyRJkqSeNqZvqM4GPgyc0e5fSSdLdAxwALAzcGWSS9rzg4FZVXVT2543q6pmd423VZKFwCrgQ1X1NWAn4K6qWtXaLAd2a+XdgFsAqmpVkrtb+9snLzTJPGAewBOftMcGvrYkSZKkTdVGE1BV1VVJnpBkV2AXOtvxZgNfqqrVwM+TfA94NnAPcEVV3bSOIfesqhVJngxclGRx6/ewqdvfXtmo6lFHVZ0JnAmw76wDe7aRJEmSxk78HapRmw8cD5xIJ2O1ri13K9c1UFWtaH9/CiwADqSTbdohyZpAcndgRSsvB/YAaM+3B+6YyktIkiRJGg8bW0B1NnASnaBqPnAJcGKSzZPsAhwGXNGj373AdmtukuyYZEYr7ww8H7iuqgq4uI0PcArw9VY+t93Tnl/U2kuSJEnq07gdSrHRbPkDqKprk2wH3FpVtyX5J+AQ4Go62+/eXlX/nuTpk/r9MskPkiwBvgX8E/B/kkzQCRo/VFXXtebvAM5O8gHgKuBTrf5TwOeSLKWTmTrp0X1bSZIkSZu6jSqgAqiq/brKBfxZu7rbLKCzja+77tWThtqPHtoWwIN71P8KOGEqa5YkSZI0nja6gEqSJEnSpstDKSRJkiRJfTFDJUmSJGlgRnEwxCiZoZIkSZKkKTKgkiRJkqQpcsufJEmSpIEo4qEUkiRJkjRdJTkmyY1JliZ554aOZ4ZKkiRJ0sBszIdSJNkc+Hvgt4DlwJVJzq2q66Y6phkqSZIkSePiYGBpVf20qn4NnA0ctyEDmqHaQDPqV+y9esoB7ZQ9MOOAoc85zva+6byRzJsHVo5k3lG69eofj3oJQ3fPAc8a2dyPG9G/9yjfeVTmX3jDqJcwdC9f8O5RL2HoHvyP20e9hKFbffbbRzLvjHMuH8m8ABc977Shz3nvJxYMfc5N2M5JFnbdn1lVZ7bybsAtXc+WA8/ZkMkMqCRJkiQNTGXkW/5ur6o5a3nWa3G1IZO55U+SJEnSuFgO7NF1vzuwYkMGNKCSJEmSNC6uBPZJsleSxwAnAeduyIBu+ZMkSZI0MFUj3/K3VlW1KslbgfOBzYFPV9W1GzKmAZUkSZKksVFV5wEDO3HMgEqSJEnSgIQas6+KxuttJUmSJGmADKgkSZIkaYrc8idJkiRpIAqonj/1NH2ZoZIkSZKkKTJDJUmSJGlgzFBtApLsmmT+Wp4tSDKnlV+VZHGSa5J8O8nOrf69SW5NsqhdL+nq/64kS5PcmOS3h/NGkiRJkjZFm1xAlWSLqlpRVcevrx3wEeCIqtofuAZ4a1eT06tqdrvOa32eQefXkp8JHAOckWTzR+VFJEmSJG3yhhZQJZmZ5IYk/5BkSZIvJDk6yQ+S/CTJwUm2SfLpJFcmuSrJca3v3CRfSfIN4DttrCXt2dZJzm5ZqHOArddM2a5tkgR4HLBiPcs8Dji7qh6sqpuApcDBj8I/hyRJkjQtFRnpNWzD/oZqb+AEYB5wJfBq4AXAscC7geuAi6rqdUl2AK5IckHrewiwf1XdkWRm15hvAu6vqv2T7A/8GKCq/jPJm4DFwErgJ8Bbuvq9NclrgIXAn1bVncBuwGVdbZa3OkmSJEl6mGFv+bupqhZX1QRwLXBhVRWdoGcm8CLgnUkWAQuArYA9W9/vVtUdPcY8DPg8QFVdQ2drH0m2pBNsHQjs2urf1fp8HHgKMBu4DfjbVt8rpK3JFUnmJVmYZOEv77yr75eXJEmSprfRZqdGkaEadkD1YFd5out+gk62LMArur5t2rOqrm9tVq5j3IcFPXSCJarq31rQ9mXgea3u51W1ugV2n+ShbX3LgT26xtidHtsEq+rMqppTVXN22nGHdSxLkiRJ0nS2sR1KcT7wh+2bJ5Ic2EefS4CTW/tZwP6t/lbgGUl2afe/BVzf2j2pq//LgSWtfC5wUpIZSfYC9gGumPrrSJIkSZrONrbfoXo/8GHgmhZULQNeup4+Hwc+k+QaYBEtAKqqFUneB1yS5D+Bm4G5rc9fJZlNJ7O1DHhD63Ntki/T+ZZrFfCWqlo9sLeTJEmSprmq8fodqqEFVFW1DJjVdT93Lc/e0KPvWcBZvdpX1QN0jjrvNecngE/0qP+9dazzg8AH1/ZckiRJktbY2DJUkiRJkjZRBSM5GGKUNrZvqCRJkiRpk2FAJUmSJElT5JY/SZIkSQPjlj9JkiRJUl/MUEmSJEkaGDNUkiRJkqS+GFBJkiRJ0hS55U+SJEnSgIQqt/xJkiRJkvpgQCVJkiRJU+SWP0mSJEkDUcCEp/xJkiRJkvphhkqSJEnSwPg7VJIkSZKkvhhQSZIkSdIUueVPkiRJ0mAU/g6VJEmSJKk/ZqgkSZIkDYyHUkiSJEmS+mJAJUmSJElT5JY/SZIkSQMSD6XYFCTZNcn8tTxbkGROV/nGJIva9YRWPyPJOUmWJrk8ycyu/u9q9Tcm+e1hvI8kSZKkTdMml6FKskVVrQCO77PLyVW1cFLd64E7q2rvJCcBpwEnJnkGcBLwTGBX4IIkT62q1YNavyRJkjRdFR5K8ahJMjPJDUn+IcmSJF9IcnSSHyT5SZKDk2yT5NNJrkxyVZLjWt+5Sb6S5BvAd9pYS9qzrZOcneSaJOcAW/exnOOAz7byfOCoJGn1Z1fVg1V1E7AUOHjQ/xaSJEmSpodhZ6j2Bk4A5gFXAq8GXgAcC7wbuA64qKpel2QH4IokF7S+hwD7V9Ud3Vv0gDcB91fV/kn2B348ac7PJFkNfBX4QFUVsBtwC0BVrUpyN7BTq7+sq+/yVvffJJnX3oHdn/QbU/l3kCRJkjQNDDuguqmqFgMkuRa4sKoqyWJgJrA7cGySU1v7rYA9W/m7VXVHjzEPAz4KUFXXJLmm69nJVXVrku3oBFS/B/wj9MxD1jrq/3tF1ZnAmQCzn7nvw55LkiRJ48pDKR5dD3aVJ7ruJ+gEdwFeUVWz27VnVV3f2qxcx7g9g5qqurX9vRf4Ig9t31sO7AGdb7KA7YE7uuub3YEV/b2aJEmSpHGzsZ3ydz7wh+17JpIc2EefS4CTW/tZwP6tvEWSnVt5S+ClwJLW51zglFY+ns42w2r1J7VTAPcC9gGuGMSLSZIkSeNgYsTXsG1sp/y9H/gwcE0LqpbRCYTW5eN0vpO6BljEQwHQDOD8FkxtDlwAfLI9+xTwuSRL6WSmTgKoqmuTfJnOt1yrgLd4wp8kSZKktRlaQFVVy4BZXfdz1/LsDT36ngWc1at9VT1AC4h6OGgta/kVncMxej37IPDBtYwnSZIkSf9lY8tQSZIkSdqEeSiFJEmSJKkvBlSSJEmSNEVu+ZMkSZI0EEWonj/tOn2ZoZIkSZKkKTJDJUmSJGlgPJRCkiRJktQXAypJkiRJmiK3/EmSJEkamHE7lMKAagPd/uA2nPWvBw993i38X26ormbuSOb92S33j2TeUfrL+z466iUM3fwLbxjZ3Ect+tBI5h3lO4/Kjkc9fdRLGLrvL1w86iUM3X5bj99/2/f+zjNGMu/Plm8+knkBnv2Xhwx9zi1vuXHoc6o//t9ySZIkSYNRMFGjXsRw+Q2VJEmSJE2RAZUkJKFnVgAAIABJREFUSZIkTZFb/iRJkiQNRDF+h1KYoZIkSZKkKTJDJUmSJGlgqsxQSZIkSZL6YEAlSZIkSVPklj9JkiRJA1P+DpUkSZIkqR9mqCRJkiQNSJjw2PRNV5K5ST7WynsmuTjJVUmuSfKSVj8zyQNJFrXrE139D0qyOMnSJB9NMl7/NUiSJEl6RKZzhur/Bb5cVR9P8gzgPGBme/ZvVTW7R5+PA/OAy1r7Y4BvDWGtkiRJkjZBQ89QJXljV3boppZFelXLDC1JclpX2/uSfDDJ1UkuS/LEVv8/klzesk8XrKmfpIDHtfL2wIr1rOtJwOOq6odVVcA/Ai8byEtLkiRJY6Do/A7VKK9hG3pAVVWfaNmhZwPLgbOA04AjgdnAs5OsCWS2AS6rqgOAS4A/aPXfB55bVQcCZwNv7zHVe4HfTbKcTrbpD7ue7dWCse8lObTV7dbWs8byVidJkiRJPY1yy99HgIuAu4AFVfULgCRfAA4Dvgb8Gvhma/8j4LdaeXfgnJZVegxwU4/xXwWcVVV/m+QQ4HNJZgG3AXtW1S+THAR8LckzoefXcz0PfUwyj87WQHbYec9H9taSJEnSNOax6UOQZC7wm8D76B3IrPGfbfsdwGoeCgD/DvhYVe0HvAHYqkff1wNfBqiqH7Y2O1fVg1X1y1b/I+DfgKfSyUjt3tV/d9ayTbCqzqyqOVU1Z5vtd1nP20qSJEmarkbxDdVBwKnA71bVBHA58MIkOyfZnE5m6XvrGWZ74NZWPmUtbX4GHNXm3JdOQPWLJLu0eUjyZGAf4KdVdRtwb5LnttP9XgN8farvKUmSJGn6G8WWv7cCjwcubqeSLwTeBVxMJ1t1XlWtL5B5L/CVJLfSOZFvrx5t/hT4ZJI/prN1b25VVZLDgP+VZBWdrNcbq+qO1udNdL7p2prO6X6e8CdJkiQ9AjVmv0M19ICqql67lkdf7NF2267yfGB+K3+dHtmjqjqLTkBEVV0HPL9Hm68CX13L2hYCs9bzCpIkSZIETLMf9pUkSZKkYZrOP+wrSZIkaZgKJjzlT5IkSZLUDzNUkiRJkgaigKrxOpTCDJUkSZIkTZEBlSRJkiRNkVv+JEmSJA1MeSiFJEmSJKkfZqgkSZIkDcwEHkohSZIkSeqDAZUkSZIkTZFb/iRJkiQNjIdSSJIkSZL6YoZKkiRJ0kAUocpDKSRJkiRJfTCgkiRJkqQpcsufJEmSpMEomPBQCkmSJElSP8xQSZIkSRoYj02XJEmSJPXFgEqSJEmSpsgtf5IkSZIGpvB3qDZZSeYm+VhX+RdJFrXr97vanZLkJ+06pav+oCSLkyxN8tEk4/VfgyRJkqRHZLpnqM6pqrd2VyR5PPAXwByggB8lObeq7gQ+DswDLgPOA44BvjXcJUuSJEmbpsJj0x91Sd7YlTW6KcnFSV7VMkNLkpzW1fa+JB9McnWSy5I8sdX/jySXJ7kqyQVr6vv028B3q+qOFkR9FzgmyZOAx1XVD6uqgH8EXjbAV5ckSZI0zQw9oKqqT1TVbODZwHLgLOA04EhgNvDsJGsCmW2Ay6rqAOAS4A9a/feB51bVgcDZwNvXMt0rklyTZH6SPVrdbsAtXW2Wt7rdWnly/cMkmZdkYZKFK+/+RZ9vLkmSJGm6GeU3VB8BLgLuAhZU1S+qahXwBeCw1ubXwDdb+UfAzFbeHTg/yWLgz4Bn9hj/G8DMqtofuAD4bKvv9V1UraP+4ZVVZ1bVnKqas832u6z9DSVJkqQxUzXaa9hGElAlmQv8JvA+egcya/xn234HsJqHvvn6O+BjVbUf8AZgq8kdq+qXVfVgu/0kcFArLwf26Gq6O7Ci1e/eo16SJEmSehrFN1QHAacCv1tVE8DlwAuT7Jxkc+BVwPfWM8z2wK2tfEqvBu2bqDWOBa5v5fOBFyXZMcmOwIuA86vqNuDeJM9tp/u9Bvj6I39DSZIkSeNiFKf8vRV4PHBxO5V8IfAu4GI62arzqmp9gcx7ga8kuZXOiXx79WjztiTHAquAO4C5AFV1R5L3A1e2dv+rqu5o5TfR+aZrazqn+3nCnyRJkvQIjGLb3SgNPaCqqteu5dEXe7Tdtqs8H5jfyl+nR/aoqs6iExBRVe+iE6j1WsOngU/3qF8IzFrPK0iSJEkSMP1/h0qSJEnSkFTBRK3riITpZ5Sn/EmSJEnSJs2ASpIkSZKmyC1/kiRJkgZm3A6lMEMlSZIkSVNkhkqSJEnSwJihkiRJkiT1xYBKkiRJkqbILX+SJEmSBmZizLb8GVBtoCfmF/zhVp8c+rw/e9oxQ59znO255NyRzLvqWXuMZN5RunvHF416CUP38gXvHtncdx9+/EjmHeU7j8r3Fy4e9RKGbtWc/Ua9hOF7y+xRr2Dojj70OSOZ9+bZJ4xkXoA7P/md4U96wguHP6f6YkAlSZIkaSAKqMqolzFUfkMlSZIkSVNkQCVJkiRJU+SWP0mSJEmDUf4OlSRJkiSpT2aoJEmSJA3MuB2bboZKkiRJkqbIgEqSJEmSpsgtf5IkSZIGovM7VKNexXCZoZIkSZKkKTKgkiRJkqQpGkpAleS8JDus5dmyJIuTLGrX85LMTLJkAPP+VZJrk1yf5KNJ0uoXJLmxa84ntPoZSc5JsjTJ5UlmbugaJEmSpHFSNdpr2IbyDVVVvWRyXQtu0m6PqKrbu57N3NA5kzwPeD6wf6v6PvBCYEG7P7mqFk7q9nrgzqraO8lJwGnAiRu6FkmSJEnT08AzVEm+luRHLTM0r9UtS7Jzyzxdn+QM4MfAHn2Mt1WSz7Qs1lVJjmj15yXZv5WvSvLnrfz+JL9P55u4rYDHADOALYGfr2e644DPtvJ84Kg1WS1JkiRJ6zdRo72G7dHY8ve6qjoImAO8LclOk54/DfjHqjqwqm5udRe3rXeX9xjvLQBVtR/wKuCzSbYCLgEOTfI4YBWdbBTAC4BLq+qHwMXAbe06v6qu7xr3M23O93QFTbsBt7T5VgF3A5PXT5J5SRYmWfiLu+/r+x9GkiRJ0vTyaARUb0tyNXAZnQzUPpOe31xVl02qO6KqZlfVc3qM9wLgcwBVdQNwM/BU4FLgsPb8n4FtkzwWmFlVNybZG9gX2J1OoHRkksPamCe3AO3Qdv1eq++VjXpYnFtVZ1bVnKqas8v22671H0KSJEnS9DbQgCrJ4cDRwCFVdQBwFZ1td91WPtJh11J/JZ0s2KF0slVXAX8A/Kg9fzlwWVXdV1X3Ad8CngtQVbe2v/cCXwQObn2W07YhJtkC2B644xGuV5IkSRpPIz6QYhSHUgw6Q7U9nUMd7k/ydFoAs4EuAU4GSPJUYE/gxqr6NZ3tea+kkw27FDi1/QX4GfDCJFsk2ZLOgRTXt/ud23hbAi8F1pwoeC5wSisfD1xUNW4/TSZJkiSpX4M+5e/bwBuTXAPcSCfQ2VBnAJ9IspjOt1Jzq+rB9uxS4KgWwF1KZ3vfmoBqPnAksJjOtr1vV9U3kmwDnN+Cqc2BC4BPtj6fAj6XZCmdzNRJA1i/JEmSNBYKmJgY9SqGa6ABVQt0Xtzj0cz293Zg1qQ+Myc3rqpla9pV1a+AuWuZ7z3Ae1p5BV3bA6tqNfCGHn1WAgetZbxfASf0eiZJkiRJkw3lh30lSZIkaToayg/7SpIkSRoP43YCgRkqSZIkSZoiM1SSJEmSBsYMlSRJkiSpLwZUkiRJkjRFbvmTJEmSNBBVMOGWP0mSJElSP8xQSZIkSRqYGrNTKcxQSZIkSdIUGVBJkiRJ0hS55U+SJEnSwIzZjj8zVJIkSZI0VWaoJEmSJA3MxMSoVzBcZqgkSZIkaYoMqCRJkiRpitzyJ0mSJGkgqjyUQpIkSZLUJwMqSZIkSZoit/xJkiRJGpgJt/xJkiRJkvoxlIAqyXlJdljLs2VJFidZ1K7nJZmZZMkA5l3dNe65XfV7Jbk8yU+SnJPkMa1+Rrtf2p7P3NA1SJIkSeNkzcEUo7qGbSgBVVW9pKru6q5Lx5r5j6iq2e36lwFO/UDXuMd21Z8GnF5V+wB3Aq9v9a8H7qyqvYHTWztJkiRJ6mngAVWSryX5UZJrk8xrdcuS7NwyT9cnOQP4MbBHH+NtleQzLYt1VZIjWv15SfZv5auS/Hkrvz/J769jvABHAvNb1WeBl7Xyce2e9vyo1l6SJEmSHubROJTidVV1R5KtgSuTfHXS86cBr62qNwO0eOXiJKuBB6vqOZPavwWgqvZL8nTgO0meClwCHJpkGbAKeH5r/wLg8628VZKF7fmHquprwE7AXVW1qrVZDuzWyrsBt7T5ViW5u7W/fer/HJIkSdL4qDE7leLRCKjeluTlrbwHsM+k5zdX1WWT6o6oqrUFLS8A/g6gqm5IcjPwVOBS4G3ATcA/A7+V5LHAzKq6sfXds6pWJHkycFGSxcA9PeZY8796r2zUw/6LaJm3eQB77PL4tSxbkiRJ0nQ30C1/SQ4HjgYOqaoDgKuArSY1W/lIh11L/ZXAHOBQOtmqq4A/AH60pkFVrWh/fwosAA6kk23aIcmaYHJ3YEUrL6dtQ2zPtwfumDxxVZ1ZVXOqas4u22/7CF9HkiRJmp6qOsemj/IatkF/Q7U9nUMd7m/b8547gDEvAU4GaFv99gRurKpf09me90rgMjoZq1PbX5LsmGRGK+9MZ0vgdVVVwMXA8W38U4Cvt/K57Z72/KLWXpIkSZIeZtAB1beBLZJcA7yfTqCzoc4ANm/b9c4B5lbVg+3ZpcDPq+r+Vt69/QXYF1iY5Go6AdSHquq69uwdwJ8kWUrnG6lPtfpPATu1+j8B3jmA9UuSJEmapgb6DVULdF7c49HM9vd2YNakPjMnN66qZWvaVdWvgLlrme89wHtaeQVd2wPb8ev7raXfT4GDe9T/CjihVx9JkiRJ6zdu+7uG8jtUkiRJkjQdPRqn/EmSJEkaUxNjdmy6GSpJkiRJmiIDKkmSJEmaIrf8SZIkSRqIwkMpJEmSJEl9MkMlSZIkaTDKDJUkSZIkqU8GVJIkSZI0RW75kyRJkjQgxcSY7fkzoNpAv97qcdz21COHPu8/X7vH0OccZ6966vNHMu9dW+4yknlHade7rhv1Eobuwf+4fWRzP3blaOYe5TuPyn5b3zDqJQzfW2aPegVDd9XfLxr1EobuOdtvM5J5v7HZ/zOSeQHm/eYFQ59zq4mVQ59T/TGgkiRJkjQwNTHqFQyX31BJkiRJ0hQZUEmSJEnSFLnlT5IkSdJAFFBjdiiFGSpJkiRJmiIDKkmSJEmaIrf8SZIkSRqMgglP+ZMkSZIk9cMMlSRJkqSB8VAKSZIkSVJfDKgkSZIkaYoe1YAqyXlJdljLs2VJFidZ1K7nJZmZZMkGzrlZko8mWdLGvzLJXo9wjBlJLmjrOnFD1iNJkiSNiwImarTXsD2q31BV1Usm1yUJkHZ7RFXd3vVs5gCmPRHYFdi/qiaS7A6s7Ldzki2AA4Etq2r2ANYjSZIkaZoaWIYqydeS/CjJtUnmtbplSXZumafrk5wB/BjYo4/xtkrymZZluirJEa3+vCT7t/JVSf68ld+f5PeBJwG3VdUEQFUtr6o7W5v7usY/PslZrXxWkv8vycXAJ4HPA7Nbhuopg/o3kiRJkqa1gpqokV7DNsgM1euq6o4kWwNXJvnqpOdPA15bVW8G6CSquDjJauDBqnrOpPZvAaiq/ZI8HfhOkqcClwCHJlkGrAKe39q/gE4gtBL4fpJDgQuBz1fVVX2s/6nA0VW1OsnhwKlV9dJeDVvAOA9g11137WNoSZIkSdPRIL+heluSq4HL6GSg9pn0/OaqumxS3RFVNbtHMAWdAOlzAFV1A3AznaDnUuCw9vyfgW2TPBaYWVU3VtVyOsHbu4AJ4MIkR/Wx/q9U1ep+XrSqzqyqOVU15/GPf3w/XSRJkiRNQwPJULWMztHAIVV1f5IFwFaTmvX9HdOaYddSfyUwB/gp8F1gZ+APgB+taVBVDwLfAr6V5OfAy+hkq7pzgBu6PkmSJEmTjNnPUA0sQ7U9cGcLpp4OPHcAY14CnAzQtvrtCdxYVb8GbgFeSScbdilwavtLkmcl2bWVNwP2p5PdAvh5kn1b/csHsEZJkiRJY2xQ31B9G3hjkmuAG+kEOhvqDOATSRbT+VZqbss8QSd4OqoFcJcCu7c6gCcAn0wyo91fAXysld8JfJNOQLYE2HYA65QkSZLUTIzi7PIRGkhA1QKdF/d4NLP9vR2YNanPzMmNq2rZmnZV9Stg7lrmew/wnlZeQdf2wKr6Np0Ar1e/+cD8HvVzJ90vABb0GkOSJEmS1nhUf9hXkiRJkqazR/WHfSVJkiSNj6qixuxUCjNUkiRJkjRFZqgkSZIkDUxNjHoFw2WGSpIkSZKmyIBKkiRJkqbILX+SJEmSBmbCQykkSZIkSf0wQyVJkiRpYDw2XZIkSZLUFwMqSZIkSZoit/xJkiRJGogqmJhwy58kSZIkqQ8GVJIkSZI0RW75kyRJkjQwY3bInxkqSZIkSZoqM1SSJEmSBqY8lEKSJEmS1A8DKkmSJEmaIrf8SZIkSRqIqmJizE6lMEMlSZIkSVO0yQdUSXZI8uau+82SfDTJkiSLk1yZZK9HOOaMJBckWZTkxMGvWpIkSZqeaqJGeg3bJr3lL8nmwA7Am4EzWvWJwK7A/lU1kWR3YOUjGHML4EBgy6qaPeAlS5IkSZpGNoqAKslpwM1VdUa7fy9wL/AbwIuBAj5QVeckORz4C+A2YDZwDfCUJIuA77b626pqAqCqlnfNc19VbdvKxwMvraq5Sc4C7qATSC0DDgV2aWO+oqr+7VH9B5AkSZK0SdooAirgbODDPJRleiVwGnAMcACwM3Blkkva84OBWVV1U5KZrTwboGWkvp/kUOBC4PNVdVUfa3gqcHRVrW5B26lV9dJBvJwkSZI0LvwdqhH4/9u7+2C7qrOO498fLwOUMDBAaIkgVDsDtAmkkJSiFhnKdFCr0ApNKXUa6xBg1Jm2IwgtOqhtHSQztIUZbKqA06JQAgyFSsNrAFHy0hZIUhDUJMpLEQEp4a3knsc/7k65hnuTw8nJ2eGe72fmDHvvs9Z61rrDHzw8a6/TJDz7JJmW5DDgOUarT/9QVSNV9RRwFzC76bK0qlZPMNZjwEHAuUAHuD3JB7uYxjVVNdLNfJPMS7I8yfJnn322my6SJEmSJqFtpUIFsBA4idFtflcBv7iJtpt8J6qqXgVuBm5O8hRwIqPVqrHp8s5vZsyNxl8ALACYMWPGcKXgkiRJ0kQKhqxAtW1UqBpXAR9nNKlaCNwNzEmyfZKpwNHA0nH6vQDstuEmyeFJpjXX2wGHAmubr59Kckjz/CNbbSWSJEmShsI2U6GqqlVJdgMer6onk1wPHAU8wGhl6eyq+nGSgzfq90ySe5OsZLQqdTvwjSQ7NU2WApc01+cANwH/BawEpmz1hUmSJEmatLaZhAqgqmaMuS7grOYzts1iYPFGzz6x0VDfm2D8hYxWvzZ+PndzMSRJkiRtWuGhFJIkSZKkLm1TFSpJkiRJb2XF6Eaz4WGFSpIkSZJ6ZEIlSZIkST1yy58kSZKk/ijoeCiFJEmSJKkbJlSSJEmS1CO3/EmSJEnqG0/5kyRJkiR1xQqVJEmSpL4ooIbsUAoTqi0W1mfHgUc9eP/1A485zPb496WtxN19h8H/u9W25fvPaXsKAzdy1dmtxd7+pItaidvmmtvywm+8u+0pDNxxHziy7SkM3JG779r2FAZuyZfvbSXu9GXt/bfQ3SPHDjzmOnYbeEx1xy1/kiRJktQjK1SSJEmS+qOGb8ufFSpJkiRJ6pEVKkmSJEl9UnQ8Nl2SJEmS1A0TKkmSJEnqkVv+JEmSJPWNh1JIkiRJkrpihUqSJElSXxRQHkohSZIkSeqGCZUkSZIk9egtlVAlmZZk4QTfLU4yq7n+dJIVSR5MsjLJCT3EujDJqiQXbum8JUmSpKFQ0OlUq59Be8u8Q5Vkh6p6AjhpM+32A74AHF5VzyeZAkx9k3HWA6cDU6vq1S2ZtyRJkqTJa6tXqJIcmOThJH/TVIuuTHJcknuTPJrkfUl2TXJZkmVJfrihopRkbpJrktwI3NKMtbL5bpckVzVVqKuBXZqQ+wAvAOsAqmpdVa1u+oytYu2dZM0Ecb4D7AosSTJna/+NJEmSpMmiOtXqZ9AGVaF6F3AyMA9YBnwC+BXgt4DPAz8C7qiqTyfZA1ia5Lam71HAoVX1bJIDx4x5JvBSVR2a5FDgB83zB4CngNVJbgeuq6obu5jjz+IAJFlXVTPHa5hkXrMWpk2b1s36JUmSJE1Cg3qHanVVraiqDrAKuL1Gz1NcARwIfAg4J8n9wGJgZ+Dnm763bkhyNnI08C2AqnoQeLC5HgGOZ3Rr4CPARUnO72KOE8V5g6paUFWzqmrWnnvu1U0XSZIkSZPQoCpUY99D6oy57zRzGAF+u6r+dWynJEcCL25i3HFrek2ytpTRStetwOXA+cB6Xk8id96o26biSJIkSdqs8neoWrII+MMkAUjy3i763A2c2rSfDhzaXE9LcviYdjOBtc31GuCI5nqTh1tIkiRJ0uZsK6f8/QXwFeDBJqlaA3x4M30uBS5P8iBwP6MVKYAdgflJpgGvAE8DZzTfzQe+neR3gDv6ugJJkiRpyFVBdTptT2OgtnpCVVVrgOlj7udO8N3p4/S9ArhivPZV9TLw8QnCHjvBXB6mqWQ1zhsvTvNsygRjS5IkSRKw7Wz5kyRJkqS3nG1ly58kSZKkSaDTwm9BtckKlSRJkiT1yIRKkiRJknrklj9JkiRJfePvUEmSJEmSumKFSpIkSVJ/VFEeSiFJkiRJ6oYJlSRJkiT1yC1/W2jHF59hv3+5cuBxHz3tmwOPOdSu/WwrYUdW/1srcdv0k312ansKA7fT1Utai/3iq9u3ErfNNbflPx9r52/dprUzT257CgN343afaXsKAzd92fpW4o7MntFKXIA1Nzw88JivvpaBx+xFgVv+JEmSJEndsUIlSZIkqW861Wl7CgNlhUqSJEmSemRCJUmSJEk9csufJEmSpP4oD6WQJEmSJHXJCpUkSZKkvijKCpUkSZIkqTsmVJIkSZLUI7f8SZIkSeqbKrf8SZIkSZK6MCkqVEnOB9ZV1fwtHGcEWMHo3+Uh4FNV9dKWz1CSJEkaAgWdTqftWQyUFar/7+WqmllV04GfAme0PSFJkiRJ267WE6okuyb5bpIHkqxMMifJEUnuSvL9JIuS7Nu0PS3JsqbttUnetpmxP5lkaZL7k3w9yfZJzkzyV2PazE1y8Tjd7wHe1d/VSpIkSZpMWk+ogOOBJ6rqsKYy9D3gYuCkqjoCuAz4UtP2uqqaXVWHMbol7/cmGjTJIcAc4JeraiYwApwKLAQ+OqbpHODqjfruAPwao9v/JEmSJHWpOtXqZ9C2hXeoVgDzk1wA3AQ8B0wHbk0CsD3wZNN2epIvAnsAU4BFmxj3g8ARwLJmnF2A/66qp5P8R5L3A48CBwH3Nn12SXJ/c30P8LfjDZxkHjAPYP+9dn/TC5YkSZI0ObSeUFXVI0mOAH4d+EvgVmBVVR01TvMrgBOr6oEkc4FjNjF0gL+rqnPH+e5q4GPAw8D19frZji831azNzXkBsADg8F/4ueE6F1KSJEmaQFFUeSjFQCWZBrxUVd8C5gNHAlOTHNV8v2OS9zTNdwOeTLIjo9v3NuV24KQk+zTj7JnkgOa764ATgVPYaLufJEmSJHWr9QoVMAO4MEkHeA04E1gPfC3J7ozO8SvAKuBPgCXAWka3Cu42Zpzzknxmw01V7ZfkPOCWJNs1Y/8+sLaqnkvyI+DdVbV0q69QkiRJ0qTUekJVVYsY/12oo8dpeylw6TjPzwfOH+f51UxQgaqqD4/zbMpmJyxJkiRpfEUrB0O0qfUtf5IkSZL0VmVCJUmSJEk9an3LnyRJkqTJwy1/kiRJkqSuWKGSJEmS1CdFx9+hkiRJkiR1w4RKkiRJknrklj9JkiRJfVH+DpUkSZIkqVtWqLbQMzvsyzff/vmBx33524OPOcxWvdLO/2n5ydvWtxK3TWfdc07bUxi4O37pgtZiH/vPf9xK3DbX3JbZXz6q7SkM3HPfuKXtKQzcvANua3sKA3f3yLGtxF1zw8OtxAV45wkHDzzmTiNrBx6zV9XxUApJkiRJUhdMqCRJkiSpR275kyRJktQfHkohSZIkSeqWFSpJkiRJfVJUeSiFJEmSJKkLJlSSJEmS1CO3/EmSJEnqiwI6HkohSZIkSeqGFSpJkiRJ/VFQHQ+lkCRJkiR1wYRKkiRJkno0KRKqJHOTXNJcH5RkcZL7kzyUZEEP430gyapmjF36P2NJkiRpMiqq0+5n0CbjO1RfAy6qqhsAksx4M52TbA+cCsyvqsu3wvwkSZIkTRIDq1AlOaOp+NyfZHWSO5OckmRFkpVJLhjTdl2SLyV5IMl9Sd7ePP/NJEuS/DDJbRueb2Rf4LENN1W1oun7sypWc39TkmPGxPvzJEuAc4GPAX+a5Mqt8beQJEmSNDkMLKGqqr+uqpnAbEYTniuAC4BjgZnA7CQnNs13Be6rqsOAu4HTmuf/BLy/qt4LXAWcPU6oi4A7ktyc5LNJ9uhiersCK6vqyKr6IvAd4KyqOnW8xknmJVmeZPm655/uYnhJkiRpOFR1Wv0MWhvvUH0VuAP4X2BxVT1dVeuBK4GjmzY/BW5qrr8PHNhc7wcsSrICOAt4z8aDN9v0DgGuAY4B7kuy02bmNAJc2+0CqmpBVc2qqllTdp/abTdJkiRJk8xAE6okc4EDgD8Dsommr1XVhjfKRnj9Xa+LgUuqagZwOrDzeJ2r6omquqyqTgDWA9Obf45d79i+r1TVyJtcjiRJkqSxiqE7lGKQ71AdAfwR8MkarcUtAX41yd7NQRCnAHdtZpjdgceb609NEOf4JDs21+9VBa0mAAACGElEQVQA9mr6rAFmJtkuyf7A+7ZwSZIkSZImkSQnN6d9d5LM6qbPIE/5+wNgT+DOJADLGT0A4k5Gq1X/uOFkvk04H7gmyePAfcA7x2nzIeCrSV5p7s+qqh8neQpYDawAVgI/2LLlSJIkSZpkVgIfBb7ebYeBJVRV9bsTfPX347SdMuZ6IbCwub4BeEPSVVVXMHrIBVX1OeBz47QpRo9DH29uUza6nzvBXCVJkiRNoCiqM/iDIfqlqh4CaApAXZkUP+wrSZIkSW2YjD/sK0mSJKkFLz7/yKJ7bzpm75ansXOS5WPuF1TVgg03SW4D3jFOvy908QrSG5hQSZIkSeqLqjq+7TlsTlUd18/x3PInSZIkST0yoZIkSZIkIMlHkjwGHAV8N8mizfVxy58kSZIkAVV1PXD9m+ljhUqSJEmSemRCJUmSJEk9MqGSJEmSpB6lqtqew1takqeBtT123xv4nz5OZ1uP22Zs1zwcsYctbpuxXfNwxB62uG3Gds3DEXtL4h5QVVP7ORn1hwlVi5Isr6pZwxK3zdiueThiD1vcNmO75uGIPWxx24ztmocjdptr1tbjlj9JkiRJ6pEJlSRJkiT1yISqXQuGLG6bsV3zcMQetrhtxnbNwxF72OK2Gds1D0fsNtesrcR3qCRJkiSpR1aoJEmSJKlHJlSSJEmS1CMTKkmSJEnqkQmVJEmSJPXIhEqSJEmSevR/p4YRXaotNYkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# colors make the correlation matrix more readable\n", "size = 15\n", "fig, ax = plt.subplots(figsize=(size, size))\n", "cax = ax.matshow(cor_table,cmap=cm.get_cmap('coolwarm'), vmin=-1,vmax=1)\n", "fig.colorbar(cax, ticks=[1,0,-1], shrink=0.8)\n", "plt.xticks(range(len(gcm_variables.columns)), gcm_variables.columns, \n", " rotation='vertical', fontsize=12)\n", "plt.yticks(range(len(gcm_variables.columns)), gcm_variables.columns)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Analysis: Average Temperature Prediction\n", "\n", "### 3.1 Support Vector Machine (SVM) VS. Ridge Regression\n", "#### 3.1.1. Ridge Regression for Climate Downscaling" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "RidgeCV(alphas=array([1.00000e-10, 1.66810e-08, 2.78256e-06, 4.64159e-04, 7.74264e-02, 1.29155e+01, 2.15443e+03, 3.59381e+05,\n", " 5.99484e+07, 1.00000e+10]),\n", " cv=None, fit_intercept=True, gcv_mode=None, normalize=True,\n", " scoring=None, store_cv_values=True)" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ridge = RidgeCV(alphas=np.logspace(-10, 10, 10), \n", " fit_intercept = True, normalize = True, \n", " store_cv_values = True)\n", "ridge.fit(station_GCM_sdf[station_predictors], \n", " station_GCM_sdf[station_predictand])" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEaCAYAAADtxAsqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3XmcXHWd7vHP01v2hXQFhJCQpDuyXJXgBAa6gxfX67iMODLuQgBBFHTUYWacmTvjcu+MXh03XAdFQS+4IcriBrLIJEEkIIQlYBYTgYQkbEkIIUl3f+ePczopiu7qqqSrTi3P+/WqV586dZZvVyr91O/8zjk/RQRmZmalaMm6ADMzqx8ODTMzK5lDw8zMSubQMDOzkjk0zMysZA4NMzMrmUPDrE5Imi0pJLWN5rJm5XBoWNVJerukZZKekrRB0i8kLZT0NklrJalg+TZJmyS9bohtLZLUn24r/3FI9X6j0SXppPQP/t9nXYtZIYeGVZWkDwNfAP4dOAiYBXwVeAPwE2Aq8D8LVns1EMAvh9nsLRExseCxfoh9P+db9758E5fUWu46ZToNeDz9aVZTHBpWNZKmAJ8Azo2IKyJie0TsjoirI+LvIuIZ4IfAqQWrngpcGhF9+7DPtZL+QdJyYHvaahlq3pGSbpL0pKR7Jf1l3jYulvQ1ST+XtB14acE+3ippWcG8D0m6Kp1+jaT7JG2T9LCk84vUOx44BTgXmCdpQZFlb5L0SUm/k7RF0pWSphUs9g5Jf5L0qKR/zlv3OEm3pL/vBklfltQx4htqTc+hYdV0AjCWpEUxnEuAUySNgz1B83rgO/ux37cBrwWm5gXPnnmAgKuBa4EDgfcDl0o6PG8bbwf+DZgELC7Y/lXA4ZLmFSx/WTp9EfCeiJgEvAC4oUitbwKeAn4E/IrnBmihU4EzgEOAPuCCgtcXAocDLwf+VdKR6fx+4ENAjuTf5eXA+0bYl5lDw6qqE3i0WIshIpYAG4E3prPeDPwhIu4sst3j02/Mg4/VBa9fEBEPRsSOYeYdD0wEPhURuyLiBuAakmAZdGVELImIgbRFlF/z08CVg8un4XEESZgA7AaOkjQ5Ip6IiDuK/C6nAT+IiH6S0HmbpPYiy383Iu6JiO3AvwBvLjh89vGI2BERdwF3AUenNd8eEb+NiL6IWAv8J889LGj2HA4Nq6bHgFwJ/QjfYe837HeRtD6K+W1ETM17dBW8/uAQ6+TPOwR4MCIG8uatA2aMsI18l7E3ZN4O/DQNE0haD68B1kn6jaQThtqApJkkh74uTWddSdIye22R/ebXtQ5oJ2k9DHokb/ppknBE0vMlXSPpEUlbSfqY8tczG5JDw6rpFuAZ4OQRlvsO8PL0j+vx7D3Ms6+GupVz/rz1wExJ+f8fZgEPj7CNfNeSBOJ8kvDYU3NE3BYRbyA59PVTkn6bobyL5P/k1ZIeAdaQhEaxQ1QzC2reDTw6Qq0AXwPuB+ZFxGTgn0gO05kV5dCwqomILcC/Al+RdLKk8ZLaJf2FpE/nLbeOpN/ge8B1EfHIMJscLbcC24G/T+s5iaQf5fulbiA95HY58BlgGnAdgKQOSe+QNCUidgNbSfoThnIq8HFgft7jTcBrJXUOs847JR2VdqB/Arg8PbQ1kklpLU9JOgJ4bym/p5lDw6oqIj4HfBj438BmksMr55F8A893CXAYpXWAnzDEdRrHllHTLuAvgb8g+Zb+VeDUiLi/1G2kLgNeAfyooN/mXcDa9DDQOcA7C1eUdDwwG/hKRDyS97gKWMWz+1fyfRe4mOQw1FjgAyXWej7JYbRtwDeAH5S4njU5eRAms/ok6Sbg/0fEN7OuxZqHWxpmZlYyh4aZmZXMh6fMzKxkbmmYmVnJHBpmZlayurjXfi6Xi9mzZ2ddhplZXbn99tsfjYjpo7nNugiN2bNns2zZspEXNDOzPSStG+1t+vCUmZmVzKFhZmYlc2iYmVnJHBpmZlYyh4aZmZXMoWFmZiWri1Nu99WKDVt56IkdIy/YRNpaRFuraG9toT392daSN90qOlpbaMufbhGtLULyGD1mza6hQ+OyW//Ed3876qcpNyUJ2ltanhs4g8+HfK2F9pbCMBLj2lt570ndPG/K2Kx/LTMrU0OHxjkndfGWY2eOvGCTiIC+gQH6BoLd/QPs7g/6+gf2TO/uH6CvP9g9MMDuvmS5XYPzCpcfiIJlkuUGl9/VN8D2Xf3pMs9+bfNTO9m2s4/PvXl+1m+JmZWpoUNjxtRxzJg6LusyrMC//ew+Llr8R97/snnMyU3IuhwzK4M7wq3qzn5JFx1tLXzphpVZl2JmZXJoWNVNnzSGdx1/GD/9/cOs2fxU1uWYWRkcGpaJwdbGl29YlXUpZlYGh4ZlYvqkMZx6wmx+eqdbG2b1xKFhmTn7JXPTvg23NszqhUPDMpObmLQ2rrzzYVa7tWFWFxwalqmzXzKXMW2t7tswqxMODctU0to4zK0Nszrh0LDMnZW2Nr50va/bMKt1Dg3L3GBr46q71rNqk1sbZrXMoWE1YbBvw1eJm9U2h4bVhM6JYzi1x60Ns1rn0LCacfaJcxnX7taGWS1zaFjN6Eyv20haG9uyLsfMhlCx0JA0VtLvJN0l6V5JH0/nz5F0q6SVkn4gqaNSNVj9OevEOYxrb+WC633dhlktqmRLYyfwsog4GpgPvFrS8cD/Az4fEfOAJ4AzK1iD1ZnB1sbVy93aMKtFFQuNSAz2aLanjwBeBlyezr8EOLlSNVh9Ovslc93aMKtRFe3TkNQq6U5gE3AdsBp4MiL60kUeAmYMs+7ZkpZJWrZ58+ZKlmk1ZtqEDk7rSVobKze6tWFWSyoaGhHRHxHzgUOB44Ajh1psmHUvjIgFEbFg+vTplSzTatBZ6ZlUF/ieVGY1pSpnT0XEk8BNwPHAVEmDY5MfCqyvRg1WXwZbG9e4tWFWUyp59tR0SVPT6XHAK4AVwI3AKelipwFXVqoGq29nnTiX8e2tfNH3pDKrGZVsaRwM3ChpOXAbcF1EXAP8A/BhSauATuCiCtZgdWywtfGzuzfwB7c2zGpCJc+eWh4Rx0TEiyLiBRHxiXT+mog4LiK6I+KvI2JnpWqw+jfY2rjArQ2zmuArwq2mHTChg0W9bm2Y1QqHhtW8dy9034ZZrXBoWM0bbG38/O4NPPCIWxtmWXJoWF1498K5TOhoc9+GWcYcGlYXDpjQwaL0TCq3Nsyy49CwunHmwjlMHOPWhlmWHBpWN/JbG/c/sjXrcsyakkPD6sq7T3RrwyxLDg2rK1PHd3B672x+fvcjbm2YZcChYXXHfRtm2XFoWN3Jb22s2ODWhlk1OTSsLp25cA6T3NowqzqHhtWlwdbGL+5xa8OsmhwaVrfOXDiXSWPa+OKv3dowqxaHhtWtKePbOX3hHH557yPct96tDbNqcGhYXTuz130bZtXk0LC6lt/auHf9lqzLMWt4Dg2re2f2zmHSWLc2zKrBoWF1b8r4ds7oncOv7t3o1oZZhTk0rCGcsTBpbfhMKrPKcmhYQ5gyLmltXHvfRu552K0Ns0pxaFjDGGxtuG/DrHIcGtYwpoxr58yFbm2YVZJDwxrK6emZVF90a8OsIhwa1lAGWxvXubVhVhEODWs4p/fOYbJbG2YV4dCwhpO0Nua6tWFWARULDUkzJd0oaYWkeyX9TTr/Y5IelnRn+nhNpWqw5nX6wtlMHtvGF3zdhtmoqmRLow/424g4EjgeOFfSUelrn4+I+enj5xWswZrU5LHtvPvEufx6hVsbZqOpYqERERsi4o50ehuwAphRqf2ZFVrUO9ja+EPWpZg1jKr0aUiaDRwD3JrOOk/ScknfknTAMOucLWmZpGWbN2+uRpnWYPa2NjZx90NubZiNhoqHhqSJwI+BD0bEVuBrQBcwH9gAfHao9SLiwohYEBELpk+fXukyrUEt6p3NlHHtfPF6tzbMRkNFQ0NSO0lgXBoRVwBExMaI6I+IAeAbwHGVrMGa2+Sx7bx74Rx+vWITyx96MutyzOpeJc+eEnARsCIiPpc3/+C8xd4I3FOpGswAThtsbfhMKrP9VsmWRi/wLuBlBafXflrS3ZKWAy8FPlTBGsz2tDauv9+tDbP91VapDUfEYkBDvORTbK3qFvXO5puL/8gXfr2Sby06NutyzOqWrwi3pjBpbDtnnTiHG+7fxF0PurVhtq8cGtY0TuuZzdTx7b4nldl+cGhY00haG3Pd2jDbDw4NayqnnnAYU8e3+ypxs33k0LCmMtjauPGBzdzp1oZZ2Rwa1nQGWxtfdGvDrGwODWs6+a2NlRu3ZV2OWV1xaFhTeuMxyQ2Xf/MH3wzTrBwODWtKh0wdx5zcBJasejTrUszqikPDmlZPVye/++Pj7O4fyLoUs7rh0LCmtbA7x/Zd/b5mw6wMDg1rWid0dSLBklWPZV2KWd1waFjTmjq+g/9xyGT3a5iVwaFhTa23O8fvH3yCp3f1ZV2KWV1waFhT6+3Ksbs/+N0fH8+6FLO64NCwpnbs7Gl0tLb4EJVZiRwa1tTGdbTy4sOmujPcrEQODWt6vV057tuwlce378q6FLOaVzQ0JLVK8hje1tB6unMALF3tQ1RmIykaGhHRD7yhSrWYZeLoQ6cwcUybD1GZlaCthGWWSPoy8ANg++DMiLijYlWZVVFbawvHz53mloZZCUoJjZ705yfy5gXwstEvxywbPV05fr1iEw8+/jQzp43PuhyzmjViaETES6tRiFmWevP6Nd4ybVbG1ZjVrhHPnpI0RdLnJC1LH5+VNKUaxZlVy/MPmkhu4hj3a5iNoJRTbr8FbAPenD62At+uZFFm1SaJ3u5Olq5+lIjIuhyzmlVKaHRFxEcjYk36+Dgwt9KFmVVbb1eOR5/axQMeAtZsWKWExg5JCwefSOoFdlSuJLNs9M5L+jV8iMpseKWExjnAVyStlbQW+DLwnpFWkjRT0o2SVki6V9LfpPOnSbpO0sr05wH79RuYjZIZU8cxu3M8S30fKrNhjXRFeAtweEQcDbwIeFFEHBMRy0vYdh/wtxFxJHA8cK6ko4CPANdHxDzg+vS5WU3o6c7x2zWPeQhYs2GMdEX4AHBeOr01IraWuuGI2DB4AWBEbANWADNIrjC/JF3sEuDkfajbrCIGh4Bd/pCHgDUbSimHp66TdH56uGna4KOcnUiaDRwD3AocFBEbIAkW4MBh1jl78DTfzZs3l7M7s312wlwPAWtWTCmhcQZwLnAzcHv6WFbqDiRNBH4MfLDMlsqFEbEgIhZMnz691NXM9ssBEzo46mAPAWs2nFL6NN4ZEXMKHiWdciupnSQwLo2IK9LZGyUdnL5+MLBpP+o3G3ULu3Pc8ScPAWs2lFL6NP5jXzYsScBFwIqI+FzeS1cBp6XTpwFX7sv2zSqlpzsZAva2tU9kXYpZzSnl8NS1kt6UhkA5eoF3AS+TdGf6eA3wKeCVklYCr0yfm9WMY2cfQHurfOqt2RBKucvth4EJQJ+kZwABERGTi60UEYvTZYfy8rKqNKui8R1tHDPrABY7NMyeY8SWRkRMioiWiOiIiMnp86KBYVbvFnYnQ8A+4SFgzZ5l2NCQ9M686d6C186rZFFmWevt7iQCblnjU2/N8hVraXw4b/pLBa+dUYFazGrGiw6dyoSOVp96a1agWGhomOmhnps1lPbWFo6f2+nQMCtQLDRimOmhnps1nJ7uHGsfe5qHn/RNnc0GFQuNIyQtl3R33vTg88OrVJ9ZZnq7OwHc2jDLU+yU2yOrVoVZDTr8oEnkJnawZNWjvHnBzKzLMasJw4ZGRKyrZiFmtUYSPV05lq5+jIig/OtbzRpPKVeEmzWt3u5ONm/bycpNT2VdillNcGiYFdHbnQwBu3il+zXMoPjFfedL8oFca2qHHjCewzrHs3S1Q8MMirc0ZgBLJd0s6b2SctUqyqyW9HTluHXN4/R5CFiz4UMjIj4EzAL+hWR88OWSfiHpVEmTqlWgWdZ6uzvZtrOP5Q9vyboUs8yNNJ5GRMRvIuK9wEzgC8CHgI3VKM6sFvR0JY3sJe7XMCutI1zSC4FPAF8BdgH/VMmizGrJtMEhYN2vYTb8dRqS5gFvA94K9APfB14VEWuqVJtZzejt7uSSpevYsaufcR2tWZdjlpliLY1fAWOAt0TECyPi3xwY1qx6u3Ps6h/gtrWPZ12KWaaKhcb/An4REXfnz5R0oqSuypZlVluOmzON9lb5EJU1vWKh8Xlg6xDzd5B0iJs1jfEdbRwz8wCWrvKgTNbcioXG7IhYXjgzIpYBsytWkVmN6u3Occ/6LTz5tIeAteZVLDTGFnlt3GgXYlbr9gwBu9qtDWtexULjNklnFc6UdCZwe+VKMqtNR89Mh4B1v4Y1sWLjaXwQ+Imkd7A3JBYAHcAbK12YWa1pb23hz+d2ul/Dmlqx8TQ2Aj2SXgq8IJ39s4i4oSqVmdWgnq5Obrh/E+uf3MEhU32U1ppPsZYGABFxI3BjFWoxq3mDt0pfsupR/tqj+VkT8ngaZmU4/KBJdE7oYKk7w61JOTTMytDSInq6cyxe9SgRkXU5ZlVXsdCQ9C1JmyTdkzfvY5IelnRn+nhNpfZvVim9XckQsKs8BKw1oUq2NC4GXj3E/M9HxPz08fMK7t+sIvL7NcyaTcVCIyJuBnx3N2s4M6eNZ9a08Sz2qbfWhLLo0zhP0vL08NUBwy0k6WxJyyQt27x5czXrMxtRb3cnt655zEPAWtOpdmh8DegC5gMbgM8Ot2BEXBgRCyJiwfTp06tVn1lJerpybNvZx90eAtaaTFVDIyI2RkR/RAwA3wCOq+b+zUZLT1cngE+9taZT1dCQdHDe0zcC9wy3rFkt65w4hiMPnsxijxtuTWbEK8L3laTvAScBOUkPAR8FTpI0HwhgLfCeSu3frNJ6uzr5zm/X8czufsa2ewhYaw4VC42IeNsQsy+q1P7Mqq13Xo5vLv4jy9Y+wcJ5uazLMasKXxFuto+Omz2Nthax2NdrWBNxaJjtowlj2jhm1lSWenwNayIODbP90NOV4+6Ht7Dl6d1Zl2JWFQ4Ns/2wcF4uGQJ2jU+9tebg0DDbD0cfOpXxHa2+D5U1DYeG2X7oaGvhuDnTPG64NQ2Hhtl+WtidY83m7WzYsiPrUswqzqFhtp96ugZvle5+DWt8Dg2z/XTE8yYxbUIHS92vYU3AoWG2n1paRE9XJ0tWewhYa3wODbNR0NudY+PWnaze7CFgrbE5NMxGQa/7NaxJODTMRsGszvHMnDbO12tYw3NomI2S3q4ct6x5jP4B92tY43JomI2Snu4c257xELDW2BwaZqNkcAhYH6KyRubQMBsluYljOOJ5k3yrdGtoDg2zUdTbneO2tU/wzO7+rEsxqwiHhtko6u3uZFffALeveyLrUswqwqFhNoqOm9NJW4vcr2ENy6FhNoomjmlj/sypDg1rWA4Ns1HW050OAbvDQ8Ba43FomI2yhd05BgJ+6yFgrQE5NMxG2fyZUxnX3upbpVtDcmiYjbLBIWAXOzSsATk0zCpgYXeO1Zu388iWZ7IuxWxUOTTMKqCnO7mliK8Ot0ZTsdCQ9C1JmyTdkzdvmqTrJK1Mfx5Qqf2bZenI501m2oQOH6KyhlPJlsbFwKsL5n0EuD4i5gHXp8/NGk5LizhhbidLVz3mIWCtoVQsNCLiZuDxgtlvAC5Jpy8BTq7U/s2y1tud45Gtz7Dm0e1Zl2I2aqrdp3FQRGwASH8eONyCks6WtEzSss2bN1etQLPR0jvYr+FDVNZAarYjPCIujIgFEbFg+vTpWZdjVrZZ08YzY+o492tYQ6l2aGyUdDBA+nNTlfdvVjWSWNid45bVHgLWGke1Q+Mq4LR0+jTgyirv36yqero72fpMH/eu9xCw1hgqecrt94BbgMMlPSTpTOBTwCslrQRemT43a1g9XTkAH6KyhtFWqQ1HxNuGeenlldqnWa2ZPikdAnbVY7zvpO6syzHbbzXbEW7WKHq6cty29nEPAWsNwaFhVmG93Z3s7BvgDg8Baw3AoWFWYX8+t5PWFrHE96GyBuDQMKuwvUPAelAmq38ODbMq6O3qZPlDT7L1GQ8Ba/XNoWFWBT2DQ8CudmvD6ptDw6wKjpmVDgHr0LA659Awq4Ixba0cO2caS3yRn9U5h4ZZlfR2dbJy01Ns3OohYK1+OTTMqqS3O7mliIeAtXrm0DCrkqMOnswB49t96q3VNYeGWZW0tIgTujpZsupRDwFrdcuhYVZFvd05Nmx5hj96CFirUw4NsyrqTW+VvsSn3lqdcmiYVdFhnckQsB433OqVQ8OsiiTR293JUg8Ba3XKoWFWZb3dObbs2M1967dmXYpZ2RwaZlV2QlcngG+VbnXJoWFWZQdOGsvzD5roW4pYXXJomGWgtzsZAnZnn4eAtfri0DDLQG9Xjmd2D3DHuiezLsWsLA4Nswz8+dxptLbI96GyuuPQMMvApLHtHH3oFBa7X8PqjEPDLCO93TmWP7SFbR4C1uqIQ8MsIz1dOfoHglvXPJ51KWYlc2iYZeTFh01lbHuLD1FZXWnLYqeS1gLbgH6gLyIWZFGHWZbGtLVy7OxpXHffRg49YBwdbS20t7bQ0dpCe1sLHa3aM6+9tYWOtvS1dLq9VXTsmR58CElZ/2rWwDIJjdRLI8JfsaypvfaFB/ORK+7m//5sxahtsyMNj/wwGdP27LAZKoRaW/aGTeF4H4V3ySocDiSe9VrxdQtnxHOXsDzvO6mbF8yYknUZe2QZGmZN763HzeLkY2awq3+A3X0D7O4PdvUNJM/7B9jVt/dnMi/2zst7fe8yUWS9fnanr+/sG+CpnX15ywZ9AwOIvcFR2GApbL8Utmg07JMy17Vn2bG7ti4AzSo0ArhWUgD/GREXZlSHWebGtrcytr016zLMSpJVaPRGxHpJBwLXSbo/Im7OX0DS2cDZALNmzcqiRjMzK5DJ2VMRsT79uQn4CXDcEMtcGBELImLB9OnTq12imZkNoeqhIWmCpEmD08CrgHuqXYeZmZUvi8NTBwE/STvC2oDLIuKXGdRhZmZlqnpoRMQa4Ohq79fMzPafrwg3M7OSOTTMzKxkDg0zMyuZCi/5r0WSNgPr9nH1HFCPtyupx7rrsWZw3dVUjzVD/dZ9eERMGs0N1sVtRCJiny/UkLSsHm+IWI9112PN4LqrqR5rhvque7S36cNTZmZWMoeGmZmVrBlCo15vhliPdddjzeC6q6keawbXvUdddISbmVltaIaWhpmZjRKHhpmZlcyhYWZmJWua0JA0V9JFki7PmzdB0iWSviHpHVnWV4ykEyV9XdI3JS3Nup5ypO/x7ZJel3UtpRrqs1JvJJ2cfq6vlPSqrOsphaSTJP1X+lk/Ket6SiVplqSrJH1L0keyrqeYYf4OlvV5r+vQSP+RNkm6p2D+qyU9IGnV4D9iRKyJiDMLNvFXwOURcRbwl1Uqe7DGcmr/r4g4B7gGuKSadRYqp+7UPwA/rG6VzzUKn5Wq2Yf3+Dki4qfp53oR8JYKljtY237XTDIM9FPAWOChStWab5Tqfj7ws4g4Aziqlmsd6rNd9uc9Iur2AbwEeDFwT968VmA1MBfoAO4Cjsp7/fK86X8E5qfTl9VB7T8EJtfLew68AngryR+u19VL3UN9Vmq1VuCFJF8m8h8H5q33WeDF9VAz0JKudxBwab2810AncCNwA3B6Ldda7LNd6ue9Lm4jMpyIuFnS7ILZxwGrIhm3A0nfB94A3DfEJh4CDgXupMqtrnJrlzQL2BIRW6tZZ6Ey654ITCD5EO+Q9POIGKhiuXuMwmelasqpNSI+CTzn0J+SUc4+BfwiIu6obMWjU3OeJ4Axlaiz0Ci91+cDH023dTnw7VqtdTTU9eGpYcwAHsx7/hAwQ1KnpK8Dx0j6x/S1K4A3SfoacHWV6xzKkLWn02dSoQ/jKBiy7oj454j4IHAZ8I2sAqOIcj4rWSv22RjK+0laeqdIOqeShRVRVs2S/krSfwLfBb5c4dqKKfe9/iXwgfQzs7aCdQ2l3Pf4OZ/tcj/vdd3SGIaGmBcR8RhwTsHM7cDpVamqNEPWDhARH61yLeUYtm6AiLi4eqWUpeTPSg0o+h4/54WIC4ALKldOScqt+QqSL3JZK7fue4BTKldOUeXWOtTfwbI+743Y0ngImJn3/FBgfUa1lKtea3fdlVdPtQ6qx5qhvuqueq2NGBq3AfMkzZHUQdIRe1XGNZWqXmt33ZVXT7UOqseaob7qrn6t1ThDoYJnE3wP2ADsJkncM9P5rwH+QHJWwT9nXWcj1e66XWuj1FxvdddKrb5hoZmZlawRD0+ZmVmFODTMzKxkDg0zMyuZQ8PMzErm0DAzs5I5NMzMrGQODasISf2S7pR0j6SrJU2twD5OknRNmescUuq4AQXrTZX0vv3dzjDbfp2k30u6S9J9kt6Tzj9ZUsVutV1QwyJJh+Q9XyspV419W31xaFil7IiI+RHxAuBx4NysC5LUFhHrI2Jf7hM0FdgTGvuxncKa2oELgddHxNHAMcBN6csnM8z4DJJG+75xi4BDRlrIzKFh1XALeXfelPR3km6TtFzSx/Pm/4uk+yVdJ+l76S2nkXSTpAXpdE7S2sIdSDpO0tL0G/tSSYen8xdJ+pGkq4FrJc0eHMRGyUiId6aPzZI+KmmipOsl3SHpbklvSHfxKaArXfYzBdsZK+nb6fK/l/TSvH1fIemXklZK+vQQ780kkhuHPgYQETsj4gFJPSQDg30m3WdX+j78u6TfAH8jabqkH6fv5W2SetP9fkzJgD03SVoj6QPF3mNJpwALgEvTfY1LF39/3vtwRDn/4NbAsr403o/GfABPpT9bgR8Br06fv4rkm7VIvrRcQzK4zAKScU3GkfwhXQmcn65zE7Agnc4Ba9Ppk4Br0unJQFs6/Qrgx+n0IpJbLkxLn88mbxCbdN5hwP3pzzbSga7Sfa1Ka33WevnPgb8Fvp1OHwH8iWT0uUXAGmBK+nwdMHOI9+qbwCaS20S8g72DEV0MnJK33E3AV/OeXwYsTKdnASvS6Y8BS0nGpMiRBFJ7qe9x+nwt8P50+n3AN7P+TPlRG49GvDW61YZxku7/h37PAAACgUlEQVQk+eN6O3BdOv9V6eP36fOJwDySP2JXRsQOgLRlUI4pwCWS5pHcGro977XrIuLxoVaSNJYk1M6LiHXp4aJ/l/QSYICkhXTQCPteCHwJICLul7SOZAhQgOsjYku6r/tIgil//AMi4t2SXkgSducDryQJnKH8IG/6FcBR0p67Y0+WNCmd/llE7AR2StqU/g4LKe89HrxN+e0kQyOb+fCUVcyOiJhP8keyg719GgI+GUl/x/yI6I6Iixh6XIBBfez9rI4dZpn/A9wYSR/K6wuW215k218HroiIX6fP3wFMB/4srX9jkX0OKlb7zrzpfoYZwyYi7o6Iz5MExpuKbC//d2kBTsh7L2dExLYi+y1WZ7Hah63bmo9Dwyoq/Zb9AeD89Fv8r4AzJE0EkDRD0oHAYuD1af/AROC1eZtZC/xZOj1c5/MU4OF0elEptUk6F5gUEZ8q2M6miNid9k0cls7fRtIaGsrNJGGDpOeTHCp6oMQaJko6KW/WfJLDWCPtE+Ba4Ly8bc0fYXfF3uOR9mUGODSsCiLi9yQD3r81Iq4lORZ/i6S7gctJ/nDfRjIOwF0kh0WWAVvSTfwH8F5JS0mO0Q/l08AnJS0h6UcpxfnAC/M6w88BLgUWSFpGEgT3p7/DY8ASJacQf6ZgO18FWtPf5wfAovTQUCkE/L2kB9LDeR9nb+h9H/i7tHO9a4h1P5DWujw99FV09LUR3uOLga8XdISbPYdvjW41Q9LEiHhK0niSb+9nR8QdWdfVSPwe2/7ycUqrJRcquZhtLHCJ/5hVhN9j2y9uaZiZWcncp2FmZiVzaJiZWckcGmZmVjKHhpmZlcyhYWZmJXNomJlZyf4buN7i+qnKYnEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = plt.gca()\n", "ax.plot(ridge.alphas.ravel(), \n", " np.mean(ridge.cv_values_, \n", " axis = 0).ravel())\n", "ax.set_xscale('log')\n", "ax.set_xlim(ax.get_xlim()[::-1]) # reverse axis\n", "plt.xlabel('Regularization Strength')\n", "plt.ylabel('CV Error')\n", "plt.title('CV Error vs Alpha')\n", "plt.axis('tight')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimal Regularization Parameter = 0.0004641588833612782\n" ] } ], "source": [ "print('Optimal Regularization Parameter = {0}'.format(ridge.alpha_))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3.1.2 Support Vector Machine for Climate Downscaling" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "svr_rbf = SVR(kernel= 'rbf', C = 1.0, gamma='auto').fit(\n", " station_GCM_sdf[station_predictors],\n", " station_GCM_sdf[station_predictand].values.ravel())" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([10.178156, 10.247221, 10.389448, 10.494384, ..., 9.024728, 9.024996, 9.054723, 9.058166])" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "svr_rbf.predict(output_GCM_sdf[station_predictors])" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import IFrame\n", "IFrame(src='https://eath.maps.arcgis.com/apps/StorytellingSwipe/'\\\n", " 'index.html?appid=ca2c84c6611140f6883f9843c06d1d56#',\n", " width=\"100%\", height=\"800px\") " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 EBK Regression Prediction VS. Geographically Weighted Regression (GWR)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Let's use EBK regression to predict average temperature\n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\"\n", "input_feature = station_fmv\n", "out_ga_layer = \"EBK_pred_ga\"\n", "out_raster = \"EBK_pred_raster\"\n", "\n", "arcpy.ga.EBKRegressionPrediction(station_fmv, \"T_MONTHLY_AVG\",\n", " \"vort850_201203_ebk_surface;specHum850_201203_ebk_surface;\"\\\n", " \"zonalSurf_201203_ebk_surface\", out_ga_layer, out_raster, \n", " None, None, 95, None, \"NONE\", \"EXPONENTIAL\", 100, 1, 100, \n", " \"NBRTYPE=StandardCircular RADIUS=15 ANGLE=0 NBR_MAX=15 \"\\\n", " \"NBR_MIN=10 SECTOR_TYPE=ONE_SECTOR\")" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Now, we use GWR to predict average temperature and compare the \n", "# result with EBK\n", "\n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\"\n", "input_feature = station_GCM\n", "out_feature_class = os.path.join(outGDB, \"GWR_out_feature_class\")\n", "kernal = \"ADAPTIVE\"\n", "number_of_neighbor = 30\n", "prediction_location = \"US_Prediction_Grid_GCM\"\n", "prediction_output = \"GWR_pred\"\n", "\n", "arcpy.stats.GeographicallyWeightedRegression(input_feature, \n", " \"T_MONTHLY_AVG\", \"specHum850;vort850;zonalSurf\", \n", " out_feature_class, kernal,\"BANDWIDTH_PARAMETER\", \n", " None, number_of_neighbor, None, None, 9.22719999999999E-02, \n", " prediction_location,\"specHum850;vort850;zonalSurf\",prediction_output)" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import IFrame\n", "IFrame(src='https://eath.maps.arcgis.com/apps/StorytellingSwipe/'\\\n", " 'index.html?appid=fe637afce2a4436db5332d08704722e2#',\n", " width=\"100%\", height=\"800px\") " ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The predicted output using GWR tool is point type. To better visualize continouous value (temperature), we convert the points result into raster with point to raster method. \n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\" \n", "input_feature = \"GWR_pred\"\n", "output_result = \"GWR_pred_raster\"\n", "assignmentType = \"MEAN\"\n", "priorityField = \"\"\n", "cellSize = \"2000\"\n", "\n", "arcpy.PointToRaster_conversion(input_feature, \"PREDICTED\", \n", " output_result, assignmentType, \n", " priorityField, cellSize)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "### 3.3 Forest-Based Classification and Regression VS. Forest-Based Classification and Regression with Distance Features " ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# use Forest-Based Classification and Regression method to predict average temperature\n", "\n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\"\n", "input_feature = \"station_GCM\"\n", "predict_location = \"US_Prediction_Grid_GCM\"\n", "output_result = \"RF_pred\"\n", "\n", "\n", "arcpy.stats.Forest(\"PREDICT_FEATURES\", input_feature, \"T_MONTHLY_AVG\",\n", " None,\"zonalSurf false;zonal850 false;zonal500 false;\"\\\n", " \"vortSurf false;vort850 false;vort500 false;specHum850 \"\\\n", " \"false;specHum500 false;specHum2m false;\"\\\n", " \"seaLevP false;relHumSurf false;meridSurf false;\"\\\n", " \"merid850 false;merid500 false;geoHgt850 false;\"\\\n", " \"geoHgt500 false;airT2m false;airFlowSurf false;airFlow850 \"\\\n", " \"false;airFlow500 false\", \n", " None, None, predict_location, output_result, None, \n", " \"specHum2m specHum2m;airT2m airT2m;geoHgt850 geoHgt850;\"\\\n", " \"geoHgt500 geoHgt500;relHumSurf relHumSurf;specHum850 \"\\\n", " \"specHum850;vort850 vort850;vort500 vort500;merid850 \"\\\n", " \"merid850;merid500 merid500; zonal850 zonal850;zonal500 \"\\\n", " \"zonal500; airFlow850 airFlow850;airFlow500 airFlow500 \"\\\n", " \"airFlowSurf airFlowSurf; vortSurf vortSurf; meridSurf \"\\\n", " \"meridSurf;seaLevP seaLevP;zonalSurf zonalSurf\", \n", " None, None, None, None, \"TRUE\", 100, None, None, 100, None, 10)" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The predicted output using Forest-Based Classification and \n", "# Regression method is point type. To better visualize continouous \n", "# value (temperature), we convert the points result into raster with \n", "# point to raster method. \n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\"\n", "input_feature = \"RF_pred\"\n", "output_result = \"RF_pred_raster\"\n", "assignmentType = \"MEAN\"\n", "priorityField = \"\"\n", "cellSize = \"2000\"\n", "\n", "arcpy.PointToRaster_conversion(input_feature, \"PREDICTED\", \n", " output_result, assignmentType,\n", " priorityField, cellSize)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the distance features to be applied for the Forest() method later would require a coordinate system defined, the cell below is to use the existing arcpy method \"DefineProjection_management\" to have the projection defined." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Start Time: Wednesday, February 13, 2019 2:52:49 AM\n", "Succeeded at Wednesday, February 13, 2019 2:52:50 AM (Elapsed Time: 1.75 seconds)\n" ] } ], "source": [ "try:\n", " # set workspace environment\n", " arcpy.env.workspace = home_dir + \"/climate_output.gdb\"\n", " \n", " # get the coordinate system by describing a feature class\n", " dsc = arcpy.Describe(\"station_GCM\")\n", " coord_sys = dsc.spatialReference\n", " \n", " # set local variables\n", " in_dataset = \"Large_Water_Body\"\n", " \n", " # set workspace environment\n", " arcpy.env.workspace = home_dir + \"/climate_input.gdb\"\n", " \n", " # run the tool\n", " arcpy.DefineProjection_management(in_dataset, coord_sys)\n", " \n", " # print messages when the tool runs successfully\n", " print(arcpy.GetMessages(0))\n", " \n", "except arcpy.ExecuteError:\n", " print(arcpy.GetMessages(2))\n", " \n", "except Exception as ex:\n", " print(ex.args[0])" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The power of GIS is that we consider space, so next, we include \n", "# the distance of water body to our analysis\n", "arcpy.env.workspace = home_dir + \"/climate_output.gdb\"\n", "input_feature = \"station_GCM\"\n", "output_result = \"RF_pred_water_body\"\n", "distance_feature = home_dir + \"/climate_input.gdb/Large_Water_Body\"\n", "\n", "arcpy.stats.Forest(\"PREDICT_FEATURES\", input_feature, \"T_MONTHLY_AVG\", \n", " None, \"zonalSurf false;zonal850 false;zonal500 false;\"\\\n", " \"vortSurf false;vort850 false;vort500 false;specHum850 \"\\\n", " \"false;specHum500 false;specHum2m false;\"\\\n", " \"seaLevP false;relHumSurf false;meridSurf false;merid850 \"\\\n", " \"false;merid500 false;geoHgt850 false;\"\\\n", " \"geoHgt500 false;airT2m false;airFlowSurf false;airFlow850 \"\\\n", " \"false;airFlow500 false\",\n", " distance_feature, None, \"US_Prediction_Grid_GCM\", output_result,\n", " None, \"specHum2m specHum2m;airT2m airT2m;geoHgt850 \"\\\n", " \"geoHgt850;geoHgt500 geoHgt500;relHumSurf \"\\\n", " \"relHumSurf;specHum850 specHum850;\"\\\n", " \"vort850 vort850;vort500 vort500;merid850 merid850;merid500 \"\\\n", " \"merid500; zonal850 zonal850;zonal500 zonal500; airFlow850 \"\\\n", " \"airFlow850;airFlow500 airFlow500 \"\\\n", " \"airFlowSurf airFlowSurf; vortSurf vortSurf; meridSurf \"\\\n", " \"meridSurf;seaLevP seaLevP;zonalSurf zonalSurf\",\n", " None, None, None, None, \"TRUE\", 100, None, None, 100, None, 10)\n" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import IFrame\n", "IFrame(src='https://eath.maps.arcgis.com/apps/StorytellingSwipe/'\\\n", " 'index.html?appid=409a3bf73fca49c1a8cc77c20900f30d#',\n", " width=\"100%\", height=\"800px\") " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Conclusions\n", "In this study, we explored different means of incorporating spatial information in climate model downscaling." ] } ], "metadata": { "esriNotebookRuntime": { "notebookRuntimeName": "ArcGIS Notebook Python 3 Advanced", "notebookRuntimeVersion": "10.7.1" }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.7" } }, "nbformat": 4, "nbformat_minor": 2 }