{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import nivapy3 as nivapy\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr\n", "import os\n", "import matplotlib.pyplot as plt\n", "import toc_trends_analysis as resa2_trends\n", "from collections import defaultdict\n", "from scipy.stats import theilslopes\n", "\n", "plt.style.use('ggplot')" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Connection successful.\n" ] } ], "source": [ "# Connect to NIVABASE\n", "eng = nivapy.da.connect()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# TOC trends - October 2018 (Part 3: Climate data)\n", "\n", "This notebook performs some further cleaning of the trends data and also estimates updated climate trends. The workflow is based on that described [here](https://nbviewer.jupyter.org/github/JamesSample/icpw/blob/master/toc_trends_oct_2016_part3.ipynb).\n", "\n", "## 1. Site elevations\n", "\n", "I have manually transferred the site elevations from the previous work here\n", "\n", " ...\\ICP_Waters\\TOC_Trends_Analysis_2015\\CRU_Climate_Data\\toc_trends_sites_for_climate_update_elev.xlsx\n", " \n", "to the new site codes here\n", "\n", " ...ICP_Waters\\TOC_Trends_Analysis_2015\\update_autumn_2018\\toc_trends_oct18_stations.xlsx\n", " \n", "The code below adds these elevations to the database, so they are stored for future reference." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Read station data\n", "stn_path = r'../../update_autumn_2018/toc_trends_oct18_stations.xlsx'\n", "stn_df = pd.read_excel(stn_path, sheet_name='Data')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## Update stations table\n", "#with eng.begin() as conn: \n", "# for idx, row in stn_df.iterrows():\n", "# # Add new vals to dict\n", "# var_dict = {'elev':row['elevation'],\n", "# 'stn_id':row['station_id']}\n", "#\n", "# # Update table\n", "# sql = (\"UPDATE resa2.stations \"\n", "# \"SET altitude = :elev \"\n", "# \"WHERE station_id = :stn_id\")\n", "# conn.execute(sql, **var_dict) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Read gridded elevation data\n", "\n", "The code below compares the measured station elevations to the pixel elevations in the 0.5 degree global topography dataset used for the climate workflow. Differences between actual elevations and pixel elevations will be used to correct temperature values later in this notebook." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.6/site-packages/ipykernel_launcher.py:14: DeprecationWarning: Dataset.sel_points is deprecated: use Dataset.sel()instead.\n", " \n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvX2cHFWV8P/t6QzJJE7SwTEh6UkENCokvq0y+ltdlzVGxQcXd37mSMTwImtUNNGNi0yMRAgiA2jkRVEDogQh8cQdFx/fgI0v+AIGdVFB3QURyUwiIZKZhGQSMtP1/FHVk57uqu7q7uqe6pnz/Xzm0923bt17+nZNnbrnnntOwnEcDMMwDCNuNI21AIZhGIbhhykowzAMI5aYgjIMwzBiiSkowzAMI5aYgjIMwzBiiSkowzAMI5aYgjIMwzBiiSkowzAMI5aYgjIMwzBiyaSxFqCOWMgMwzCM+JAoVWEiKSh27txZ0XltbW3s2bMnYmnGLzZe4bGxKg8br/DEeazmzp0bqp6Z+AzDMIxYYgrKMAzDiCWmoAzDMIxYYgrKMAzDiCWmoAzDMIxYYgrKMMaYlp4eZnV0MKe9nVkdHbT09Iy1SIYRCyaUm7lhxI2mzZuZ8dGP0jQ4CMCkvj5mfPSjAAx2do6laIYx5tgMyjDGkOS6dSPKKUvT4CCt3d1jJJFhxIe6zqBE5GbgdGC3qi7yyq4G3go8A/wJOE9V+71ja4DzgWFglare6ZW/GbgWSAI3qar9NxuNyY4dvsXJCjeVG8Z4ot4zqK8Cb84ruxtYpKovAf4XWAMgIicDZwILvXNuEJGkiCSBzwOnAScDy7y6htF4zJvnWzwccqe9YYxn6qqgVPUe4Km8srtUdcj7eB/Q7r0/A9iiqodV9c/AI0CH9/eIqj6qqs8AW7y6htFwDK9fT6alZVRZpqWF/V1dYySRYcSHuDlJvBv4uvc+jauwsvR6ZQA78spf5deYiKwAVgCoKm1tbRUJNWnSpIrPnYjYeIWnaflyMkBi3TrX3DdvHpn165m2bBnTxlq4GGLXVnjGw1jFRkGJyFpgCLjNK/KLdOvgP+vzjVSuqhuBjdk6lQZOjHPQxThi4xWetrY29ixZAkuWjD5g4+eLXVvhifNYhQ0WGwsFJSLn4DpPLFbVrLLpBXIN9O1AduU4qNwwDMMYJ4y5gvI88i4C/lFVD+Yc+hZwu4hsAOYCC4DtuDOrBSJyAtCH60jxzvpKbRiGYdSaeruZbwZOBdpEpBf4BK7X3mTgbhEBuE9V36eqD4mIAr/HNf19QFWHvXY+CNyJ62Z+s6o+VM/vYRiGYdSehONMmESzjiUsrA82XuGxsSoPG6/wxHmsvDWokhl1LZKEYRiGEUtMQRmGYRixxBSUYRiGEUtMQRmGYRixxBSUYTQwlkvKGM+M+T4owzAqo6Wnx3JJGeMam0EZRoPS2t1tuaSMcY0pKMNoUIJyRlkuKWO8YArKMBqUoJxRlkvKGC+YgjKMBmV/V5flkjLGNeYkYRgNStYRorW7m+TOnQzPncv+ri5zkDDGDaagDKOBGezsNIVkjFvMxGcYhmHEElNQhmEYRiwxBWUYhmHEktBrUCLyLCAF9Kvq07UTyTAMwzBKKCgRWQS8F/g/wHNxE0w5IvIY8D3gS6r6u1oLaRhG/Wnp6TEPQWNMCTTxeenZbwd2Ae8C2oBjvNflQB9wm4hsqYOchmHUkWycv0l9fSQcZyTOnwWjNepJsRnU7ar6f33K9wI/9/6uEJHTayKZUYA90Rr1olicP7vmjHoRqKAClJNfvW9HJ44RhEWuNuqJxfkz4kA5ThL/ALwceFZuuap+KmqhjELsidaoJ8Nz5zKpr8+33DDqRSg3cxG5HvgG8DrgpJy/F9VONCMXe6I16onF+TPiQNgZ1FnAIlWt6m4oIjcDpwO7VXWRV3Ys8HXgeOAxQFR1r4gkgGuBtwAHgXNV9dfeOecAH/ea/aSq3lKNXI2APdEa9cTi/BlxIOxG3R3A4Qj6+yrw5ryyLmCbqi4AtnmfAU4DFnh/K4AvwIhC+wTwKqAD+ISIzIxAtlhjT7RGvRns7GT39u3s6u1l9/btppyMuhN2BnU+cKPnev5E7gFVvSdsZ6p6j4gcn1d8BnCq9/4W4EfARV75JlV1gPtEJCUic7y6d6vqUwAicjeu0tscVo5GxJ5oDcOYaIRVUK/AndG8DshdqXeA+VXKMFtVdwGo6i4RmeWVp3Fnbll6vbKg8gJEZAXu7AtVpa2trSIBJ02aVPG5kbJiBZkVK8h4H6d5f3EjNuPVANhYlYeNV3jGw1iFVVCfAt6qqv9VS2HySPiUOUXKC1DVjcDGbJ09e/ZUJEhbWxuVnjsRGcvxCrtXLC57ysKOVVzkHWsa4dqKC3G+b80NuXYedg3qABDalFcmT3imO7zX3V55LzAvp147sLNIuTHBCRv9oNGiJDSavOMR+w3GhrAzqHXANSKynqMKBABVzfifEppvAecA3d7rHTnlH/RCKb0KGPBMgHcCn8pxjHgjsKZKGYxxQNi9YnHfU5b/pJ44cCDW8k4E4n7NjFfCKqibvdf35pQlcE1rybCdeU4WpwJtItKL643XDaiInA88Diz1qn8X18X8EVw38/MAVPUpEbkMuN+rtz7rMGFMbMLuFYvznjK/iCG+9mviIe9EIc7XzHgmrII6IYrOVHVZwKHFPnUd4AMB7dzMUaVpjDH5T/xcfjksWVJ3OcLuFYvznjK/J3W/RVeIh7wThThfM+OZUApKVf9Sa0GMxsT3if+CC2i58sq6mz72d3WNkgX894qFrTcWBD2R53sHxUXeiUKcr5nxTLF0GxtE5LhiJ4vIcSKyIXqxjEbB94n/4EFau7vrLstgZycDV13FUDqNk0gwlE4zcNVVBYoybL2xIOiJPDNzZizlnSjE+ZoZzxSbQf0PsF1E/gD82Pu8H2gFXoC7lvRC4JM1ltGIMXGzzQ92doa6aYStV2+CntT3rV8fS3knEnG9ZsYzgTMoVf0S8DzgRu/1IuBzwEeBE4EvAs9X1ZvqIKcRU4Ke+M02XxnZJ/XhmTNxcE17zuTJYy2WYYwJRdegVPUIbhTzb9RHHKPR8Hvid6ZONdt8lSQOHRpZc0r291vuL2NCEnajrmH44mebH77hBruRVkGxPTeGMZEInbDQMILIt823tbVBTEOsNAJxW9czjLHCZlCGETNsXc8wXExBGUbMsNxfhuES2sQnIjNw3cqflVuuqj+IWijDmMhY7i/DcAmloETkXODzwNO4cfGyOLgu54ZhRIjtuTGM8DOoy4G3q+r3aimMYRiGYWQJuwY1CbirloIYhbT09DCro4M57e3M6uiw3DNGw1LqWo76Wq+0velr1jBn/nzmpNPMmT+f6WvKz+RTbt+59WcvXMjsRYsK3y9axOyFC0fKj1uwwJUxnWb2okW09PQU9Nu0eXPFMgYdr/c9KXnJJZeUrLR169YMsGTr1q0/Wbp0aVD0/7hzyf79+ys6cerUqRw8eLB0xQjJBmFNPvUUCaBp/34m/+hHDLe3M3TSSXWVpVzGYrwalYkwVqWu5XKu9TDjVen/zvQ1a5i2aRMJxyEBJByH5t/8hqY9ezj8hjdE8l1L1j90iCZvk3ax94kjR1wZvc9T7rqLyXfdRbK/f6TfxF13+fZb6e/RtGsXrddeG8k9qbW1FeDSUvUSjlNa34jIDuA44Bngb7nHVHV+WZKNHc7OCveRjEXq5FkdHb7h/YfSaXZv315XWcolzqmmIV6pu+M+VlFQ6lou51oPM16V/u/MmT+fxPBwQbmTTLLr8ceL9llp30H1o8Kv30p/DyeZ9B2fSu5JXsr3oEwyI4Rdg3pXWb0bVWObNWuDX3oQCyNUW0pdy1Ff6xW353PzLVoeQd+1/n/2a7/S3yNoHGr5HcLmg/pxzSQwfLEEabXBUnfXn1LXctTXesXtJZP+N+Fk6KThZfcdVD8q/Pqt9PcIGp9a3pNCOUmISLOIXCoij4rIIe/1UhE5pmaSTXBss2ZtsJlp/Sl1LUd9rVfa3oGzziJ/wcPxymvVt1/9SnCSSTLNzaPLAoI2V/p7HDjrrLrfk8I6SXwG+AfgQqAb+ClwLrBw6dKld9ZMumhpKCeJoZNOYri9nebf/pbE008znE6z79JLG+IpP84L/1O3bKHJ5zoYTqc58J731F+eGI9VVJS6lsu51sOMV6X/O4ff8Aaa9uyh+cEHwXEgmeTA8uXsu+KKyL5rqfqZVAqnpYXE4cOj38+ciTNlykg5TU1w5AjgJrMcuPJKDr/pTaP6dTZsYOD008uWMej4gZUrI7snRe0k0Qu8VFX/llPWBvxGVdNlSzc2NJSTRCMT5/HKX4MC9ylwrLKjxnms4oiNV3jiPFZRO0kENVSyA8OIExZGyDAah7AKaivwf0XkUuBx4LnAxwGtlWDGxKaWruAWRsgwGoOwCuqjuArp88BcYCewGfhkVIKIyL8B/4q7Lvk74DxgDrAFOBb4NbBcVZ8RkcnAJuAVuPuy3qGqj0UlS72I036cOGGu4IZhQHg382eAdd5f5IhIGlgFnKyqgyKiwJnAW4DPquoWEfkicD7wBe91r6o+X0TOBK4E3lEL2WqF3YSDaRRXcHvAMIzaEqigROR1qnqP9/71QfUiTLcxCWgRkSPAVGAX8Hrgnd7xW4BLcBXUGd57gG8AnxORhKo2TBimRrkJjwWN4ApuDxiGUXuKzaBuABZ5778cUCeSdBuq2icin8Zd3xrEDUz7K6BfVYe8ar1A1mMwDezwzh0SkQHg2cAolxURWQGs8Oq5qcgrYNKkSRWfG0Sxm3DUfdWbqsdr3jzwCy0zb15sxqb56qtJ+DxgpK6+mmkrVoRupxbX1njGxis842GsAhWUqi7KeX9CLYUQkZm4s6ITgH5cp4zTfKpmZ0h+3oMFsydV3QhszB6v1OWyFu6as4rs5o6ra2hYqh2vlgsv9HcFv/BCBmMyNnN27PA/sGNHWd89zq7AccTGKzxxHqu5IaNPhI0kcUdAeVSx1t8A/FlVn1TVI0AP8PdASkSySrQd1zkD3NnUPE+GScAM4KmIZKkLFikimMHOTgauuoqhdBonkWAonR6zfUpBFAtdYxhGNIT14vungPJTI5LjceDVIjIV18S3GPgl8EPg7biefOcAWUX5Le/zvd7xHzTS+hPYfpxSxN0VfH9Xl+8szx4wDCM6iiooEVnvvT0m532WE4G/RCGEqv5CRL6B60o+BPw3rmnuO8AWEfmkV5ZdC/sycKuIPII7czozCjnqTdxvwkYw9oBhGLWn1AxqnvfalPMe3PWeHRz1pKsaVf0E8Im84keBDp+6h4ClUfVtGJVgDxiGUVuKKihVPQ9ARH6uqjfWRyTDMAzDCL9R90YAEWkF2sjxolPVR2sjmmEYhjGRCaWgROQk4HbgpbjmvQRH3brDZ/MyDKNmWGQLY7wRys0cN3rDD3Fj4u0DZgJfwvWkMwwjj5aeHmZ1dDCnvZ1ZHR209FS/I6NYm9nIFpP6+kg4zkhkiyj6NYyxIqyCeilwkar2AwlVHcBNXnhZzSQzjAalFsqiVJvFQmcZRqMSVkEdArL5hPeIyHzv3GfXRCrDaGDCKIuWnh5mL1pE8+TJzEmnmb1wYVEFNv3ii4u22QjxCw2jXMIqqJ8A4r3/BvA94MdAVIFiDQOojWksTD/T16xh9qJFzEmnQymMYpRSFi09PaRWrya5dy8J3AXdZH8/M1av9u2zpaeHpv5+/za9cFkW2cIYj4T14pOcjx8DHgRacXMyGUYk1CtCuF8/0zZtGhXgMaswKul7uEicRXBnWIkjRwqONx054hvNvrW7Ozh1ddL1UbLIFsZ4JOE4pSMEicjLVPWBOshTS5ydFZo74hx0MY5UOl6zOjp8b+xD6TS7t2+PQrSi/fhRSd/5ChC8YLdePME57e0kAv7vnESCXb29o8qK1gd2ed9lInjx2f9ieOI8Vl6w2MDnrixhY/HdLSJP4rqa3257n4xaUK91lHLaq6TvUmGQgmZY2WN+ZYH10+mR9xbZwhhvhF2DOg7Xa+9FwAMicq+IrBSRWbUTzZho1GsdpZz2Ku17sLOT3du3s6u3l93bt49SHPu7unCamwvOyTQ3+5rkyq1vGOOFUApKVYdV9Tuq+i5gNnAtbhTxgKQ4hlE+9UpB4tePnwGtVgpgsLOT/g0bGJ45E8freziVYmDDBt8ZULn1DWO8EHYGBYCITAFOB94BvBLXu8+oM/XydKs31eaBCjsufv0cOPvsuiqAwc5OnnjwQY4cPsyuvj6eeOihon1l6+/q6wtV3zDGA2GdJN4CvBP4Z+D3uPmZtqjqX2srXqSMCyeJUgvwcWAsxiuqcam3o0Gcrq1GwMYrPHEeq6idJD4NbAZerqp/qkIuo0qKbQKNi4IaC6IYl3q5uRuGEY6w+6BOrrUgRjgaJWLAyEwkpDu3L4kEOA7OtGkkDhwYKXZ8ZkaB41JG/6b8DSNehI1mPhlYBywDnq2qM0TkjcALVPVztRTQGE2pTaBxoGnz5gJzW0V45udc5QSQGBxkxoc+BFDadTuRoKWnJ5SCaRTlbxgThbBOEtcAi4CzOOrw9BDw/loIZQRTL0+3akiuW1e9cipBUyYzKrbd/q4unEShSTvhOAUx8PIdKbJlBKzHxkn5G8ZEIqyCehvwTlW9F8gAqGofkC56lhEZ2ZtoatUqnMmTXY+zCjzd6sKO+uw+yJ3ZDHZ2BiqY3Bh4+RHBU6tXM2P1arfM59x85T9ePSgNI46EdZJ4Jr+uiDwH+FvkEhkF5C/eJ/v7ybS00H/ddfFSTFnmzYPHH695N/kzm+F0umQMvPyZXeLIEV/F5Hjt5XrxmROFYdSXsDOorcAtInICgIjMAT6H625u1JhGy/UzvH59gRkyajJNTQVmzVLmz7LWkhKJgggQjfY7GEajE3YG9THgKuB3wFTgYeBG4NIayWXkENfF+5aeHmasXu0bmduZNo3M1KkkDh4sudkhcCdeGV58UF0MvHz81p3i+jsYxnglrJv5M8CHgQ97pr09qlp6h28ZiEgKuAnXGcMB3g38D/B14HjgMUBUda+IJHDDLb0FOAicq6q/jlKeOBFHz72Wnh5mrFwZOAVPHDiAk0xy+LWvZfLPfjYqGrcDOFOnkhgcjHwzbLGAqX4pKZzmZhzcVBdZgpxO4vg7GMZ4JlBBiciJRc5rFXFTREUY2fxa4Puq+nYROQZ3pvYxYJuqdotIF9AFXAScBizw/l4FfMF7HZfEMddPa3d3SftwYniYyT//eUGqiAQwPHMmux9+uGby+RE0w/Ir81NycfwdDGM8U2wG9Qjuw24xC40DJKsVQkSmA68DzoWRGdszInIGcKpX7RbgR7gK6gxgkzeLu09EUiIyR1V3VStLHClluqoVxcL+hDZrZTK+xcm+Pua0txe069cnRPfdg2ZYYdobq9/BMCYqoWLx1RoReRmwETfO30uBXwEfAvpUNZVTb6+qzhSRbwPdqvpTr3wbcJGq/jKv3RXACgBVfcUzzzxTkXyTJk1iaGioonMblabNm0lecAGJgwdHypypUxm+4QYyy5bRvGABiYg89bLtAoV9HnMMOM6oda5cORqdiXhtVYONV3jiPFbHHHMMRBiLDwARmQekVfW+CuUqJsffAStV9Rcici2uOS+IIM/gUajqRlzFB+BUGjgxzkEXa8WstWtHKQrA/bx2LXuWLKHlwguLrkGVQ7bdkfe5x3weKnLlaHQm4rVVDTZe4YnzWM0NuW4b6v4iIvNF5GfAH4H/8sreLiI3VSzhaHqBXlX9hff5G7gK6wnPpT3r2r47p/68nPPbAXOlipBSHmuDnZ0MXH89Gc/JIPsXhIObzjyoTnLnzvIy3fb1FWyUtU20hjG+CPsA/CXgO0ArkLW13A1E8gjrpe3YISIv9IoW45r7vgWc45WdA9zhvf8WcLaIJETk1cBAI68/xfHGGia77WBnJ3997DH6r78ep6Wl6Hx9OJ1mV2/vqBTl+e2W4w2X4OhG2Wy4ovwoEdljhmE0JmEVVAfumk8G70FZVQeAGRHKshK4TUR+C7wM+BTQDSwRkYdxlWF2R+R3gUdxHTluBC6IUI66Etcbazkx//w2sAadV6xd30y3zc1kfNKdZ8lulLVNtIYx/gi7BvUE8Hzgf7MFInIyEFk8G1V9ADdLbz6Lfeo6wAei6nssiWuKh3I81oJMc37hggCcyZNxvO+cmTmTfevXF0Rs8PXiC4iXV8w0aJtoDaNxKSdh4bdF5Apgkogsw92jZI+nVRLn6ATFNr3mEriBNZ1m9/btI5/9st4mDh0K1edgZyezOjqKbpS1TbSGMb4IZeJT1ZuBjwJLgR3A2cDFqnpbDWWbEIRZ6wmi3mtXQf2FNQdWa4Yr1zxokcgNo7EJ7Wauqv8J/GcNZZmQVBqdoN6RtcP0lzXNMW8eAxdeWCBHtbPFMGbHoGN+8qdWrmT6xRez77LLbLOtYcSQWGzUrRPOzgrNZrXeT1AsYkMQQeauoTyzWlSU01/QeNVb5jB9g/tAMFY5teK8VyWO2HiFJ85j5e2DKrlRN4p9lkaVDHZ2snv7dnb19hakeAii3mtXgf319YU2lY1lNuBi42LefoYRT0xBNSjVrF1F2V8CQrvFD3Z2MnDVVQyl03XPBlxqXOLglGIYxmhMQTUo9Z6N+PWXpZwZSCWzxSgoJj+Yt59hxJFi6TbWh2lAVddFJ44RlrD7lCpZ3yrWX2rlypJ7kZo2b2bW2rUl+yxXtlL1ix3Pvk5ft46mvXtHfQc/xV7rqOqGYZQm0ElCRL4SpgFVPS9SiWpH1U4SUd3s64XfvqNqHQIC9yLNnMkTDz5IS08PqYsuGhX01a9PP9kcIJNK+XrVlfou5XzXMIoubGLDak2UcV7IjiM2XuGJ81iFdZIwL74QtLW1cWDjxshv9rWmFl5zLT09pHzSvGeamxnYsIHW7u5QfZbrVVfqu0T5XYvJFkX7ucT5JhJHbLzCE+exilxBichJwNuB2ar6QS+w62RV/W01gtaRqhRU04knjpmLdKXMaW8vyGYLblTxXb29Fbc7e+FCkv39BeVD6TTJnTv9+8SNLJGdtQSFLcptK3dcS32XKL9rUFt+VDuWcb6JxBEbr/DEeawidTMXkaXAPUAaN4oEuJHNN1QoX8MR55BEQdTK069pYMC3PKt8fEkkRgXEJVH82kzu3Dkq8gNN/pdqtr8ov2s555hzhWHUjrBefOuBJar6PmDYK/sNbvbbCUG93bqjoBaefi09PUWVxf6uLpypU0eVO1AwI0k4Dk4RJZWZMWNUlPfE8HBBLqmwUdLLJWxU9Xrt4TKMiUpYBTULVyHB0bx0pXLUjStK3QDjGOct6n1HWeeBxPBwwbHsWAx2djJ8ww1un7gmsEA15DgMz5zpq3hIJAri9iUAJ5n0/S6VfNeg38yvrf4NGxjYsGFM9nAZxkQlbCy+XwHLgU05ZWcC8Vx8qQHF3LrrHRevHMJGJA9DUN4nJ5lk4KqrANfBIDs+mVTKd60qSzbaeUtPD9MvvpimbN2mJpr27vU/KZMJXPMp57uW+s2KRVU3DKM+hJ1BrQI+KSI/BqaJyJ3AZcC/1UyyGBK0yTROyfKqmcmVOjdwvS2TAShIvNhURDnlm8cShw+TwJ0lNR04EDjrisqkGqffzDAMf8Km2/gj8CLg88DHga8AL1bVh2soW8NQLweKUgqk0uy8LT09zF60iNTKlaPOTa1cyeyFC0fOD1IOmVSK1KpVviY5P7IzrmIK3refCNd8GtHpxTAmGqFMfCLybFX9G6B55c9T1T/VRLIGIjBhX4QOFGHTXZSTnbelp4cZ3qZaP2WSAJL9/SP9+KUGcZqbSQwMhHfLBpwpU0bJkAyx5yh7XmrVKlq7u6veJF2P38wwjOoIa+J7UEROyy0QkfcDv4hepMajHnHxgpRP6sMfHpnhlDMryG64bQpQTvn9ZJVcvvNAZto0mjwTXxiyJrzU6tVMX7PGdboIdWKC5N69R2eGq1cze9Giip1SDi1eXOBFaF55hhEvwjpJvBu4SUTuwN37dD0wF3h9rQRrJMLGxauGIOWTGB4emeGUMyto7e4uiAYRpv9854E57e2h28glceQI0267zdcjMB8nkSiYoTUdOQKeI0W5TiktPT1M3bp1VJtOIsHBpUvNCcIwYkTYNajvAS8GXgv8D/A34JQGiiJRc2odpbuY6Sk7wylnJlfuWktN9oGFUU4AIcyH5Tg4+M1GE47DlG3bQp1vGEZ9CLsG9Szg08AM4LO4M6pzgY1RCiMiSeCXQJ+qni4iJwBbgGOBXwPLVfUZEZmM6/L+Clxl+Q5VfSxKWcYSv2Cmfus/uSR37ixrJhc02/KjmOlrf1dXYITzkiSTJZVUJpXCmTYtlKxhla45SBhGYxB2Deo3QDPwElX9d1zT3koR+U7E8nwI+EPO5yuBz6rqAmAvcL5Xfj6wV1Wfj6swr4xYjjEjyBMPYOCqq3CSSd/zsjOZYjO5XC/Apv7+orusnZy/hDc78VvnGezs5MDZZxe0VWrO4zQ3c+Css4rmaMo0N7PvsstK5nLKEnY214hRQQxjIhJWQa1R1eWqOgCgqg8Ap+Ca+yJBRNqB/wPc5H1O4CrCb3hVbgHe5r0/w/uMd3yxV7/hKeWJ13/NNWU7ZLT09DB74cJRbuT5e43yFUoi76+Yy/q+K66g//rrA6NHOLjKJqvwhmfOpH/DBvZdccUop4vhmTMZTqWORmrYsGFkzWtUvVQKp4qwQ2OZet4wjPCEMvGpqvqUHQJWRyjLNcBHcYPQAjwb6FfVIe9zL26wWrzXHZ4cQyIy4NWPZ+jeHErlIiplfirXIcMvt5EfYbR7MZf1rCLxS1WRAIZnzfKN+h42+kN+vWpyc9XDqcUwjOopllF3o6qu8N5vCqqnqmcHHQuLiJwO7FbVX4nIqV6x3z3TCXEst90VwApPTtra2iqSb9KkSRWfm0vT5s0kc5L5TerrI3XRRbS2tpJZtsytNG8ePP544cnz5h2VYcUKMitWkHXunub9+dF89dUkQmyCDUuyr6/oWBRTsFGM4QhljEE+TZs3k7x1VSeMAAAgAElEQVT6ati5E2bOJHnoEKlVq0hdfTXD69cf/S3qQFTX1kTBxis842Gsis2g/pzzvtabcV8D/LOIvAWYAkzHnVGlRGSSN4tqB7J3v15gHtArIpNwnTeeym9UVTdy1JHDqTQ3SlR5VWatXTsq0yzgfl67lj1LlgDQcuGF/okRL7yQwQpkmLNjR3VC55NMFh2LWUVc3eOQm2Yk4G12fJ/KuWwef5ym97+f/fv31202FeecPXHExis8cR6ruSHXewMVlKpekfP+0ghkCkRV1wBrALwZ1L+r6lkishU3SeIW4BzgDu+Ub3mf7/WO/0BVYx9ZPYz3WNTmp3K89cI1WNzrbn9Xl2/K97is75QKq1TMjGkYRn0Ju1EXEXk9sAx3g+5OYIuq1nrjyEXAFhH5JPDfwJe98i8Dt4rII7gzpzNrLEckhN1IW20E8tz1mUwq5XriVdzaaDIzZ46KWJ6vPAc7O2ltbYW1a2O5vhPGldzczQ0jHoRK+S4iq4Eu3CCxfwHmA+cBV6nqZ2oqYXRUlfI9iqmyn8NCpqUl0rxCfn1EpaAczxOvKScChZ/8cTYt+Dlx5JOfbr6WxHms4oiNV3jiPFaRpnwHPgK8XlUvUtUbVLUL1wX8I5WLOH4IE2V8VkcHqVWrcKZMGe1KHXHSO98oCT71su7gQTiJhFsnmcSBo3H38sIjNVqKilJ7qmpljoxjQkvDiDthFRTAI3mfH2UCZdQNolSKi/zjyb17SRw+TP9119UkJFKYyODgZaedMiUwo23/ddexq6+PXY8/zq6+PvZ3dQXmd6qnSSy7p2tOOs2cdJrZixaVdbMf7Ozk4NKlBco5uz+rFllyK02DYhgTnbAmvvcApwKXcNSD7mLgx8DN2XqqGj6sdf2piYkvyGSUNROVOh41c+bPDxWAFdyb8q6+vsA9RSPlfX3gE7A1S/53Kce0UM5+ppaeHmasXl0wi3Oam+n3NvWGIerfpNR3KNZf5tFHY2uGiSNxNlvFjTiPVVgTX1gniS95r8sYvaRxlncs4ZX7x+EZx5TyzKtnMsPW7u5QAVjz8XPKKFjLClBO1ZjEwuS4yqW1u7tAOYEbGb0cz7sof5Mw36FYf3F+ojOMsSasie+EnL8TAz6fWAsB406puG71iPs2yoRUxnnO5MmB5rIwWW4dqMokVm7a9WIKpBzlEuVvEuY7WOw/w6iMsKGO/lJrQRoVvyjjubOKUsejoJQyCfLiSxw+TNPhwyOfk3v3klq5kub77w91wx9Op6taryk2s/AzmxXb01XOzT7K3yTMbKxYf2GjXxjGRKQcJwnDB78ss7mzCr/jB5cupbW7u2qPrqxnWDHHiGIu5kFp3qfdeiuZVKpo31Eo2SClkkmlfJ0KDi1eTCYvSCy4a1DlyFLqN4viO+SWR9mfYUwkQjlJjBPqug+qmONBFHuhwgaBrXQP1HAqRdPAgK9jhJNM0n/NNYHyhh2voLFwJk8m6eMxOJROs7+ri+kXXzziUZiZOZN969eP2c2+2t8zzgvZccTGKzxxHquonSSMMii2cF4qnUZYwqwRAaGSAvrRNDDAgeXLmXbrraOUVJQbi4PCOqVWrfKtn03KGKeZh0VGN4zaYTOoEJT7JFLMrTi5c6f/rCSRYFdvb+g+5rS3B7p9Z8m0tHBw6VKmbt0aTpnlybp7+/aK0lpU++RWb9f8sSTOT7lxxMYrPHEeq6pnUCLyE0JsxFXV15Ul2QSg2MJ52Hh8UGgmPLR4MVO2bSvpwJD90ZwpUzhyyikMnHIKqQ9/OPT+qNz1paAZSzX5mEpRD8cSwzDiTzEniZtwg7J+GfgRrhv5T4CvAffgupb/sMbyxYZyQtUUWzgPk821paeH444/flQG3El9fUzbtGnkc7HZUzYLbnLvXlKrVtF8//2QKb7jJhv6KMwCfq0jI5hTgWEYED6SxH3A+ar6UE7ZycDNqvrqGsoXJRWb+GbdfTdN739/6IXwUgvnxWYfLT09zPjQh2gqoVDKwUkkyMyY4et4kGV45kyeePDBUO2VMsHF2bQQN2ysysPGKzxxHquonSROojBp4Z+BF5UnVmOSXLeuICttqfTnELxwXmyhv7W7O1LlBLizrUSCTEtL8FpUGWuR9YqOYRjGxCasgvox8FURuZijsfguwTX5jX8CstIWuyFX6m1Wq5t8U38//dddR2rlSt/HlqaBgYKyoJleOetohmEYlRJ2o+653utDwAHgd7jTs/NqIFP8mDfPt7gWN+RibVbjbzk8d66rXNLpUP0GrTNNX7OGxIEDvhHQwzoxWOoJwzDCEEpBqepTqnomMAWYA7So6jJVjaeBM2KG168v6dgQFfu7usg0+f8s2Yi85ZIfeinMdwnarzXt1ltJ9vePzMLKTVNhqScMwwhL6FBHInISsBa4WFUzIvJCEXlJ7USLD5lly+rmVTbY2cnAtdcWVUTFEg2CG/onKCliWA+5IFNjvvdgAnCmTg09FuUGiDUMY+ISSkGJyFJc1/I0cLZX3ApsqJFcsWOws5Pd27ezq7e3JokG8/sKNMWl0xxYvrww4Z6XAXconaZ/wwaeeOihQFnDfJdyzJflrJvVwsHCTIaGMT4JO4NaDyxR1fcB2d2evwFeWhOpJhBBN9diprh9V1xB/3XXjZoFZTPgRqU8/foPmrmVo8yiTj1hJkPDGL+EVVCzcBUS5AQqoLp1+3FJOU/zfjfX1OrVzF64kNSqVW5K9iKmutxZEFDxLKKlp4fZixYdzQu1cCFAgSnwwPLlVa/FhV0DC4uZDA1j/BLWzfxXwHJgU07ZmcD4CoxWJUFBYpvvv38kRFHWXRvwDT+UOHJkZENtcu/ekSeApr17mb5uHamVK48GgPVeM6kUTfv2kfD2T03q62PG6tWAf2bafJlTq1eTyMlUm+zvZ8bq1Qxs2FAQ++7IKadUFeIo6uCqtifLMMYvYSNJvAi4C3dz7qtxQx+9AHijqj5crRAiMg9X+R0HZICNqnqtiBwLfB04HngMEFXdKyIJ4FrgLcBB4FxV/XWJbmoeLDYowkJ+ygunuRkHfNOXR8lwKsUTDz1UtE6QzFB5cNZ67mBv9MCycd7tH0dsvMIT57EKG0kirJv5H3GjRnwe+DjwFeDFUSgnjyHgI6p6Eq4C/IAXSqkL2KaqC4Bt3meA04AF3t8K4AsRyVEVgZ5v+Z+PHKm5cgJ3c24ps19UadThqHmzecqUujkrRGUyNEcLw4gfofNBqepBQGshhKruAnZ57/eLyB9wPQbPAE71qt2CO3O7yCvfpKoOcJ+IpERkjtfOmFEsJflYkZUnNydVrjktqjTqxXJg1drjEaozGY6V7IZhFCfQxDdW6TZE5Hhcl/ZFwOOqmso5tldVZ4rIt4FuVf2pV74NuEhVf5nX1grcGRaq+opnnnmmIpkmTZrE0NBQyXpNmzeTvOACEgcPVtRPPXDmz+fIw0cnvk2bN5NcsYJE3tg4zc0M33gjmWXLQrXbvGABiccfL9lfHBlL2cNeW4aLjVd44jxWxxxzDFQZLPamnPfPA96NO4v5CzAfOAe4uXIRCxGRZwH/AXxYVfeJSFBVvy9WoExVdSOwMXu8UntsaFvukiW0XHml+zTf11dRqvWoCEz1vmPH6O+yZAktn/kM09eto2nvXgAyqRT7LruMwSVLIOSYzQmIV1jQXwwZS9njvE4QR2y8whPnsZob0joTqKBU9Zbsey/dxpvy0m3cjqugPlG5mEcRkWZc5XSbqmYXAJ7Imu5EZA6w2yvPBqzN0g7Ewm0rGyTWL+WGH1mtGqRxg5Sc33lOIgGOw3A6TeLAAd/0Gn5muzCBbUslKKxFANlaJkXMxYLfGkY8CbsPqqbpNjyvvC8Df1DV3OgU38KdqeG93pFTfraIJETk1cDAWK8/5TPY2cnBpUtLhiWCYCVU7MwEQCIxEkHCSSY5sHz5yGbdfZdd5rvR9tDixSG/QY7jQDpNatWqoptho97fVM8NuFHLbhhGNIRVUNl0GwtEpEVEXoCrUKJKt/Ea3H1WrxeRB7y/twDdwBIReRhY4n0G+C7wKPAIcCNwQURy+NK0eXNJDy8/L7Ap27YVzXwLIYywxc71MusmgMTwMNM2bWL2okW09PT4KsiE4zB169ZQN/lRCoLCGHz5m2GjzoJbzw24lsHXMOJJ2H1QxwI3AJ1AEtctvAdY2UARzSvaB9XS00PqootGOT7kZ9MNyqCbGBwsqoCKmfCqIStfa3d3xXuEiu2PyuIkEuzq7S0oj8L2Pae93Ve5B/XZqMR5nSCO2HiFJ85jFXYfVCgFlUVEmoDnAE+qarRpX2tPRQoqzEbQwA26iUTJGVStGEqnSe7c6X+Txw06W2xtJ0hB5Pfhp+ii+Mdo9A24YYnzTSSO2HiFJ85jVXXKdxE5sch507Iedqr6aLnCNRJhQukEbmgdI+UEjCgf31lQIlHV/iio/RrN/q4u31mprQsZxsSh2BrUI8DDOa/Z97mf473BJQJKRd8utp4zlm7m2ZlRgaMEpdeTIDiaeTalR63XaGxdyDCMYm7mI8pLRM4D3gBcgrsP6rnAOtzwQ+Oa/V1dvmtQ+7u6jgZaHcOZkh9OIjHKbDf94otp8rLgBinN/Flg1EFdKyGM+7thGOOXsKGOLgMWqGrW3vKwiLwX+F/gq7UQLC4MdnbS2toKa9cW3KhndXSMigIeGxxnVFqO1u5uEj57onKpdH+UYRhGrQiroJpwI4r/IafsubgefeOezLJl7FmypKA8tikdEokRV3MoLaet7RiGEUfCKqjPAj8Qka8AO3CjOJzrlU9Y6hkctljEiXwSjkNrd/eIggqSM+vNd2jx4qO5psgJdWSzJ8MwxpCw6TauBs4DZgP/jJu36d2qelUNZYs1LT09NOUkFKw5ySSZmTPDV89RSEGREvqvv55DixczbdMmknv3jqxRZRMWWsoJwzDGknLSbXwf+H4NZWkYwsbZi5LE8DCJp5/GaW4ete4VuNk3edT6GuTwADDt1lt9z286cmTULKzW1CvunmEYjUOxfVBrVfVy7/36oHqquq4WgsUZvzA89aDpyBF3JjQ0NOI5GGjyy0sl7+fwMKujo6gHYr3W2Cwfk2EYfhSbQbXnvJ8XUCde/tV1YiydI0qFT8oynE6XrFPqe9QrmnexuHumoAxj4lJsH9T7YSS80a3Az1T1cL0EizNxzJybS1ivvGLfI9PcXDfPvjDROgzDmHiUdJLwYu7dYcrpKH5OB/UiMD9UMll2xIWgSBOZqVMZ2LChbrOXUtE6DMOYmIRNt3GPl3fJIC8Mz1gLg+eRd8019F93HQCpVasC04Lk4hdOqP/66/nrww/X1bRm+ZgMw/AjrBffX4DvicgduPugRu7LE9FJAo46HcyZP7/AISEKwqTiyO5jyt7IK3E0iEO0iDiEVTIMI36EVVAtwH9679uLVZwItPT0jMS3g+jzOjnNza47eU78P19yYu7N6uhoaEeDOChKwzDiRSgFparn1VqQRmH6mjVM27SpppHKM9Om0TQwULJebsQIczQwDGO8UVJBiUizqh7x3r+W0etWP1fVoVoJFzdaenoCN7ZGSVN/v7vRNoTpMKuAMjNmkPQLCNvUNCoun2EYRqNQVEGJyPuBvweWe0V3AX/z3k8FPgp8uWbSxYSmzZuZtXYtyb6+yJRTKbNgIuS6ljN1Kse+4x0j5ka/dmzTq2EYjUgpL76zgU/nfD6sqvNUdR6wGPjXmkkWE1p6ekhecAGTylBOxTz7HIorp6BjQW0mDhxg8k9/WlQ2v4SEhmEYcaeUgjpBVX+T8/n3Oe9/AxRLCz8uaO3uLu2skEfRmVGRY5W4rBdLQphLsq+PWR0dzGlvD+WCbhiGMdaUUlDPEpFp2Q+q+pqcY1OBaYWnjC9q4WRQiZmwatNiIuHOAh1nxAXdlJRhGHGmlJPEg8AbgW/6HHsz8FDkEpWBiLwZuBY3ceJNqhq5HcuZMoVEnQLDZmbOBMfxd3aoAieRKAgK20gu6IZhTExKzaCuAW4Qkbd5MfkQkSYR+Rfgc97xMUFEksDngdOAk4FlInJy1P0kDh2KuklfnOZm9q1fz77LLiPT3FxdW7ix9LLRIQiIWG4u6IZhxJmiCkpVt+A6SXwNOCQiO4FDwCZgg6purr2IgXQAj6jqo6r6DLAFOCPyXoqkowg8pcy6Q+k0/V7su8HOTgY2bCgZRinomAMcOPts/vrYY+zq7WX39u2Bkc0t1p1hGHGm5D4oVf2MiNwI/H9AG66b+b2qWnonaW1J44ZdytILvCq3goisAFYAqCptbW31ky4EztSpDN9wA5lly5hGzoLeihVkVqwgOWWKr4J0AI49Fp56qrDR+fM55sYbGfVNL78c54ILRjl7OFOnwuWX12RMJk2aFLuxjis2VuVh4xWe8TBWYSNJ7APurLEs5VLSG1tVNwIbs8f27NlTdifHtbSUvwZVZJOtk0iA44zE0BtcsgQC5JoVkA4jk0qx79JLC7L6ZlpaGLjwQgbz21uyhJYrryyMdVek72poa2ujkrGeiNhYlYeNV3jiPFZzQ1pvQqd8jyG9jE6k2A5EvqgycNVVpD70IRKZjO/x/H1LDsDwcIFjgoOnWC67LLRjwv6uLl8llNtG2ACrFuvOMIxGo5EV1P3AAhE5AegDzgTeGXUng52dTP/d72ja6E7ECpRRHiPHHadwtpSjIFp6enyVS375waVLmbJtm68SMqVjGMZ4pmEVlKoOicgHcU2PSeBmVY3c7b2lp4emr3zF155Yam9SwnEYSqfZvX17QZt+qTGa77+fqVu3jiqfunVr6ASEhmEY44mEU4GXWoPi7KzArXpWR0dV6d2dRIJdvb2h2nSSSd8YfH5KLs7E2fYdN2ysysPGKzxxHitvDapk/IGwGXUnLNXuFfJz5Q5sM8CxwvYrGYYxETEFVYJq9go5XkLB0G0mk5HLYBiG0aiYgirB/q4unAoiOziJBAeWL/ddOzq0eLHrQJFDpqWFA2edRaalpaDcT8kZhmGMd0xBlWCws5PhG29kOJUaSZVR6m8onab/uuvYd8UVBe219PQwdevW0S7oiQQHly5l3xVXMHDVVW4UCS9MkTlIGIYxUWlYL756klm2jD1LlhR43/nW9WY8QUqltbu74PyE4zBl2zb2Ya7jhmEYWWwGVQZ+yiWfUskBgxwezBHCMAxjNKagyiCsEilWL8jhwRwhDMMwRmMKqgzCKpFi9fZ3dZkjhGEYRghMQZWBn3LJp5SyGezsNEcIwzCMEJiTRBmMCtDa13c0arn36hdzL6gdU0iGYRjFMQVVJqZcDMMw6oOZ+KqgpaeHWR0dzGlvZ1ZHBy09PbE/1zAMo1GwGVSFBEUkB0rOsMbqXMMwjEbCZlAV4rcnqtQeqLE+1zAMo5EwBVUh1Wy4HatzDcMwGglTUBVSzYbbsTrXMAyjkTAFVSHVbLgdq3MNwzAaCXOSqJBRe6J27mR47txQe6DG8lzDMIxGwlK+hyDOqZPjiI1XeGysysPGKzxxHitL+W4YhmE0NKagDMMwjFhiCsowDMOIJWPuJCEiVwNvBZ4B/gScp6r93rE1wPnAMLBKVe/0yt8MXAskgZtU1XapGoZhjDPiMIO6G1ikqi8B/hdYAyAiJwNnAguBNwM3iEhSRJLA54HTgJOBZV5dwzAMYxwx5jMoVb0r5+N9wNu992cAW1T1MPBnEXkE6PCOPaKqjwKIyBav7u/rJLJhGIZRB8ZcQeXxbuDr3vs0rsLK0uuVAezIK3+VX2MisgJYAaCqWdfGiqjm3ImIjVd4bKzKw8YrPI0+VnVRUCLyX8BxPofWquodXp21wBBwm3fMz0fewd8s6buZS1U3AhvLFjgPEfmlqr6y2nYmCjZe4bGxKg8br/CMh7Gqi4JS1TcUOy4i5wCnA4tVNatseoF5OdXagexO26BywzAMY5ww5iY+zyPvIuAfVfVgzqFvAbeLyAZgLrAA2I47s1ogIicAfbiOFO+sr9SGYRhGrYmDF9/ngFbgbhF5QES+CKCqDwGK6/zwfeADqjqsqkPAB4E7gT+4VfWhGstYtZlwgmHjFR4bq/Kw8QpPw4/VRIrFZxiGYTQQcZhBGYZhGEYBpqAMwzCMWDLmThJxx8IqFSIijwH7cUNQDanqK0XkWNw9bMcDjwGiqntFJIE7fm8BDgLnquqvx0LueiEiN+N6pe5W1UVeWdnj43m3ftxr9pOqeks9v0c9CBirS4D3AE961T6mqt/1jk3Y8GciMg/YhLtlJwNsVNVrx/O1ZTOoIlhYpaL8k6q+LGefRRewTVUXANu8z+CO3QLvbwXwhbpLWn++ihueK5eyxse76XwCdxN6B/AJEZlZc8nrz1cpHCuAz3rX18tylNNED382BHxEVU8CXg18wPue4/baMgVVnA68sEqq+gyQDatkFHIGkH0KuwV4W075JlV1VPU+ICUic8ZCwHqhqvcAT+UVlzs+bwLuVtWnVHUvbsxKvxt5QxMwVkGMhD9T1T8D2fBnE+L/VFV3ZWdAqrof14s5zTi+tkxBFSdNYVildEDdiYQD3CUiv/LCSQHMVtVd4P4jAbO8chtDl3LHZ6KP2wdF5LcicnPO072NlYeIHA+8HPgF4/jaMgVVnKBwSxOd16jq3+GaED4gIq8rUtfGsDhB4zORx+0LwPOAlwG7gM945TZWgIg8C/gP4MOquq9I1YYfL1NQxSkWbmnCoqo7vdfdwDdxTSxPZE133utur7qNoUu54zNhx01Vn/A25WeAGzmaxWDCj5WINOMqp9tUtccrHrfXlnnxFed+LKzSKERkGtCkqvu9928E1uOGpjoH6PZe7/BO+RauuWYL7qLsQNYcMcEoa3xE5E7gUznmrTfi5Uob74jInJxr5F+AB733Ezr8meeV92XgD6q6IefQuL22TEEVQVWHRCQbVikJ3FyHsEpxZzbwTREB9/q5XVW/LyL3Ayoi5wOPA0u9+t/FdXN9BNfV9bz6i1xfRGQzcCrQJiK9uB5T3ZQxPqr6lIhchvuQBLBeVcM6EzQMAWN1qoi8DNfs9BjwXnDDn4lINvzZEF74M6+difB/+hpgOfA7EXnAK/sY4/jaslBHhmEYRiyxNSjDMAwjlpiCMgzDMGKJKSjDMAwjlpiCMgzDMGKJKSjDMAwjlpibuWFUgRdy5s9As5ftOex5HwNOVNV/rZVsOX39CPiaqt5U677y+v0i0Keql9Wg7ZNx48y9smTlwnNnAz8CXqaqh6OWzYgOU1BGzfHSc8wF5qrqnpzyB4CXAieo6mNjI93YoKqfyr6vVMnFCRE5F/hXVX1ttkxV31fDLi8DPl3Jiar6hIj8EDfC9/WRSmVEipn4jHrxZ2BZ9oOIvBhoGTtxChGRhIjY/0TM8cL5/BPwn1U0cxveBmAjvtgMyqgXtwJnc/SJ9Rzc5GufzFYQkcnA5YAAk3Hj/P2bqg56YVluxQ3ZMgn4GfA+Ve31zj0XWAc8B9gDfFxVb/OS3z1fVd/l1TuenNmKZ/76GW40g78DXiwiTwIbcHfhZ4CvAJ9Q1WEv99CVwLnAPo4GMvVFRC4CVgHTceOdXaCq2/Lkuser3u9F6FiiqveKyLuBC3ET1G0HVqjqXwL6ebUn88nAX4APqeqPAur6tuuZ5J5W1X/PqXsH8GNV3SAiXbiJBGfhRsNeq6rfFJGTgC8CzSLyNG4Sy5SIfBXoVdWPe229B7gIOBb4Ke7vt9M75gDvBz4CtAG3Ax9UVb9IAkuAX6vqoRw5H8PNCbUcN9DsFtwoC18FXosb9Xupl14C7/OJIvLcoDE1xh57WjTqxX3AdBE5ybvJvwP4Wl6dK4EX4Eaxfj5uCoB13rEmXEXxXGA+MAh8DkbiA14HnKaqrcDfAw8QnuW45p5W3Jv7LbihdJ6Pm9LgjUB2reg9uBlgXw68Enh7UKMi8kLgg8Apnlxvwg3dk082GnxKVZ/lKae34d5gO3GV7k+AzQH9pIHv4Cr7Y4F/B/5DRJ7jU7dYu7cD7/BivuE9FLwR92YP8CfgH4AZwKXA17y4eX8A3gfc68mf8un39cAVuA8fc3DHeUtetdOBU3DNvuKNlx8vBv7Hp/z/x1VeLwDeCnzP+65tuNfPqmxFz5T6iNeXEVNsBmXUk+ws6sfAH3EDewIjgTDfA7wkGxdMRD6Fe9Nco6p/w43inK1/OfDDnLYzwCIRedwLNFpOQNqvZmO3eQvop+Eqi0HggIh8FleBfQn3xnmNqu7w6l+BO/vyYxh3JniyiDxZ5jrbe4ErvJt/diw+FvDE/y7gu9nMs8DdIvJL3BlgfirvwHZxlZWDq4TuwVW+9+ZEr9+a087XxU2/3sHR4KTFOAs3Rl425fgaYK+IHJ8zLt2q2o87k/wh7oPK933aSgF/8ym/XlWf8Nr/CW4a+f/2Pn8TWJxXf7/XlhFTTEEZ9eRW3BvfCbjmvVyeA0wFfuWZucCNUp0EEJGpwGdxM39mozC3ikhSVQ+IyDtwZw5fFpGf4abG/mNIuXKTtz0XaAZ25cjRlFNnbl79QPOQqj4iIh8GLgEWelGkV2dv+CV4LnCtiOSaEBO4s8r8Pp8LLBWRt+aUNTNagZds1zPzbcFdK7wHNyL4yCxXRM4GVgPHe0XPwp2dhGEu8OvsB1V9WkT+5n2fx7ziv+bUP+i178de3NluPk/kvB/0+ZzfXivQX0pwY+wwBWXUDe8G+GfcJ/vz8w7vwb2JLFTVvoKT3bWJFwKvUtW/etGu/xsv+Zqq3gncKSItuKauG3FnAgdwFV+W43zazl3n2AEcBtoCPOp2MTqXzny/75pFVW/HTRExHXcGdiWuSTGo/1w5LlfV24q1n1P3VlV9T8i6xdrdjJstuRt3ve9fALwZ1o24s5B7vfW4Bzia/K5U1OmduMoRr71pwLPJmUWXwW9x1zArRkQm4Zpwf1NNO0ZtsTUoo96cD7xeVdi+bJYAAAICSURBVA/kFuYkp/usiMwCd21FRLLrEK24CqxfRI7FTcuAV2+2iPyzd9M7DDyNa14Ddy3qdSIyX0RmUCLvjWcevAv4jIhMF5EmEXmeiPxjtgqwSkTavTWarqC2ROSFIvJ6z/njkCf/sE/VJ3FNlCfmlH0RWCMiC722ZojIUp9zwZ3lvFVE3iQiSRGZIiKniki7T92i7XomsSeBm4A7PZMbwDRcJfSkd955wKKcdp8A2kXkmAAZbwfOE5GXeePxKeAXFW4vuBv4OxGZUsG5WTqAx8xBIt6YgjLqiqr+SVV/GXD4ItyF6/tEZB/wX7izJoBrcN3S9+A6XOSuTTThzrB2Ak8B/whc4PV3N/B13KfuXwHfDiHm2cAxuHmH9gLfwF3YB1eJ3on75P1roMevAY/JuLl69uCar2bhLtqPQlUP4nov/kxE+kXk1ar6TdzZ1hZvLB7EXRsrwFsPO8Nr+0ncWdKF+Px/h2x3M/AGXKWSPe/3uB6L9+Iqoxfjej9m+QHwEPBXEdlDHqq6DbgYdx1xF66n3Zl+36cU3jrTD3C/c6WchausjRhj+aAMw2g4xI0kcQvQEeCKXuzcWbiOOi/PdVU34ocpKMMwDCOWmInPMAzDiCWmoAzDMIxYYgrKMAzDiCWmoAzDMIxYYgrKMAzDiCWmoAzDMIxYYgrKMAzDiCX/D9K8xJ1ETQnKAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Read elev data\n", "nc_file = r'../../CRU_Climate_Data/elevation/elev.0.5-deg.nc'\n", "ds = xr.open_dataset(nc_file)\n", "\n", "# Shift negative in stn_df lons to match lons in elev file\n", "lon_shift = stn_df['longitude'].values.copy() # Copy lon values\n", "lon_shift[lon_shift<0] = lon_shift[lon_shift<0] + 360 # Shift\n", "stn_df['lon_shift'] = lon_shift # New column \n", "\n", "# Extract elevation data\n", "dsloc = ds.sel_points(lon=stn_df['lon_shift'].values,\n", " lat=stn_df['latitude'].values,\n", " method='nearest',\n", " dim=stn_df['station_id']) # Use stn ids to index points in result\n", "\n", "# Convert to df\n", "elev_df = dsloc['data'].to_dataframe()\n", "\n", "# Rename cols\n", "elev_df.columns = ['grid_lat', 'grid_lon', 'px_elev_m']\n", "\n", "# Reset index and tidy\n", "elev_df.reset_index(inplace=True)\n", "del elev_df['time']\n", "\n", "# Join to stn_df\n", "stn_df = pd.merge(stn_df, elev_df, how='left', on='station_id')\n", "\n", "# Get cols of interest\n", "stn_df = stn_df[['station_id', 'station_code', 'station_name', \n", " 'latitude', 'longitude', 'lon_shift',\n", " 'grid_lat', 'grid_lon', 'elevation', 'px_elev_m']]\n", "\n", "# Plot\n", "fig = plt.figure(figsize=(6,4))\n", "plt.plot(stn_df['elevation'], stn_df['px_elev_m'], 'ro')\n", "plt.xlabel('Measured site elevation (m)')\n", "plt.ylabel('Gridded pixel elevation (m)')\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Extract and aggregate climate data" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/opt/conda/lib/python3.6/site-packages/ipykernel_launcher.py:37: DeprecationWarning: Dataset.sel_points is deprecated: use Dataset.sel()instead.\n", "/opt/conda/lib/python3.6/site-packages/dask/array/numpy_compat.py:28: RuntimeWarning: invalid value encountered in true_divide\n", " x = np.divide(x1, x2, out)\n" ] } ], "source": [ "# Read all the netCDF files and combine into a single dataset\n", "nc_path = (r'../../CRU_Climate_Data/netcdf/v3.26/*.nc')\n", "ds = xr.open_mfdataset(nc_path,\n", " chunks={'lat': 90, 'lon': 90})\n", "\n", "# Excel file for output\n", "out_xls = r'../../update_autumn_2018/results/cru_climate_summaries.xlsx'\n", "\n", "# Create writer object\n", "writer = pd.ExcelWriter(out_xls)\n", "\n", "# Define variables and period of interest\n", "var_list = ['pre', 'tmp']\n", "per_list = ['ann', 'jja', 'jas']\n", "\n", "# Dict mapping vars and pers to freqs and stats\n", "attrs_dict = {('pre', 'ann'):['A', 'sum'],\n", " ('pre', 'jja'):['Q-FEB', 'sum'],\n", " ('pre', 'jas'):['Q-MAR', 'sum'],\n", " ('tmp', 'ann'):['A', 'mean'],\n", " ('tmp', 'jja'):['Q-FEB', 'mean'],\n", " ('tmp', 'jas'):['Q-MAR', 'mean']}\n", "\n", "# Loop over data\n", "for var in var_list:\n", " for per in per_list: \n", " # Resample to desired resolution\n", " if attrs_dict[(var, per)][1] == 'mean':\n", " ds2 = ds.resample(time=attrs_dict[(var, per)][0]).mean()\n", " else:\n", " ds2 = ds.resample(time=attrs_dict[(var, per)][0]).sum()\n", "\n", " # Extract point data\n", " ds2 = ds2.sel_points(lon=stn_df['longitude'].values,\n", " lat=stn_df['latitude'].values,\n", " method='nearest',\n", " dim=stn_df['station_id'])\n", "\n", " # Convert to df\n", " df = ds2[var].to_dataframe()\n", "\n", " # Add year col (and month col for seasonal)\n", " df['year'] = df.index.get_level_values('time').year \n", " if per != 'ann':\n", " df['month'] = df.index.get_level_values('time').month\n", " \n", " # Tidy\n", " df.reset_index(inplace=True)\n", "\n", " # If seasonal, get just season of interest\n", " if per == 'jja':\n", " # Totals for JJA are associated with month=8\n", " df = df.query('month==8')\n", " elif per == 'jas':\n", " # Totals for JAS are associated with month=9\n", " df = df.query('month==9')\n", "\n", " # Write to output\n", " df.to_excel(writer,\n", " '%s_%s' % (var, per),\n", " index=False)\n", "\n", "# Save and tidy\n", "writer.save()\n", "ds.close()\n", "ds2.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Climate trends" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Define variables and periods of interest\n", "var_list = ['pre', 'tmp']\n", "per_list = [[1990, 2016], [1990, 2004], [2002, 2016]]\n", "\n", "# Excel file of climate data\n", "clim_xls = r'../../update_autumn_2018/results/cru_climate_summaries.xlsx'\n", "\n", "# Output summary stats\n", "out_csv = r'../../update_autumn_2018/results/icpw_climate_stats.csv'\n", "\n", "# Output folder for plots\n", "plot_fold = r'../../update_autumn_2018/results/climate_plots'\n", "\n", "# Produce plots? (For testing)\n", "make_plots = False" ] }, { "cell_type": "code", "execution_count": 7, "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", "
station_idpre_ann_1990-2004_sigpre_ann_1990-2004_slppre_ann_1990-2016_avgpre_ann_1990-2016_sigpre_ann_1990-2016_slppre_ann_2002-2016_sigpre_ann_2002-2016_slppre_jas_1990-2004_sigpre_jas_1990-2004_slp...tmp_jas_1990-2016_slptmp_jas_2002-2016_sigtmp_jas_2002-2016_slptmp_jja_1990-2004_sigtmp_jja_1990-2004_slptmp_jja_1990-2016_avgtmp_jja_1990-2016_sigtmp_jja_1990-2016_slptmp_jja_2002-2016_sigtmp_jja_2002-2016_slp
0381030.092460-24.8636361277.0851850.045360-10.3000000.03766723.4500000.234955-5.3500...0.0645830.9211590.0047620.1245420.12963011.5147160.0333940.0564100.8818350.008333
1381040.047761-21.4800001298.7259260.013896-10.9809520.16585710.4857140.198211-6.1500...0.0615380.8818350.0095240.1507500.08333312.5187900.0269880.0611110.6548540.011111
2381050.047761-21.4800001298.7259260.013896-10.9809520.16585710.4857140.198211-6.1500...0.0615380.8818350.0095240.1507500.08333313.2027900.0269880.0611110.6548540.011111
3381060.092460-24.8636361277.0851850.045360-10.3000000.03766723.4500000.234955-5.3500...0.0645830.9211590.0047620.1245420.12963011.4667160.0333940.0564100.8818350.008333
4381070.766525-2.6076921178.8555560.867547-0.7083330.552615-6.1888890.428480-2.9875...0.0577780.691463-0.0476190.2981050.13333314.6780250.1082170.0476190.9604830.008333
\n", "

5 rows × 43 columns

\n", "
" ], "text/plain": [ " station_id pre_ann_1990-2004_sig pre_ann_1990-2004_slp \\\n", "0 38103 0.092460 -24.863636 \n", "1 38104 0.047761 -21.480000 \n", "2 38105 0.047761 -21.480000 \n", "3 38106 0.092460 -24.863636 \n", "4 38107 0.766525 -2.607692 \n", "\n", " pre_ann_1990-2016_avg pre_ann_1990-2016_sig pre_ann_1990-2016_slp \\\n", "0 1277.085185 0.045360 -10.300000 \n", "1 1298.725926 0.013896 -10.980952 \n", "2 1298.725926 0.013896 -10.980952 \n", "3 1277.085185 0.045360 -10.300000 \n", "4 1178.855556 0.867547 -0.708333 \n", "\n", " pre_ann_2002-2016_sig pre_ann_2002-2016_slp pre_jas_1990-2004_sig \\\n", "0 0.037667 23.450000 0.234955 \n", "1 0.165857 10.485714 0.198211 \n", "2 0.165857 10.485714 0.198211 \n", "3 0.037667 23.450000 0.234955 \n", "4 0.552615 -6.188889 0.428480 \n", "\n", " pre_jas_1990-2004_slp ... tmp_jas_1990-2016_slp \\\n", "0 -5.3500 ... 0.064583 \n", "1 -6.1500 ... 0.061538 \n", "2 -6.1500 ... 0.061538 \n", "3 -5.3500 ... 0.064583 \n", "4 -2.9875 ... 0.057778 \n", "\n", " tmp_jas_2002-2016_sig tmp_jas_2002-2016_slp tmp_jja_1990-2004_sig \\\n", "0 0.921159 0.004762 0.124542 \n", "1 0.881835 0.009524 0.150750 \n", "2 0.881835 0.009524 0.150750 \n", "3 0.921159 0.004762 0.124542 \n", "4 0.691463 -0.047619 0.298105 \n", "\n", " tmp_jja_1990-2004_slp tmp_jja_1990-2016_avg tmp_jja_1990-2016_sig \\\n", "0 0.129630 11.514716 0.033394 \n", "1 0.083333 12.518790 0.026988 \n", "2 0.083333 13.202790 0.026988 \n", "3 0.129630 11.466716 0.033394 \n", "4 0.133333 14.678025 0.108217 \n", "\n", " tmp_jja_1990-2016_slp tmp_jja_2002-2016_sig tmp_jja_2002-2016_slp \n", "0 0.056410 0.881835 0.008333 \n", "1 0.061111 0.654854 0.011111 \n", "2 0.061111 0.654854 0.011111 \n", "3 0.056410 0.881835 0.008333 \n", "4 0.047619 0.960483 0.008333 \n", "\n", "[5 rows x 43 columns]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get list of sites\n", "stn_list = stn_df['station_id'].unique()\n", "\n", "# Dict to store output\n", "data_dict = defaultdict(list)\n", "\n", "# Loop over data\n", "for var in var_list:\n", " for tm in ['ann', 'jja', 'jas']:\n", " # Open the climate data\n", " clim_df = pd.read_excel(clim_xls, sheet_name='%s_%s' % (var, tm))\n", " \n", " # Loop over stations\n", " for stn in stn_list:\n", " # Filter the climate data for this station\n", " stn_clim_df = clim_df.query('station_id == @stn')\n", " \n", " # Set date index\n", " stn_clim_df.index = stn_clim_df['time']\n", " del stn_clim_df['time']\n", " stn_clim_df = stn_clim_df.sort_index()\n", " \n", " # Correct temperatures according to lapse rate\n", " if var == 'tmp':\n", " # Get elevations\n", " stn_elev = stn_df.query('station_id == @stn')['elevation'].values[0]\n", " px_elev = stn_df.query('station_id == @stn')['px_elev_m'].values[0]\n", " \n", " # If pixel elev is negative (i.e. in sea), correct back to s.l.\n", " if px_elev < 0:\n", " px_elev = 0\n", " \n", " # Calculate temperature difference based on 0.6C/100m\n", " t_diff = 0.6 * (px_elev - stn_elev) / 100.\n", " \n", " # Apply correction\n", " stn_clim_df['tmp'] = stn_clim_df['tmp'] + t_diff\n", " \n", " # Loop over time periods\n", " for per in per_list:\n", " # Truncate\n", " df = stn_clim_df.truncate(before='%s-01-01' % per[0],\n", " after='%s-12-31' % per[1])\n", "\n", " # Only need averages for 1990-2016\n", " if (per[0]==1990) and (per[1]==2016):\n", " # Calculate long-term averages\n", " key = '%s_%s_%s-%s_avg' % (var, tm, per[0], per[1])\n", " val = df.mean()[var]\n", " data_dict[key].append(val)\n", " \n", " # Calculate Sen's slope and add to dict\n", " sslp, icpt, lb, ub = theilslopes(df[var].values, \n", " df['year'], 0.95)\n", " \n", " key = '%s_%s_%s-%s_slp' % (var, tm, per[0], per[1])\n", " data_dict[key].append(sslp)\n", " \n", " # Calculate MK signif and add to dict\n", " res = resa2_trends.mk_test(df[var].values, str(stn), var)\n", " sig = res[3]\n", " \n", " key = '%s_%s_%s-%s_sig' % (var, tm, per[0], per[1])\n", " data_dict[key].append(sig)\n", " \n", " # Plot\n", " if make_plots:\n", " plt.plot(df['year'], df[var].values, 'bo-')\n", " plt.plot(df['year'],\n", " sslp*df['year'] + icpt,\n", " 'k-')\n", " plt.title('%s %s at station %s (%s-%s)' % (tm, var, stn, per[0], per[1]),\n", " fontsize=20)\n", "\n", " # Save\n", " png_path = os.path.join(plot_fold, \n", " '%s_%s_%s_%s-%s.png' % (stn, tm, var, \n", " per[0], per[1]))\n", " plt.savefig(png_path, dpi=150)\n", " plt.close()\n", "\n", "# Build output df\n", "df = pd.DataFrame(data_dict, index=stn_list)\n", "\n", "# Reorder columns\n", "cols = df.columns\n", "cols = sorted(cols)\n", "df = df[cols]\n", "\n", "# Tidy\n", "df.index.name = 'station_id'\n", "df.reset_index(inplace=True)\n", "\n", "# Save\n", "df.to_csv(out_csv)\n", "\n", "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Join to chemistry data" ] }, { "cell_type": "code", "execution_count": 8, "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", "
project_idproject_namestation_idstation_codestation_namenfc_codetypecontinentcountryregion...tmp_jas_1990-2016_slptmp_jas_2002-2016_sigtmp_jas_2002-2016_slptmp_jja_1990-2004_sigtmp_jja_1990-2004_slptmp_jja_1990-2016_avgtmp_jja_1990-2016_sigtmp_jja_1990-2016_slptmp_jja_2002-2016_sigtmp_jja_2002-2016_slp
04390ICPW_TOCTRENDS_201838085Tr18_CA01Ontario, Algoma Region, Batchawana LakeLNACanadaONT...0.0750000.8427050.0083330.5857350.05238115.7941230.0992090.0466671.0000000.000000
14390ICPW_TOCTRENDS_201838086Tr18_CA02Ontario, Algoma Region, Wishart LakeLNACanadaONT...0.0750000.8427050.0083330.5857350.05238116.4481230.0992090.0466671.0000000.000000
24390ICPW_TOCTRENDS_201838087Tr18_CA03Ontario, Algoma Region, Little Turkey LakeLNACanadaONT...0.0750000.8427050.0083330.5857350.05238116.5261230.0992090.0466671.0000000.000000
34390ICPW_TOCTRENDS_201838088Tr18_CA04Ontario, Algoma Region, Turkey LakeLNACanadaONT...0.0750000.8427050.0083330.5857350.05238116.5441230.0992090.0466671.0000000.000000
44390ICPW_TOCTRENDS_201838089Tr18_CA05Quebec, Lac VeilleuxLNACanadaQuMaVt...0.0566670.1376460.0794870.298105-0.04166714.2848400.2105390.0190480.0198720.095833
\n", "

5 rows × 286 columns

\n", "
" ], "text/plain": [ " project_id project_name station_id station_code \\\n", "0 4390 ICPW_TOCTRENDS_2018 38085 Tr18_CA01 \n", "1 4390 ICPW_TOCTRENDS_2018 38086 Tr18_CA02 \n", "2 4390 ICPW_TOCTRENDS_2018 38087 Tr18_CA03 \n", "3 4390 ICPW_TOCTRENDS_2018 38088 Tr18_CA04 \n", "4 4390 ICPW_TOCTRENDS_2018 38089 Tr18_CA05 \n", "\n", " station_name nfc_code type continent country \\\n", "0 Ontario, Algoma Region, Batchawana Lake L NA Canada \n", "1 Ontario, Algoma Region, Wishart Lake L NA Canada \n", "2 Ontario, Algoma Region, Little Turkey Lake L NA Canada \n", "3 Ontario, Algoma Region, Turkey Lake L NA Canada \n", "4 Quebec, Lac Veilleux L NA Canada \n", "\n", " region ... tmp_jas_1990-2016_slp tmp_jas_2002-2016_sig \\\n", "0 ONT ... 0.075000 0.842705 \n", "1 ONT ... 0.075000 0.842705 \n", "2 ONT ... 0.075000 0.842705 \n", "3 ONT ... 0.075000 0.842705 \n", "4 QuMaVt ... 0.056667 0.137646 \n", "\n", " tmp_jas_2002-2016_slp tmp_jja_1990-2004_sig tmp_jja_1990-2004_slp \\\n", "0 0.008333 0.585735 0.052381 \n", "1 0.008333 0.585735 0.052381 \n", "2 0.008333 0.585735 0.052381 \n", "3 0.008333 0.585735 0.052381 \n", "4 0.079487 0.298105 -0.041667 \n", "\n", " tmp_jja_1990-2016_avg tmp_jja_1990-2016_sig tmp_jja_1990-2016_slp \\\n", "0 15.794123 0.099209 0.046667 \n", "1 16.448123 0.099209 0.046667 \n", "2 16.526123 0.099209 0.046667 \n", "3 16.544123 0.099209 0.046667 \n", "4 14.284840 0.210539 0.019048 \n", "\n", " tmp_jja_2002-2016_sig tmp_jja_2002-2016_slp \n", "0 1.000000 0.000000 \n", "1 1.000000 0.000000 \n", "2 1.000000 0.000000 \n", "3 1.000000 0.000000 \n", "4 0.019872 0.095833 \n", "\n", "[5 rows x 286 columns]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Read chem data\n", "csv_path = r'../../update_autumn_2018/results/toc_trends_wide_format.csv'\n", "wide_df = pd.read_csv(csv_path, keep_default_na=False)\n", "\n", "# Join climate to trends\n", "clim_df = pd.merge(wide_df, df, how='left', on='station_id')\n", "\n", "# Save output\n", "out_path = r'../../update_autumn_2018/results/toc_trends_wide_format_with_climate.csv'\n", "clim_df.to_csv(out_path, index=False, encoding='utf-8')\n", "\n", "clim_df.head()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }