{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# inter-run comparisons of N contained in phytoplankton\n",
" - S3\n",
" - Sentry Shoal\n",
" - Central Node\n",
" - JDF"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import netCDF4 as nc\n",
"import datetime as dt\n",
"import subprocess\n",
"import requests\n",
"import matplotlib.pyplot as plt\n",
"import cmocean\n",
"import numpy as np\n",
"import os\n",
"import re\n",
"import dateutil as dutil\n",
"from salishsea_tools import viz_tools, places\n",
"import glob\n",
"import pickle\n",
"import matplotlib.dates as mdates\n",
"import matplotlib as mpl\n",
"mpl.rc('xtick', labelsize=14) \n",
"mpl.rc('ytick', labelsize=16) \n",
"mpl.rc('legend', fontsize=16) \n",
"mpl.rc('axes', titlesize=16) \n",
"mpl.rc('figure', titlesize=16) \n",
"mpl.rc('axes', labelsize=16) \n",
"mpl.rc('font', size=16) \n",
"\n",
"mpl.rcParams['font.size'] = 16\n",
"mpl.rcParams['axes.titlesize'] = 16\n",
"mpl.rcParams['legend.numpoints'] = 1\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"plist=['Sentry Shoal','S3','Central node','Central SJDF']"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/eolson/anaconda3/envs/py38/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3050: DtypeWarning: Columns (3) have mixed types. Specify dtype option on import or set low_memory=False.\n",
" has_raised = await self.run_ast_nodes(code_ast.body, cell_name,\n"
]
}
],
"source": [
"df1=pd.read_csv('/ocean/eolson/MEOPAR/obs/ONC/turbidity/nearSurface/search3928586/BritishColumbiaFerries_Tsawwassen-DukePoint_Turbidity-ChlorophyllandFluorescence_20140804T234330Z_20150604T070614Z-clean.csv',\n",
" skiprows=78,header=None,\n",
" names=('TimeUTC','CDOM','CDOMQC','Chlorophyll_ug','ChlQC','Turbidity_NTU','TurbQC','Lat','LatQC','Lon','LongQC'))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/eolson/anaconda3/envs/py38/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3050: DtypeWarning: Columns (1,3,5) have mixed types. Specify dtype option on import or set low_memory=False.\n",
" has_raised = await self.run_ast_nodes(code_ast.body, cell_name,\n"
]
}
],
"source": [
"df2=pd.read_csv('/ocean/eolson/MEOPAR/obs/ONC/turbidity/nearSurface/search3928586/BritishColumbiaFerries_Tsawwassen-DukePoint_Turbidity-ChlorophyllandFluorescence_20150604T070624Z_20160307T160206Z-clean.csv',\n",
" skiprows=78,header=None,\n",
" names=('TimeUTC','CDOM','CDOMQC','Chlorophyll_ug','ChlQC','Turbidity_NTU','TurbQC','Lat','LatQC','Lon','LongQC'))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"df=pd.concat([df1.drop(df1[df1.TimeUTC<'2015'].index),df2.drop(df2[df2.TimeUTC>'2016'].index)],ignore_index=True)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"dts=[dt.datetime(int(r[0:4]),int(r[5:7]),int(r[8:10]),int(r[11:13]),int(r[14:16]),int(r[17:19])) for r in df['TimeUTC']]"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"df=df.assign(dts=dts)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"df['Lat']=pd.to_numeric(df['Lat'],errors='coerce')\n",
"df['Lon']=pd.to_numeric(df['Lon'],errors='coerce')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>TimeUTC</th>\n",
" <th>CDOM</th>\n",
" <th>CDOMQC</th>\n",
" <th>Chlorophyll_ug</th>\n",
" <th>ChlQC</th>\n",
" <th>Turbidity_NTU</th>\n",
" <th>TurbQC</th>\n",
" <th>Lat</th>\n",
" <th>LatQC</th>\n",
" <th>Lon</th>\n",
" <th>LongQC</th>\n",
" <th>dts</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>2015-01-01T00:00:08.461Z</td>\n",
" <td>16.5912</td>\n",
" <td>1</td>\n",
" <td>1.2463</td>\n",
" <td>1</td>\n",
" <td>4.5548</td>\n",
" <td>1</td>\n",
" <td>49.094471</td>\n",
" <td>8</td>\n",
" <td>-123.426450</td>\n",
" <td>8</td>\n",
" <td>2015-01-01 00:00:08</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>2015-01-01T00:00:17.515Z</td>\n",
" <td>16.948</td>\n",
" <td>1</td>\n",
" <td>1.2342</td>\n",
" <td>1</td>\n",
" <td>4.7864</td>\n",
" <td>1</td>\n",
" <td>49.094835</td>\n",
" <td>8</td>\n",
" <td>-123.427538</td>\n",
" <td>8</td>\n",
" <td>2015-01-01 00:00:17</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2015-01-01T00:00:27.689Z</td>\n",
" <td>16.7696</td>\n",
" <td>1</td>\n",
" <td>1.2221</td>\n",
" <td>1</td>\n",
" <td>4.632</td>\n",
" <td>1</td>\n",
" <td>49.095233</td>\n",
" <td>8</td>\n",
" <td>-123.428714</td>\n",
" <td>8</td>\n",
" <td>2015-01-01 00:00:27</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>2015-01-01T00:00:37.874Z</td>\n",
" <td>16.7696</td>\n",
" <td>1</td>\n",
" <td>1.3673</td>\n",
" <td>1</td>\n",
" <td>4.632</td>\n",
" <td>1</td>\n",
" <td>49.095742</td>\n",
" <td>8</td>\n",
" <td>-123.429833</td>\n",
" <td>8</td>\n",
" <td>2015-01-01 00:00:37</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>2015-01-01T00:00:48.048Z</td>\n",
" <td>16.948</td>\n",
" <td>1</td>\n",
" <td>1.2826</td>\n",
" <td>1</td>\n",
" <td>4.632</td>\n",
" <td>1</td>\n",
" <td>49.096250</td>\n",
" <td>8</td>\n",
" <td>-123.430950</td>\n",
" <td>8</td>\n",
" <td>2015-01-01 00:00:48</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" TimeUTC CDOM CDOMQC Chlorophyll_ug ChlQC \\\n",
"0 2015-01-01T00:00:08.461Z 16.5912 1 1.2463 1 \n",
"1 2015-01-01T00:00:17.515Z 16.948 1 1.2342 1 \n",
"2 2015-01-01T00:00:27.689Z 16.7696 1 1.2221 1 \n",
"3 2015-01-01T00:00:37.874Z 16.7696 1 1.3673 1 \n",
"4 2015-01-01T00:00:48.048Z 16.948 1 1.2826 1 \n",
"\n",
" Turbidity_NTU TurbQC Lat LatQC Lon LongQC \\\n",
"0 4.5548 1 49.094471 8 -123.426450 8 \n",
"1 4.7864 1 49.094835 8 -123.427538 8 \n",
"2 4.632 1 49.095233 8 -123.428714 8 \n",
"3 4.632 1 49.095742 8 -123.429833 8 \n",
"4 4.632 1 49.096250 8 -123.430950 8 \n",
"\n",
" dts \n",
"0 2015-01-01 00:00:08 \n",
"1 2015-01-01 00:00:17 \n",
"2 2015-01-01 00:00:27 \n",
"3 2015-01-01 00:00:37 \n",
"4 2015-01-01 00:00:48 "
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df.head()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-123.42645043"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['Lon'][0]"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'lon lat': (-123.558, 49.125),\n",
" 'NEMO grid ji': (450, 258),\n",
" 'GEM2.5 grid ji': (138, 144)}"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"places.PLACES['S3']"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.5086"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"111*.0226"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"llon=places.PLACES['S3']['lon lat'][0]-.01\n",
"ulon=places.PLACES['S3']['lon lat'][0]+.01\n",
"llat=places.PLACES['S3']['lon lat'][1]-.01\n",
"ulat=places.PLACES['S3']['lon lat'][1]+.01\n",
"iidfnd=(df.Lon>llon)&(df.Lon<ulon)&(df.Lat>llat)&(df.Lat<ulat)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/eolson/anaconda3/envs/py38/lib/python3.8/site-packages/pandas/plotting/_matplotlib/converter.py:103: FutureWarning: Using an implicitly registered datetime converter for a matplotlib plotting method. The converter was registered by pandas on import. Future versions of pandas will require you to explicitly register matplotlib converters.\n",
"\n",
"To register the converters:\n",
"\t>>> from pandas.plotting import register_matplotlib_converters\n",
"\t>>> register_matplotlib_converters()\n",
" warnings.warn(msg, FutureWarning)\n"
]
},
{
"data": {
"text/plain": [
"(datetime.datetime(2015, 3, 6, 0, 0), datetime.datetime(2015, 3, 16, 0, 0))"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABLsAAAD8CAYAAACfFKY8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdf3Cc5X33+8+1u1rtroxsMcXxKgJbro1TK7WlRMFx+wRTCDak/LIBYyQSwzMPwtacmpIhHNGZZHrsM0koeVKmGVugdDqBp4YopqIkk/YZO4ZEMQ0FE0nEogfiSDIYr33SmmOD9cuSrvOHvYtWuven9qf0fs14kL/3dd/3JbO72v3qe30vY60VAAAAAAAAMBu48j0BAAAAAAAAIFNIdgEAAAAAAGDWINkFAAAAAACAWYNkFwAAAAAAAGYNkl0AAAAAAACYNTz5nsBM/cEf/IFdsmRJvqcBAAAAAAAwa7zxxhv/aa29LN/zSEfRJ7uWLFmiw4cP53saAAAAAAAAs4Yx5li+55AuljECAAAAAABg1iDZBQAAAAAAgFmDZBcAAAAAAABmDZJdAAAAAAAAmDVIdgEAAAAAAGDWINkFAACAnAqFQlq3bp1OnjyZ76kAAIBZiGQXAAAAcqqlpUWdnZ1qaWnJ91QAAMAsZKy1+Z7DjNTX19vDhw/nexoAAABIwO/3a3h4eFrc5/NpaGgoDzMCAACxGGPesNbW53se6aCyCwAAADkR65esxf7LVwAAUFhIdgEAAAAAAGDWINkFAAAAAACAWYNkFwAAAHKiv79fZWVlUbGysjINDAzkZ0IAAGBWItkFAACAnAgGg7rsssuiYgsXLtSiRYvyNCMAADAbefI9AQAAAMwNTrsx9vf3y+/3sxsjAADIGCq7AAAAkBOp7sYYCoW0bt06nTx5MpvTAgAAswzJLgAAAOREf3+/li1bFhVbvnx5zJ5du3bt0qFDh7Rz584czA4AAMwWJLsAAACQE8FgUGNjY5Ikr9crSRobG5vWs8vv98sYo9bWVk1MTKi1tVXGGPn9/pzPGQAAFB+SXQAAAMiZuro6bd26VZ/+9Kd17733qra2dtqYvr4+NTQ0RJJbfr9fjY2N6u/vz/V0AQBAESLZBQAAgJzp6OhQIBBQd3e3/H6/Ojo6po0JBoMqLy+PNK0fGhpSeXk5uzYCAICkmFgNQYtFfX29PXz4cL6nAQAAgAScdmOUJJ/PF7UbY7LjAABA9hhj3rDW1ud7HumgsgsAAAA5EV6e6PP5JF1IXjktTwyPCwQCkqRAIMAyRgAAkDSSXQAAAMiJ8PLEcNXW8PCw4/LE8LjBwUFJ0uDgIMsYAQBA0kh2AQAAICeMMXryySejYuGdFqc6depUZMdGr9erkydP5mSOAACg+HnyPQEAAADMDV6vV6Ojo47xyaYmv0ZHR/XCCy/IGKNi7zcLAACyj8ouAAAA5MTAwIDKy8ujYvPnz9exY8eiYvv374/06worKyvTwYMHsz5HAABQ/NiNEQAAADlTVlamwcHBSJVWIBDQuXPnosawGyMAAPnHbowAAABAEjZs2KCamhpJUk1NjTZs2DBtTF9fn1yu6LepLpcr7d0YQ6GQ1q1bR98vAADmCJJdAAAAyAm/368XXnhBvb29staqt7dXL7zwgvx+f9S4YDCopqamqNgDDzyQ8m6M4STXo48+qkOHDmnnzp0z/h4AAEDhYxkjAAAAciIUCunhhx/WCy+8oKGhIfn9fm3atEnf+c53ohJZTrszhqXy3tXj8Wh8fHxanOWQAAAkxjJGAAAAIIFgMKjy8vJIomloaEjl5eUpV2wl4vf7ZYxxTHQ1NjamvRwSAAAUB5JdAAAAyAm/368nn3wyKtba2jptGWNXV5fj+T09PUndp6+vTw0NDY7H9u7dq+rq6qSuAwAAihPJLgAAAOREOAkVCAQkSYFAwLHSqq6uzvH81atXR76O13R+6dKlevbZZx2vUVVVRWUXAACzHMkuAAAA5MzLL7+soaEh+Xw+DQ8Pq7y8XNbaqMSV1+t1PHdy/Ctf+Yo6Ozu1devWaeOGh4dj3v/48eNUdgEAMMuR7AIAAEBOtLS0KBQK6corr9Srr76qbdu26eTJk9q1a1fUbomjo6OO54+Ojkb6cf3sZz+TJO3fv1/GGPl8vsi4rq4uzZs3L+Y8pibD4lWJAQCA4sNujAAAAMgqv98ft9pqqpKSEp0/f35a/ODBg7ruuutinjf5fW15ebk+/PBDx3Fer1cjIyORv2/dulXPPPOMtm7dqh/84AdJzxMAgNmM3RgBAACAGGL9ctXj8aihoSHSoN7v96uxsVHvvvuuPB5P1NiSkhJde+21KikpiXkfY0zkWvEqu8KVY+EqsWeeeUaS9PTTT0ddAwAAFCeSXQAAAMiLsbExlZeXa2hoSJI0NDSk8vJyLVq0SG63WxUVFXrsscdUUVEhl+vC29Zjx47Fvebw8LCMMQqFQgnvHysJV+wrHwAAmOuytozRGLNB0v8paaWkCkm/l/Rvkv7aWvvWpHEVkh6XdJskv6RfSXrIWvubZO7DMkYAAIDCZoxJabzP54skwKaKtyRy3rx5euWVV/S5z30uZt8v6UKV2OjoqHw+X9RyxrDS0tKUll0CADAbsYzR2aWS3pD0f0haL+lRSTWSXjXGLJYkc+Gdz48l3SDpLyTdLqlE0svGmKoszg0AAAAFIBAIRJJhxhg1Njaqv78/5vh4SShjjFatWqWBgQGVl5fHHBfuB+aU6IoXBwAAxSFryS5r7XPW2q9Za5+31v7CWvu/JG2SdImkOy4Ou0XSf5P05Yvj//fFmEvSI9maGwAAAHLnueeei3lscHAwsmzQWqu9e/cqGAymdZ8PP/xQxhhVVlbq7NmzMcf19PRIurBro9vtjjrmdrsjxwEAQHHKdc+u/7r43/D2OrdIOmGtfTk8wFp7RtJPJN2a47kBAAAgC+67776kx/r9/hklmzwej7xeb9wxa9askSStXbtW4+PjUcfGx8cjxwEAQHHKerLLGOM2xniNMcslPSXppKQfXjxcI+mIw2m9kq4wxsTeRgcAAABFId7Sw5UrV0b9fenSpVq1alXM8YkSWW63O26/Lp/PF1kmSYN6AABmp1xUdv27pBFJ70haJelaa+3/e/HYpZI+cDjn9MX/Vjhd0BjTZIw5bIw5/Pvf/z7T8wUAAEAGxUtQvfXWW1HjTp8+HXOspLiJLEkaGBhQSUlJzOPDw8ORZZL9/f1atmxZ1PHq6moNDAzEvQcAAChsuUh2fVnS5yU1SDor6YAxZsnFY0aS06/O4m7ZY61ts9bWW2vrL7vssgxOFQAAAJmW7G6Mo6OjCoVCKe/eOFl1dbVcrthvcd1utw4ePChJqqys1NGjR6OO9/f3p90zDAAAFIasJ7ustf9hrf13a+1zkq6TNE9Sy8XDp3WhumuqcEWXU9UXAAAAZrF4lWD79++Pe+7w8HDc3RTLy8t17bXXJryP3+9PMEsAAFCoctqg3lr7/0k6KilcL96rC327plop6V1r7Ue5mhsAAACyo7+/X2VlZUmPj1fZdf3118c91+v1xk1iDQ4OJnWfeH3GAABAYctpsssY8wlJn5L0u4uhH0v6pDFm3aQx5ZJuvngMAAAARa6yslLnzp3L2PVKS0tjHjPGaGBgQB6PJyru8XgUCoWikljxKsASNcIHAACFy5N4SHqMMS9I+rWkN3WhV9eVkh6SNCbpf14c9mNJv5L0j8aYr+nCssVHdaFn199ka24AAAAoXsPDwzGrskZGRhQMBuX1ejU2NqaSkhKdP39eXq9XixYtSvoex44dy9R0AQBAjmWzsutVSbdJelrSTyV9VdIvJNVaa9+RJGvthKSbJB2QtEfSC5LGJf2Ztfa9LM4NAAAAOdLV1TWt0kqSgsGg3G53VKysrCyp3RBjVV6F4xs2bFBzc7Nef/11NTc3a8OGDdPGxuv/VV1dnXAOAACgMBlrnTZDLB719fX28OHD+Z4GAAAA4vB4PBofH09qrM/n09DQUMIxsZYh9vT0aNWqVUndK17frmJ/nwwAwEwYY96w1tbnex7pyGnPLgAAAMxNCxcuTHpsMkmm/v7+aVVhYQ0NDUndJ96Oi/TsAgCgeJHsAgAAQNadOHFCXV1dmjdvXlR8amXV8uXLk1rGuHTp0piVYr29vTLGxK3akqS+vr6YYxKdCwAAChfJLgAAAOREbW1twjFjY2NJNZLv6+tTVVVVzOquJUuWqKenJ+41gsGg7rnnnmlxt9udVMINAAAUJpJdAAAAyJlLLrlENTU1am9vlzR9yWJ/f39SVVXBYFA33XSTrLXy+XzTjpeVlSXVt2vv3r3TYuPj4zSoBwCgiJHsAgAAQM6cOHFCR44c0ebNmx13aJQUMz7VqVOntG3bNg0PD087Fl7KmMjx48dVVVUVFauqqlJ/f39ScwAAAIWH3RgBAACQF6FQSJ/85CejqrtcLpfef//9pJYyhs1kR0W/3++YLEtmR0gAAGYzdmMEAAAAUhQMBqcloyYmJlJKdElSV1eXYzxRzy4pdjKs2H8hDADAXEayCwAAAHnh9/tTiseydu1ax/iaNWsSntvf369ly5ZFxZLdERIAABQmkl0AAADIi76+PjU0NESSW36/X42NjSn3y5pJdVYwGNTY2Jikj5dDJrsjZDpCoZDWrVunnp4erVu3TidPnszKfQAAmMtIdgEAACAvgsGgysvLI72xhoaGVF5enrVEUyx1dXVasWKFrLVasWKFamtrs3avXbt26dChQ2psbNShQ4e0c+fOrN0LAIC5igb1AAAAyItMNYcPhUJavny5zp07F4mVlZXp6NGjCRNnuWpQH+s+2bofAAAzRYN6AAAAIEXhZYyBQECSFAgE0lrGGAwGddlll0mSvF6vJGnhwoVJVYjlqkH91CWbk6XzPQMAgNhIdgEAACAvwssYh4eH5fP5NDw8nPYyxrq6Ot15553yer3avHlz0ksRc9WgfuqSzcn27t2rYDCY0fsBADCXefI9AQAAAMxdp06d0rZt29TU1KS2tjaFQqG0rtPR0aEVK1boo48+Und3t95+++2kzpvaoN5am7UG9adOnYp5LFyRBgAAZo5kFwAAAPKmo6Mj8vXu3bvTukZ4F8Wwd955JxJLZjliXV2dysrK9NZbb6mmpkZXXnllWvNIZY5TjY6OZvyeAADMVSS7AAAAUNS8Xq9jsiiZaqmpjeN7e3vV29srv9+f04bxPT09ObsXAACzHT27AAAAUNRiVU0lqqaSMtckP5GSkpK4x9esWZPR+wEAMJeR7AIAAMCcNblxvMvl0tDQUNpN8uNxueK/7Z5cXQYAAGaGZBcAAADmtFOnTmnlypWy1mrlypU6efJkxu8Rr1Js3rx5LGMEACCD6NkFAACAghcKhbRlyxa1t7dntOoqVz27KisrYx776KOPtHr16qSa6QMAgMSo7AIAAEDB27Vrlw4dOqSdO3dOO9bf36+ysrKoWFlZmQYGBhJeN9yzy+12S5LcbndWenYlapafqKcXAABInin23yDV19fbw4cP53saAAAAyIKplVdhPp8vqvKquro6KrlVXV2tvr6+jF1/pnw+n0ZGRhyPlZSU6N133814nzAAAGbCGPOGtbY+3/NIB5VdAAAAKFjJVl7V1dWppqZGxhhdeeWVOnv2bFK9t2L94jfTvxCOleiSpPHxcRJdAABkEJVdAAAAKFjJVF7FGuNyuTQ+Ph73+qFQSFdffbWOHj0aiS1fvlydnZ0ZTUAZYxKOKfb35QCA2YXKLgAAACALkqm8Cld/TTUxMSFjjPx+f8zrB4NBjY2NSfq4r9bY2FjGK626urpiJryWLFnCbowAAGQQuzECAACgqAWDQZWXl8c8nqhiqq6uTl/60pfU1NSktrY2hUKhTE9RdXV1MY8NDAxozZo1Ge0RBgDAXEZlFwAAAIreqVOntHXrVl1++eVR8eXLl0/blTEUCmndunWRnl4dHR3avXu3Vq9erd27d6ujoyPj8+vq6pLL5fzW2+fzZXz3RwAA5jKSXQAAAChYsRq7T413dHToBz/4QaSRfbwliS0tLers7FRLS0sWZuystrZWn/rUp6bFjTG67777aFAPAEAGkewCAABAwfJ4nLtuxIrX1dWpublZr732mpqbm1VbWxs55vf7ZYzRM888I0l6+umnE/b0yqQPPvhANTU1+tznPhdJxq1cuTKpXSMBAEDy2I0RAAAABcvn8zlWd5WWljruwChdWKa4ZcsWtbe3R1VMpXOtTEtmd0kAAAoBuzECAAAAWZDsMsbJduzYoc7OTj344IMzvlamhXeODAQCkqRAIKDGxkZ6dgEAkEEkuwAAAFCwurq6In24wtxut3p6eqaNDS9TfP755yVJP/rRj3K6TDEZ4Z0jh4eH5fP5NDw8rPLycnp2AQCQQSS7AAAAULDq6uo0Pj4eFRsfH9fq1aunjY21FDEcD/fJmipWPFtOnTqlbdu26dVXX9W2bdvo2QUAQIY5d/YEAAAAikxXV5e+8IUv6KOPPorE5s2bp1deeUXShZ0PncSKZ0tHR0fk6927d+f03gAAzAVUdgEAAKBgdXV1Tdt5saSkxHEZ4+SdF8OMMVq1alXW5gcAAAoPyS4AAAAUrLVr12psbCwqdv78ea1Zs2baWL/fH1XVJUkffvhhpGdXf3+/ysrKoo6XlZVpYGAgs5MGAAB5xTJGAAAAFKxEfbhSGVtZWTnt2Llz5xQMBmWtncEsAQBAIaGyCwAAAAWrq6tL8+bNi4rNmzfPcRljogb0hdKgHgAAZBfJLgAAABSstWvXTlua+NFHHzkuYxwdHXW8Rjie6DgAAJgdSHYBAACgYPX19cnn80V2TDTGyOfzqb+/f9rY/fv3O17j4MGDki5UiTlxqhIDAADFi2QXAAAAClYwGNS9994bSXIZY3Tfffdp0aJF08auX7/e8RrXXXedJKmurs7x+OrVqzM3YQAAkHckuwAAAJBVoVBI69at08mTJ9M6/9SpU9q2bZteffVVbdu2LeZ1PB7nvZfC8a6uLrlc0W9/XS4XlV0AAMwyJhs7zxhj7pB0t6R6SQslvSupQ9I3rbUfThpXIelxSbdJ8kv6laSHrLW/SfZe9fX19vDhwxmcPQAAADKpublZTz31lB544AHt2bMna/fx+XwaGRmZFi8tLdXw8LD8fr/jjo0+n09DQ0NZmxcAAMXIGPOGtbY+3/NIR7aSXa/qQoLrRUnHJdVJ+mtJ/4+kP7HWTpgLjRc6JVVL+pqkDyQ9KqlGUq219ngy9yLZBQAAUJhynVxKlOwKhUJaunSpRkZGZK2VMUalpaXq7+93XBYJAMBcVszJrmwtY7zZWrvZWrvXWvsLa+0TknZIWiPpmotjbpH03yR92Vr7nLX2f1+MuSQ9kqV5AQAAIEf6+vrU0NAgt9stSXK73WpsbHRsLp8J/f39CgQCUbGysjINDAxISq3/FwAAKF5ZSXZZa3/vEH794n8/efG/t0g6Ya19edJ5ZyT9RNKt2ZgXAAAAcmfp0qV69tlnNT4+LkkaHx/X3r17VV1dnZX7BYNBLVy4UJLk9XolSQsXLoxKZiXb/wsAABQv5y6e2bHu4n//4+J/ayQdcRjXK+krxph51tqPcjIzAAAAZFxfX5+uuuoqHT/+cXeKqqoqvf7663HOmpm6ujp96UtfUlNTk9ra2hQKhaKOd3R0RL7evXt31uYBAADyJys9u6bdxJhPSuqS1GOtvf5i7B1Jv7bWbpky9n9I+r6kK6y178W4XpOkJkm64oorPnvs2LFsTh8AAABpoCE8AADFi55dcRhj5ulCo/oxSfdNPiTJKdNmEl3TWttmra231tZfdtllmZkoAAAAMsop0RUvnmuhUEjr1q1jKSMAALNMVpNdxhifpB9LWippw5QdFk9LutThtIqL//0gm3MDAABAdoX7ZiUbz7WWlhZ1dnaqpaUl31MBAAAZlLVklzGmRNI/SbpK0pestb+ZMqRXF/p2TbVS0rv06wIAAChuo6OjKcVzxe/3yxijZ555RpL09NNPyxgjv9+f13kBAIDMyEqyyxjjkrRX0nWSbrXWvuow7MeSPmmMWTfpvHJJN188BgAAAGRcrJ61uehlCwAAsi9blV27Jd0p6TuSzhljPj/pT9XFMT+W9CtJ/2iM2WKM2XAxZiT9TZbmBQAAgBx57rnnHOP79u3L8UwAAMBckpXdGI0xA5IWxzj8f1lr//riuEt1ISF2mySfLiS/vmqt7Un2XvX19fbw4cMzmi8AAAAyz5jY+w7ls4oqFApp+fLlOnfuXCRWVlamo0ePatGiRXmbFwAAhYTdGKew1i6x1poYf/560rjT1tr/bq291FobsNZel0qiCwAAAIWrq6vLMd7Tk9+3e5WVlVGJLkk6d+6cgsFgVu/L7o8AAORGVndjBAAAwNxVV1fnGF+9enWOZxItX7tEsvsjAAC5QbILAAAAWdHV1SWXK/rtpsvlyntlV6zllfGWXc4Euz8CAJBbJLsAAACQFWvXrtXExERUbGJiQmvWrMnTjGYunaWI7P4IAEBukewCAABAVvT19cnn80Uqpowx8vl86u/vz+u8RkZGUopP9uCDD6qzs1M7duzIyf0AAEDqSHYBAAAgK4LBoO69995IkssYo/vuuy/vOx6WlJQkjHd3d2vBggV68803JUmlpaUyxmjfvn2SpH379rEUEQCAAkWyCwAAAFlz6tQpbdu2Ta+++qq2bdtWEDsRnj9/PmH8rrvu0pkzZ7R582ZJ0ujoqOM5w8PDCe+Xr4b4AADMVabYewXU19fbw4cP53saAAAAKBKZakQ/b948vfLKK1q1alXccaFQSJWVlY7xfFe5AQAQizHmDWttfb7nkQ4quwAAADCnxKu0SqXa6qOPPtLq1asTjgsGg9OWOwYCARJdAABkCckuAAAAzCkDAwPTqruMMTp27FjM5YqxJJMcM8ZoaGgoKjY4OJixCjMAs1s6u8ACcx3JLgAAAMwpwWBQPp9PkuTxeCRJfr9fixYt0v79+1O6VjLJsa6uLi1evDgqtmTJEvX09KR0LwBzU0tLizo7O9XS0pLvqQBFg2QXAAAA5pwbbrhBzc3NOnz4sJqbm7VhwwZJ0i233JLxe9XW1qqsrCwqVlZWlrDXF4C5ze/3yxijZ555RpL09NNPJ9wFliow4AIa1AMAAKAgdHd365prrlFnZ2feEkGhUEhVVVWamJhIarzH44m5u2OY3+933LXR5/NNW94IAGE+n08jIyPT4qWlpTF3gm1ubtZTTz2lBx54QHv27Mn2FDHL0aAeAAAAmKF77rlHZ86cUUNDQ97mEAwG1dTUlPT49957L+GYWL9cLvZfOgPILqdEV6x4uAqstbVVExMTam1tTVgFBsxmJLsAAACQV8YYGWPU29srSert7Y3E8jGXJ598Munxyeyo2N/fr2XLlkXFli9froGBgVSnBwCO+vr6pv2ioLGxUf39/XmaEZBfJLsAAACQV4XUwD2Z3RUlye12a+PGjUmNDQaDGhsbi4qNjY0llSgDMHeVlJQkHa+srNSzzz4bFdu7d6+CwWBW5gYUOnp2AQAAYEZCoZC2bNmi9vb2tBI48Sq4cv1eNdlqsng9c6aiZxeAdKTSs6uQXkcxe9CzCwAAAHNWS0uLOjs71dLSku+pzFhXV1fCMakuQaRnF4BixM6OKGZUdgEAACAtmapYKqSKhFjfk5Nk51ZI3x+A4pHKa4fL5XJ8PTHGJL277FR33nmnnn/+eW3evFnt7e1pXQPFjcouAAAAzDmzsWIp2URXsr29pAvVYm63Oyrmdrvz0pMMwOz0/vvvO8ZPnDiR8rXCOzs+//zzkqQf/ehH7OyIokOyCwAAAHk1NRGUKJ5NySaxUtkpcu3atRofH4+KjY+Pa82aNSnNDQBSlU4fxVhJ/2R/GQAUApJdAAAASItT4+R48Vjee+89x/jx48dTntNMJZPEcrvdKfXs6uvrk8/ni1zbGCOfz6f+/v50pwlgDojVQ9CpKnTXrl2SpJqaGnV3d6u5uTnpHWOnipX0T6WiFcg3kl0AAABIS6Y+EC1dutQxXl1dnfKcZqq/vz9mRVn4+7riiitSqpYIBoO69957I0kuY4zuu+++tCouAMwdtbW1jkugV61aFfl7eMlha2urJKm3t1e1tbX6h3/4B3V0dKR131hJ/8lxmtej0JHsAgAAQFoGBgbk8XiiYh6PR8eOHUvpOn19faqqqop8qHO73aqqqspL5dPSpUunLTkMz+m1115Tc3OzamtrU77uqVOntG3bNr366qvatm0bHxABJGXhwoWqqalRe3u7ampqtHDhwqjjfX19amhoiCTjvV6vGhsbZ/T6mUzV7q5du3To0CHt3Lkz7fsA2cRujAAAAEhbWVmZBgcHI38PBAI6d+5cytfZvn272tra5PV6NTo6qgceeEB79uzJ5FSTEgqFdNVVVykUCml8fFxut1vBYFCvv/46lVgACtL27dv15JNPRv19Jq+f8ZZz+3y+jOzCi+LAbowAAACYc4wxUYkuSRocHEypeXtYoVQ+BYNB3XTTTbLWyufzyVqrm2++mUQXgIJkjIlKdElSa2trWq/DyQhXkk0200oyIBuo7AIAAEBa4n2YKub3mJs2bVIwGFRTU5Pa2toUCoXS7n1TKA4cOKAbb7xR+/fv17XXXpvv6QDIkGy8DqebKCvm1304K+bKLpJdAAAASMtsTXbNRuXl5frwww91ySWX6OzZs/meDoAMIdmFbCrmZBfLGAEAAJCWrq4ux3hPT0+OZ4JYjDEyxujDDz+UJH344YeRGAA4iffaXlJS4ngsVhzIF5JdAAAASEtdXZ1jfPXq1TmeCQDMTY8++qhj/Otf/3ra16ytrY3sjhvm8Xi0atUquVzOKYRYcSBfeEQCAAAAaTpw4IA8Ho9eeumlfE/FkdfrTSkOoLh885vfdIzv3LlzRtdduHChampq1N7erpqaGl122WWSFLMR/cDAwIzuB2QaPbsAAACQlgMHDmj9+vVRMWOMfvazn82ZJuh+v1/Dw8Py+/3TdqYsBKFQSJWVlY5xdpgEil/4NWgqn8+noaGhor8f8sMS7YYAACAASURBVIueXQAAAJhzrr/+esfeT3Mh0RXuexX+0Dc0NFSQvbCcEl2SFAwGczwTANnQ19enhoYGBQIBSVIgEFBjY2PMCqyZckp0xYsD+UKyCwAAAGnzer2qqKjQY489poqKCpbHFRiWMQKzWzAYVHl5uYaHh+Xz+TQ8PKzy8vKElZvd3d1asGCB3nzzzRzNFMgtkl0AAABI2/DwsE6fPq1HHnlEp0+f5rf7BWZ0dDSlOIDic+rUKW3btk2vvvqqtm3bppMnTyY855577tGZM2fU0NCQ0r3YhRfFgp5dAAAAQIriLVcspPfXxTJPYK4LhULasmWL2tvbY1ZlJTMmkZm+JtCza26hZxcAAAAwh5w4ccIxHgqFcjyT+FjGCBSHlpYWdXZ2qqWlZUZjEunq6tLixYujYkuWLEm6MoueXSgWVHYBAAAAaXCqkDh48GBBNej3+XwaGRmZFi8tLeXDKZCmTFRYhSVTKZXJaqqZVnZ1d3frC1/4gj766KNIbN68eXrllVe0atWqlOaCwkdlFwAAADDHbNy4Uc3Nzeru7pbH45Ek3XLLLXmeVbRYO7INDAzkdiLALLJr1y4dOnRIO3funPG1YiWYJseTGZMrtbW102LGGBJdKDgkuwAAAJC2ubyjV0dHh/bs2aPa2lqNjY1Jks6dOydjTNzqiVwKBoNasmSJpI+XLlZXV8+4GgWYi/x+v4wxam1t1cTEhFpbW2WMkd/vT/uaTpWXU+P9/f1atmxZ1PHly5fnLWl9ySWXqKamRu3t7aqpqdG8efPyMg8gHpJdAAAASFu6O3ohd+rq6tTc3KzXXntNzc3NUZUZczlZCaSqr69PDQ0NcrvdkiS3263GxsaYFZSZUllZqaNHj0bFfvvb3yoYDKb83A3PPdm4kxMnTujIkSPavHmzjhw5ErOHIZBP9OwCAABAytjl74LS0lKNjo5Oi3u93pgVG4Vk8eLFevfdd7V48WKWNgIJZLJ3VqxrSRdeX7u7uyNLA71er86fP+84tqamRkeOHEn6vqFQSJWVlY5xKj4xFT27HBhjqowx3zPG/MoYM2iMscaYJQ7jKowxf2+M+U9jzDljzM+MMX+crXkBAABg5ma6o1e6QqGQ1q1bp5MnT2b1PslySnTFixeK8FLLd999V5J07Nixglp+CRSiTPbOirdBhLVWq1evjvz92LFjKikpcRzb29ub9nPX5fo4HUCiC7NNNpcxLpO0WdIHkn7pNMBceEb+WNINkv5C0u2SSiS9bIypyuLcAAAAMAN1dXU6duxYVGxgYCDqA1o2ZLIxdCZ0dXU5xrOd9MumQkomAoUkmf5ayerq6kq611UwGNT4+HjM4/F+0XDgwAF5PB699NJLkqQHH3xQkvSHf/iH+vWvf63m5mZt3LgxxdkDhS9ryxiNMS5r7cTFr/+HpO9LqrbWDkwac6ukf5Z0rbX25Yux+ZL6Jf2jtXZHovuwjBEAACD3cr2MMZPLhzKpWJdzHjhwQBs2bIiaozFGN910k37605/qgQce0J49e/I4Q6DwZPL5nkwl1uRrbtq0Se+8847eeustlZSURFWPxlvKuGDBAp05cybmPfL9GorCxjJGB+FEVwK3SDoRTnRdPO+MpJ9IujVbcwMAAMDMxFpSEys+U+HG0JPlojF0It/61rcc448//niOZ5Ka9evXT/twbq3VT37yk4ztMgfMNqm+7s102XX4fJ/PpxdeeEG9vb2y1kYSXR6PRzU1NTp9+vS0c8NLG+MluqT4yymBYpbv3RhrJDmloHslXWGMYQ9TAACAAjS510sy8ZmqrKzUs88+GxXbu3evgsFgVu6XrEcffdQx/rWvfS2l6xRaLzKXy1UQyUSgkMRqEh8rfv/996uzs1NNTU3Tjnm93rj3amtriyzbvuuuu9TQ0KBAICBJCgQCamxs1HvvvRdzN8Rs/eIBKBb5TnZdqgs9vaYKp6YrnE4yxjQZYw4bYw7//ve/z9rkAAAAgHj279/vGD948GBK12lpaVFnZ6daWloyMa0Zm5iY0N69e1VdXZ3vqQAFI1aCamrc7/fLGKOf/vSnkqSf/OQn0yolEy1jbGpqUmtrqyYmJvTMM8/o2Wef1eDgoIwxGhoaUnl5edym8vGa2gNzQb6TXUaS0+LmuM98a22btbbeWlt/2WWXZWdmAAAAQALXX3+9Kiqifz9bUVGha6+9Nqnzwx+Kn3nmGUnS008/XVDLB1niBHxsYGBAHo8nKubxeKZt1hHreZPK82nqsu2qqiqtWLFCkrRy5cqEVaCJmtqHud3upOcEFJN8J7tO60J111ThdwxOVV8AAABAwRgcHFRFRYUee+wxVVRUaHBwMOlzYzW1znZz+0RLqABMt3TpUo2NjUXFxsbGMloBWVpaKknTlm0fP35cb7/9tqy16u3t1QsvvJAwKX7rrbeqpqYm0r/LydTkHTBb5DvZ1asLfbumWinpXWvtRzmeDwAAAJLQ398f6R8TVlZWpoGBgfxMKI+Gh4d1+vRpPfLIIzp9+nRK1RsjIyMpxTMlmZ3gAETLZHK6v79fZWVlUbFkX0P9fn/Cnnp+vz+qqX2sOc7F12zMDflOdv1Y0ieNMevCAWNMuaSbLx4DAABAAVq6dOm0CqZz587R4wkAklBZWalz585Fxc6dO5fUphvJ9OwK72A7ual9OLkWruysrq6Oew2gmGU12WWMucMYc4ekz14M3XgxFk5u/VjSryT9ozFmizFmw8WYkfQ32ZwbAAAA0peJnjSpSLYx9GyTrV0anSrzAMSXbCVmrF1av/71r0e+nulrWmtra9xljMFgUOXl5RoeHpbP59Pw8LAWLlyo5uZmvfbaa2publZtbW1S9wKKkclmPwBjTKyL/8Jae83FMZdK+o6k2yT5dCH59VVrbU8y96ivr7eHDx/OwGwBAACQrHjL4LLx/tLn8zl+0CwtLS3qJuqJ/h2bm5v11FNP6YEHHtCePXsyeu/q6uqES5h6enq0atWqjN4XKFalpaUaHR2dFvd6vdNen5ye25NfG0OhkD71qU/p7NmzkVh5ebnefvvthNVdfr9fmzZt0ne+8524lVmbNm1SeXm5fvOb3+iP//iPdfbsWXV0dMS9dqq6u7t1zTXXqLOzk9eKWcgY84a1tj7f80hHViu7rLUmxp9rJo05ba3979baS621AWvtdckmugAAAJAfua606u/vn7ZrmNvtntX9Zowxam1t1cTEhFpbWzO+S2NdXZ22bt2qT3ziEzHHrFmzJmP3A4qd026MkvT666/HPOdP//RPHePBYFDnz5+Pio2NjWnRokXq6uqKO49kljFKUkdHhwKBgLq7uxUIBDKe6JKk22+/XWfOnNGmTZsyfm1gJrJa2ZULVHYBAADkXj4qrcKVSF6vV6Ojo6qurlZfX19W7pUr8Sq7GhoaonZka2xsTFjJkcn7h4+fOHGCvj7ARYFAQENDQ/J4PJGdGbdu3aof/OAHKV0nUVVnMptI+Hw+DQ0NOR7z+/2Or8XxzklFrqt7kR9UdgEAAABZVldXN+v6zZSUlMQ8NjnRJUl79+5Nqnl1Krq6uuJWi1lrtXPnzozeEyhmN9xwg9xudyTRJUlPP/10xisvg8Ggampq1N7eriuvvFI+ny9y/WR2Yww3qHe5Lnzkd7vdCc/B3JOtvpCFgMouAAAApCwUCmnZsmVROzKWlZXp6NGjVAGlIBQKqbKyMqVzUnn/HgqFtGXLFrW3tzv+f4lV/eEkUxUhQLHLRGVrd3e3PvOZz0Q9n40x6u7udux9tX37dj355JNRf4/Xx4/KLsQT/tnwiU98Qvv27YtZnVjMlV0kuwAAAJCW2bisMB/cbrcmJibkdrs1Pj4ul8sll8sVVTkS5vF4pvX5iSdRg/tklkpJ2VlCCRSrTC3jnvr8M8ZoYmJi2rh0ElfZTkaR7CpuLpfL8f/T1MdUMSe7WMYIAACAtMzGZYX5cOutt6q5uVlvvPGGmpubdeuttzomuiTFjE/l9/sz2uA+G0sogWLllOiKF49l8lLFmpqamMnkqUsSXS5XwiWJ2d5EZP/+/Y7xgwcPZuT6yI7wz4ZYCcnZlKiksgsAAAAoMDOtmgiFQnr44YcTNrhPtrIrlXsDs92BAwe0fv36afGDBw/q2muvzfj9CrGyy+kesSrTUDiSec2fsrS2aCu7SHYBAAAABSKZHlrJvH9P9oOu1+tNaVlksX92ADIh10v40lk2mYs5+nw+lZaWamRkJPLfbO3Gi8yYS8kuljECAAAABSLRB8Wenp6M3i+8LApA8rq6uhzjmX5+zkQulhkODw/r8ssv18jIiC6//HISXUUgUbJr8jLXUCgkSSuyO6Ps4acbAAAAUCScdmlz4na7k4rH6/mT7DWBuaa2tlYrV66MitXU1CT9/ExVf3+/AoFAVKysrEwDAwMxz7nlllsc43/+53+e1hxCoZDWrVunkydPSrqQNDHGqLe3V5LU29sbiaFwvf/++3GPj46ORr7etWuXJM3L7oyyh2QXAAAAUARSaRA/Pj6eVDwYDGrJkiWSPv6N/uWXX+547vHjx5O+PzCb+f1+vfXWW1Gx3t7etDeAmKy7u1sLFizQm2++GYktXbpUg4ODUePOnTun6urqmNfp6+ubVrnpcrlSSnBPtmvXLh06dEg7d+6UdKG6bfHixVFjlixZUlDVbZgumZ8jkzc4KWYkuwAAAIAicOLEiai/HzhwQB6PRy+99NK0sXv27HG8Rltb27TY1F016+s/bs8yuUoj1k5xwFwT3h0xXO3odrsT7o6YrHvuuUdnzpxRQ0ND1P2qqqqixlVVVcW9XzAYVFNTk1wul3w+n1wulx544IGUn8exdnatq6vTsWPHosYODAxozZo1KV0fubdx40Z98YtfdDy2b9++yON7ajVhsaFBPQAAABDDgQMHdOONN2r//v1Z2WVtqlSaSpeWlmp0dFRer3da8+qZNqfetGlT5MNyW1ubQqGQOjo6Ep4HzAXp7I6YSLznrM/nS+t+k5/Hd911l95++201NTXpqaeeSnpe4Z1df/jDH2piYkIul0t333239u3bF7XkLSxe03wUFqfHXPjnw/bt29XW1qaJiQlrrS3KIimSXQAAAEAM4V3QcvUB7oc//KHuvvvuafF9+/bpjjvukJRcIivXu8UBc0k6uyMm0t3drdtuuy2qWmrJkiV68cUXtXr16pjnJft8nvyakOicUCikLVu2qL29XdXV1Y7fU2lpqS6//HIdPXo0Elu+fLk6OzupAi0SbrdbHo9H27dvV2trq8bGxiJL3cOJ0j179rxlra3J81TTUpQZOgAAACCbwo2Wwx9oR0ZGctJ8ecuWLY7xcKIrWV1dXfJ4PFExj8cT6acztdk0gNy6++67ZYzRl7/8ZUnS2rVrYy4LnLxD3mRer1c7duyQMUZf/epXHcc4vW4lei1raWlRZ2enWlpaYibGRkZGohJdkvTb3/42bh8xFJbx8XGNjIzoiSee0MjISFRPx46ODu3evVuS0itVLABUdgEAAABT5LMyKt5v21OZm8fjiTrP7XZrbGxMknT11Vfrl7/8pdatW6ef//znmf8mgFksVmVXSUmJAoGAOjs7E+7MOLXSKhQKaenSpRoZGZG1VsYYlZaWqr+/X9ZaXXHFFZHnr3Th+X3TTTfpn//5n6OuE+8+U00dH2t55lTLly/Xvn37dNNNNykUCml8fFxut1vBYFCvv/46lV2ziDHmDWttfeKRhYfKLgAAAKCAxPtteyoWLlyompoatbe3q6amRgsXLow0m/7lL38pSfrFL34hY0xGdpED5gqnRJcknT9/XmfOnInb3y9WpVVlZaXuvfdeGWPk8/lkjNF9992nRYsWKRgMqqSkRJIiFZtjY2NRia5Y105Fson83/72t6qtrdVNN90UeX0aHx/XzTffTKILBYPKLgAAAGCKQu551d3drbq6umnxtrY23X///XHPLeTvCygWqSSUpj6v4p27ceNGBYNBLViwQN/85jdVW1urrq4uSdHN5mtra5O+ZyrP+VAopKuvvnra8sRUzKRJPwoPlV0AAADALLJ//37H+MGDB3M8k+lifdD9y7/8y4Tn7tmzxzHe1tY2ozkBc0ms14dkxEo+GWMifZK++c1vSrqQ2A4LH4vXrF5SpAIs0VydXsuCwWDUUsl4PB6PGhoaIt+PMUaNjY3q7+9P6nwg26jsAgAAABzE25Y9HybvkPaZz3xGoVAo5tgTJ05Exk5eVkRlF5AZyVZ3OVVPVVZWThsXCoUUDAaTuo7L5Yr7fE2mmizV/l6p4LVk9qCyCwAAAJhlSktLVVFRoccee0wVFRUqLS3N63zuuusudXZ26q677tKJEyfijn3wwQfV2dmpHTt25Gh2wNwy+fUhFbESWqn0unr//fdTumeyr2VdXV1avHhx0ted2uvP7/dHdnwF8o3KLgAAAKBAdHd365prronazS3ZHdLiCffRobILyLxYffR6enocd2Wc3H+rra1NoVBIHR0dKT0/Y409ePBg3Ab58aRa2bVy5Uq99dZbkb/X1NToyJEjad0b6XP6uZEpVHYBAAAAmLF77rlHZ86cUUNDQyQ200TX5Gt4vV7H47HiABJzSnRJitlfa3L/rd27d6ujo0PSxzstTuUU37hxo+PYdBJdO3bsiPTccrmiUwTGmGmxsHCiK7zj6+nTp1O+N2Zu/fr1OnPmjK6//vqMXO+v/uqvZIzRN77xjYxcL1+o7AIAAADyLFO9cmLxer0aGRmJ2y8olWVUAD72wx/+UHffffe0+OOPP66HH3446euk+vz0+XwKBAJqaWnRt7/9bQ0ODqaVHJ/p60+x5xSKVbYqdadct2gru0h2AQAAAHnW3d2t2267TceOHYvElixZohdffDHh7mvJCr/vDwQCGhoaksfj0djYmAKBgM6dO5eRewBzVWlpqUZHR6NiVVVVeu+991K6TjjR4HK5NDExISl7ySQa0he3TCe7Yl3PWpvd38ZkCcsYAQAAgDyrra1VWVlZVKysrEyrVq3Snj17HM9pa2tTMBhUTU2NSkpKkr7XDTfcoObmZh0+fFjNzc3asGHDjOYOQBobG5sWO378uIwxKSWVNm7cqObmZv36179Wc3NzzOWKueR2u/M9BSBlVHYBAAAABaCyslKXXnqpvvGNb2jnzp06ffp0ZNdFpw/L+/bt0x133CHpwjKnhx9+WM8++2zM6xf7+36g0BXbBhCtra1qbm5OOC4UCsXcRVIqzO9tLqCyKz4quwAAAIACcOLECR05ckSbN2/WkSNHIomuMJfLFdW4/s4774x8XVlZGTfRlQmhUEjr1q3TyZMns3ofoFjFqrBMpfIyl7Zv357UuEWLFqmrq8vxWE9PTyanhDimvgbHq7gzxuihhx7SggUL9Oabb8a9bqrVh8WCZBcAAABQ4Ky1mpiYmJbQCn9IibXUMezgwYMznsOuXbt06NAh7dy5c8bXAmajWLsWxooXg/AyytraWsfjq1atyuV0Zr3u7u6YCaoHH3xQnZ2d2rFjhyQl7Af3xBNP6MyZMwmXwhZqMnamWMYIAAAAFIFES1aytYTK7/c77vDm8/k0NDSU9nWB2abYljGGud3uSDP8qfHJvcjiLbUOC+9MOXmZNZJXVVWl999/P2pzg3ivweG4MSapx5jTGJ/Pp5GRkTin2KLM1pLsAgAAAIpAuskur9cb74NMQuF+YC+88IKGhobk9/u1adMmfec739GiRYvSvi4w2xRrsitWsqO0tNQxyRJPOHHmcrk0Pj6eqSnOeukuI9y4caOCwaCamppiVt9N5vQ4THDvN6y19WlNLs88+Z4AAAAAgMROnDihysrKafFQKCQpdsJrJokuSQoGgyovL49UcQ0NDam8vJxEF4CIqa89ExMTkVghJ/qKXUdHR+Tra665Rj//+c8dx3k8Hr3xxhuOx0pKSnT+/Pmk48WiKMvRAAAAgLkm1m5ok5NOwWBQNTU1am9vV01NTdwd1JLl9/v15JNPRsVaW1vl9/tnfG0AyLS5upnGyy+/HPPY2NiYVq9erQMHDsjj8eill16KHJuNveYkkl0AAABA0di4caOam5vV3d2t5ubmaY2HE+3omI6+vj41NDTI47mwKMTj8aixsVH9/f0zvjYwm7BjYWG4//771dnZqaampnxPJafuvvvuhGM2bNig8fFx3XDDDZFYrNfygYGBTE0tL+jZNcuFQiFt2bJF7e3tlJoDAAAgZYXcoL61tVXNzc1qa2vT/fffn9e5AJJz/6NC/8wdCoW0bNkyDQ4ORmJlZWU6evRo0p8hYzW5z2XvrkJ+rUrkC1/4gg4dOjQtvm/fPs2fP1/r16+fduz222/X888/H/l7un2/Jje6d4gXbc8uKrtmObaIBgAAwEz09fWpqqpKbrdb0oUPtVVVVQVR2dXc3CxJc66CA4UrG0uJsy0YDGrhwoWSLmxoIUkLFy5MqVji+PHjjvH3339/5hNMUqxm+sk02Xda3ucknSWSO3bskDFGX/3qV6cd8/v9MsY4Jrok6Y477nBMdEnSP/3TP0m6kORKN9ElxU7GFnqSNhGSXbNUaWmpjDFqbW3VxMSEWltbZYyRz+fL99QAAABQRILBoG666aZIdcb4+LhuvvnmvK4acPpwN9MPfEAmZGMpcS7U1dWpublZr732mpqbm5Pa2W+yZHoKForbbrtNxhjdcccdkdiNN96o8fFxXXfddTLG6Pbbb5cxRvfff7+MMfr+978vSVq/fr06Ozu1YcOGyLnd3d1asGCB3nzzTcf7fe9735Mk/e3f/u20Y4kSSole03jNi23OJ7viZVmzIRvN8sLXbG9vj2SjlyxZ4jj2lltuyWqzvvBcLr/8chljVFVVFfOJv3LlShljFAwGo15AAAAAUFhOnTql+fPnS5Lmz58/4/eSTh82Y0n0QRLAzHV0dGj37t1avXq1du/eHbXLX7KcegqGiy4K6bPeiy++KOlCZVQ4ST51qWX4+//7v/97SReqR40xOnLkiCTpzTfflDFGfr9fa9as0ZkzZ/S5z30u6hrJJOVnultuIgsWLIh7/Otf/3rMOWR7bllnrS3qP5/97GdtIt/61resJPv4449POyYp6o/b7baSbG1t7bSxS5YssZLs0qVLE94zlu3bt1uXy2W3b9+e9jXCTpw4Ya+++mrr9/unfR+J/rhcrmnX6+rqsvPnz7dVVVVWkr3yyiunjUn0b1BeXu54P4/HM22s07jJ31coFHL8u5PwvD75yU/a+fPn256enmn/TvHOBwAAgLN47ykzcc1EfD6flWT9fn9ac+K9IJA/mXi9SNatt97q+Jpw++23T5tLtv8kup/Tv0825pDMPLxer+Oxi/HDtgDyPun8yfsEZvpn/vz5kf8ZTU1N1km6D6qwP/qjP0rqB/vUpNrkH6zhH9JT/5SWlk5L7KxZs8Z+/vOfj/kDOXxdl8uVkSdB+MnvdL2wVatWxf03iPX9zeRPWVmZtdbaioqKSOyqq66Kum84Qef0xFy7dq2VFPm+KioqHP89AQAAEFu8D0KpSua9dzJj9+/f7xg/ePDgtOt85StfsZLs1q1b0/n2AaQhled6tu+bzJxy/Sfss5/97IyvtW/fvsjXK1asSOn7TuLfpWiTXUW/G6MxJuobmPz95GL9ajL3mz9/vn7+859r/fr1+v3vf+84ZsGCBfrggw905513RnZUCO+KcOONN+pf//Vf9dBDD+m73/2uPB5Pzna0SIa1NqdrhX0+nxYtWpT2VqjFsBsHAABAoQiFQvrUpz6ls2fPRmLl5eV6++23U+7HE+8949TPJYnGJtr1rph3ZgOKWXiX1FhSyUHs2LFD3/ve9yKfhevr6/XGG2/oqquu0r//+79PGx9+XXC5XJHdIauqqnT8+HEFAoGoHSfzpaSkRKOjo5Iyk7OYaU6nu7tbdXV10+I9PT1avXp10e7GOOuSXWG5SMAsXbpUv/vd7+Z0U7jwv8HKlSv1H//xH/meTtKK/XEPAACQS4FAQENDQ/J4PBobG1MgENC5c+dSvk6s981ut1tjY2Nxx0xmrZXP51MgEFBLS4u+/e1va3BwMCq55fP5HHvOlJaWJrU7G4D0uN3uSKLJydTPYqtXr9abb74ZeX258sor5Xa7E36+tNbq8ssvj7kTZPg+hfZ5PRQKxd2lc3Kizinu9Xq1fft2tba2amxsLCOFMOXl5frwww8jf7/kkkt09uxZGWOKNtk1axvU5+IB/bvf/S7r9yh04X+Dt956K88zSV5NTU2+pwAAAFBUbrjhBjU3N+vw4cP64he/qMHBwchqhGQk2ikx1odVJwcPHpQkDQ8P6/Tp03rkkUd0+vTpaQmsZJsuZ2MDKWAuCj/P4yW63G535OvwpmbhzSfCCe933nknqUIKY0zc145C3KHV5XIlrIidnLya/O81Pj4ua61GRkb0xBNPaGRkJGMrvubNm6eamhq1t7erpqZG8+bNy8h18ynvyS5jzOXGmOeNMWeMMWeNMR3GmCtmcs1LLrkkU9NLytq1a3N6v0Ln9XpTPmfhwoVZmImz8A4aAAAASM7kndp+9rOfSZLuvPPOjF1/8oe/EydOxB177bXXZuSef/InfxLZPbyzs1Of/vSnM3JdYK7as2dPwjGTkzOpJLlzoba2Nqrq7MYbb4x8PbkX1KpVqxzPT2ZZ96233ipJ+ou/+AvH4w899FDU/cbGxj5uuJ5FJ06c0JEjR7R582YdOXIk4etwMchrsssYE5D0kqRPSdoq6cuSlkt62RhTlu51J/cTSJa1VsFgMCqbGau0cOqD7d/+7d/SnWrG1dbWzuj8yZnjeE8qp2Ph2MjIyLTmcPFYa3Xq1KmMPYF9Pl9GrgMAAICPOVVJhGPGmEiVxpe//OVp58b75bDf75d0oTePMUbXXXed47iKiorIPb1er4wxWrFihePYu+++O+H38qtfQ9CPfgAADN9JREFU/UqSIlUo//Vf/yVjTGQ+AFKzffv2hGP27dtXkBVXktTV1SXp48+1//Iv/+L4ebanp8fx/FAoFPf6jz/+uDo6OiRJf/d3f+c45rvf/W6q00YMee3ZZYx5UNJ3Ja2w1h69GKuW9FtJj1hrE/6fdurZNblfV3jdbzh+8RxJ0uLFi3Xs2LGoY1NVVlbq0ksv1Te+8Q3t3LlTp0+fdsxyFsKTNVZTT6fvM3ystrZW3d3djueHud1ueTyeGa0LTuXffPL4yR577DF9+9vf1gcffBDzvKnfX7i5f6L7AQAAIL5U3u86ve9Kp+F8d3e32traFAqF1NHREbMpfXV1tQYGBjLWU5f3jUB64j33wksc4/XZyrXwpnBSas/78Pc59fx0NuFI5jN5vhRzz658J7sOSvJZa/90SvwXkmStXZfoGvX19fbw4cNZmmFmhB/Ebrc7kiSK9wbA6/VGdmeIxSmJN5uQrAIAACgs6SSQcr1TeqbwvhNIn9Nz3ev1RhVNzPT1oKamRr29vTGPx0os5epz9Gz5PFvMya589+yqkeTUQKlX0soczyVrkl1vG28Z4NQ/58+fz8na3XxJVDoKAACA3HruuedmdH74Pd2yZcsyNKPsiLVECUBqAoFA5OtkmqlXVFSotLR0Wnuh8LHHHnssMubIkSNxPy9PlevP0Xyezb98V3aNSvqutbZlSvz/ltRirfXEOK9JUpMkXXHFFZ8NL4sDAAAAkD2pVGM0NTXpqaeeysh1/uzP/ixhH65M4UMpAFxAZdfMOP00ifvTz1rbZq2tt9bWX3bZZVmaFgAAAIDJXC6XvF6vHnzwwYQ7cMdKdE3m8Tj+bnvadbZs2ZL0HJ04VYpIUmlpaeTreBtUAQCKS+KfLtn1gaRLHeIVF48BAAAAKBCTlyE98cQTaffdidVPJx6Xy+W4aZJTf1xJWrFihd5++21JitpgavPmzWnNGQBQPPKd7OrVhb5dU62U9FaO5wIAAAAgBZla8pfMdaYm2jI9BwDA7JHvZYw/lvR5Y8zScMAYs0TSn148BgAAAAAAACQt38mu70sakPSiMeZWY8wtkl6U9J6kxIv8AQAAAAAAgEnymuyy1p6TdK2kdyT9L0l7JfVLutZa+1E+5wYAAAAAAIDik++eXbLWvivp9nzPAwAAAAAAAMUv38sYAQAAAAAAgIwh2QUAAAAAAIBZwxT7Vr3GmA8lvZ3veSDv/kDSf+Z7EphzeNxB4nGAC3gcQOJxgPzhsQeJxwEyb4W19pJ8TyIdee/ZlQFvW2vr8z0J5Jcx5jCPA+QajztIPA5wAY8DSDwOkD889iDxOEDmGWMO53sO6WIZIwAAAAAAAGYNkl0AAAAAAACYNWZDsqst3xNAQeBxgHzgcQeJxwEu4HEAiccB8ofHHiQeB8i8on1MFX2DegAAAAAAACBsNlR2AQAAAAAAAJJIdgHA/9/e+cd6VdZx/PWm0jCoxcgu/VDYciuBLbOaNsCvLP+IfqzRpsytxpbLSbNMSddsis0lrbLhVnPzD1n0g7Zmi1wmiiJsICnlEjdM1oUlIBRlQAEGe/rjeb7d0+Hc7znf7929t3PO+7U9u99zns/5PM9zvu/zvvc+55cxxhhjjDHGmAbhyS5jjDHGGGOMMcYY0xhKJ7skfU3SM5KOSvqLpF9JmpeLkaRVkg5IOiFps6S5uZgvSHpS0quSgqTZBW3tTXXZsrpCH+dLeiq1vV/SHZKUqZ8l6SeSdks6I2ltWc7MtiskDUs6KWmnpIW5+qWSHk37JkjqVM1dJ9qsA0mLJG1IOYOk5QUxawv6/HSV/KY3DdHeUkkbU/+PSdoh6VMVx28Pot06sAeN0BAdXCFpm6QjKWa3pJUVx28/oN06sB9MHk3QXS52gaTTknZVHL/9h3brwP4zPjRBU5I6BXmDpPdWyD1u3lLlyq4O8APgI8Bi4DTwuKQZmZhbgVuAG4EPAYeBxyRNz8ScB2wEVpW09w1gVqbc3StY0puBx4BDqe0vAV8Fbs6EnQv8FVgN7ChpP5v7GmAN8E3gEmAb8IikCzJhb0rrbz47Q6Po0FIdANOAXcCXgRM94h7P9XlJH22Y0elQf+1dATwBfJzoJb8GfpE384Lc9qAROrRUB9iDsnSovw6OA/cBi4CLU867JK0oyW0/GKFDS3WA/WAy6VB/3XVj3wr8ENhU0oduvP1nhA4t1QH2n/GiQ0M0BczN5X6pJPf4eksIoa9CFPkZ4JNpWcBB4PZMzFTgGHB9wfYfBAIwu6BuL7Cyz/7cABwFpmbWfR3YT3rbZC7+YWBtxdw7gAdy614C7imInZnG1el3n9axtEkHue2OA8sL1q8FHp7s76UNpe7ay8T8FvhuSW57kHWQj7cHNVMHDwE/LcltP7AO8vH2g0ksddZd0tqdxH+Kd1XIbf+xDvLt2H+sqf9qijhhF4CZfeYeV28Z5Jld04lXhP09Lc8BhoiziACEEE4AW4izk/2yUvGS7uck3S7pnJL4y4Gtqc0ujwLvAGYP0D4Aqd1LyYwrsZHBxtU0WqGDPlkg6bCkP0p6QNL5E9Ru22iK9qYzMoazsAeV0god9EkbPaj2OpB0SerbU6MltR+U0god9Ekb/WCiqaXu0tWDQ5RczZGJt//0phU66BP7z9iopaYSz0o6KGmTpCt7JZ0Ib3n9ANusAZ4DtqflofTzUC7uEPDOPnPfB/weOAJ8mHi72Rzguh7bDAEvF7TdrRvusw9dZgKvo3hcHx0wZ5Noiw6q8hvi2ZFh4kF/N/CEpEtDCKfGue22UXvtSfoi8C5gXY+89qDetEUHVWmrB9VWB5JeBt5G/FvsrhDC/T3y2g960xYdVKWtfjDR1E53kuYTr+S5LIRwRsWPccpj/+lNW3RQFfvP2KmdpohXnt0APAOcA3wW2CSpE0LYMkrecfeWvia7JN0LLAAWhBDO5KpDPrxgXU9CCPdmFv8g6SjwM0m3hRCOSHoBuDDVbw0hfKxH20XrC0nPS3kks+p64MkeufsaV9Nokw5CCD+u2Of1mcXnJe0E9hGfzfNQlRymnCZoT9JngG8Dy0II+9I6e1AftEkH9qDRaYAOFhJvVbgM+Jak4RDCOvtBf7RJB/aD/x/qqDtJ5wLribcwFZ4Etv/0R5t0YP+ZGOqoqZT3ReDFTP12xQfkrwS2TJa3VJ7skvQ9YBlwZQjhT5mqV9LPIeDPmfXnc/YsXb90HyL+HuLs4xLgDWld9zK6VxiZ7cy2TR/tPwu8P7N8CDhFvFe2KPdYx1VbWqiDgQghHEhnay8aNIf5X5qgvTTBsQ74XAhhQ6bKHlSRFupgIJruQU3QQeafjOclvZ34zJR12A8q00IdDETT/WCiqbHuZhFfhvCgpAfT+inEF72dTjm3Yv+pRAt1MBD2n+rUWFO9ci9Lnyflb5tKz+yStAa4FlgcQtidqx4m7oCrMvFvJJ6p2jbG/nV3yEGAEMK+EMKeVPanuu3AwtRml6uAA8QHsJUSQjiRybsnhHAshPAasJPMuDK5xzquWtJGHQzaYUkziZeVHhw0hxmhCdqTdDXwI+LDPH+ebcQeVI026mDQDjfZg5qggwKmEN8YbD+oSBt1MGiHm+wHE03NdbcfmJ9ydcv9wJ70eZv9pxpt1MGgHbb/VKPmmuqVu5t3crwllD8h//vEp+8vJs66dcu0TMxtKWYpMI94aeQBYHomZigN+FriZWlL0vKMVH858JW0bg5wNfFg/GVJ/95C/PLXp7aXpr7ckovrHsxbgA3p88Ulua8BXiPew/o+4v2zx4ELMzEzUq5OGtd1aXmobN/WqbRcB9My2/0LuCN9viBT/53U99lJC9uJ9zZPr7J/XZqtPeJZjX8TX9WcHcOMktz2IOug6zH2oObo4EbgE8Sz3BcBn08xq0ty2w+sg+7xbj+w7gbSXcE2q6j2Nkb7j3Vg/7Gmev1Ouwn4NPH32VzgntSHpSW5x9Vbquz8MEpZlYlROkgOAieJb5KZV3AQFeVZnuo/ADwNvEq8ZG532ua8Cn2cT5y8OJn6cCdnv1q1qO29FXKvIM5YniLOPC7K1S8v2z9NKG3WQebAype1qX4q8Y0Uh4kH6z7ia3ffPdnfWxNKE7QHbB6l7c0VctuDWq4D7EFN08FNwAvAP4F/AL8jHudTKuS2H7RcB9gPrLsx6K4gfhUVJjlSrP2n5Tqw/1hTo2kKuJV4deAJ4G/EW2GXVBz/uHmLUgJjjDHGGGOMMcYYY2pPpWd2GWOMMcYYY4wxxhhTBzzZZYwxxhhjjDHGGGMagye7jDHGGGOMMcYYY0xj8GSXMcYYY4wxxhhjjGkMnuwyxhhjjDHGGGOMMY3Bk13GGGOMMcYYY4wxpjF4sssYY4wxxhhjjDHGNAZPdhljjDHGGGOMMcaYxvAfeqthJ8sZ1jsAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1440x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig,ax=plt.subplots(1,1,figsize=(20,4))\n",
"ax.plot(df.loc[iidfnd,['dts']],df.loc[iidfnd,['Chlorophyll_ug']],'k*')\n",
"ax.set_xlim(dt.datetime(2015,1,1),dt.datetime(2015,5,1))\n",
"dt.datetime(2014,12,31)+dt.timedelta(days=65),dt.datetime(2014,12,31)+dt.timedelta(days=75)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"with nc.Dataset('/ocean/eolson/MEOPAR/NEMO-forcing/grid/mesh_mask201702_noLPE.nc') as fm:\n",
" tmask=np.copy(fm.variables['tmask'])\n",
" umask=np.copy(fm.variables['umask'])\n",
" vmask=np.copy(fm.variables['vmask'])\n",
" navlon=np.copy(fm.variables['nav_lon'])\n",
" navlat=np.copy(fm.variables['nav_lat'])\n",
" e3t_0=np.copy(fm.variables['e3t_0'])\n",
" e3u_0=np.copy(fm.variables['e3u_0'])\n",
" e3v_0=np.copy(fm.variables['e3v_0'])\n",
" e1t=np.copy(fm.variables['e1t'])\n",
" e2t=np.copy(fm.variables['e2t'])\n",
" e1v=np.copy(fm.variables['e1v'])\n",
" e2u=np.copy(fm.variables['e2u'])\n",
" A=fm.variables['e1t'][0,:,:]*fm.variables['e2t'][0,:,:]*tmask[0,0,:,:]\n",
"\n",
"t0=dt.datetime(2015,1,1)\n",
"\n",
"#stm=np.shape(tmask)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/data/eolson/results/MEOPAR/SS36runs/calcFiles/comparePhytoN/ts_hindcast2015_SentryShoal.nc\n",
"/data/eolson/results/MEOPAR/SS36runs/calcFiles/comparePhytoN/ts_hindcast2015_S3.nc\n",
"/data/eolson/results/MEOPAR/SS36runs/calcFiles/comparePhytoN/ts_hindcast2015_CentralNode.nc\n",
"/data/eolson/results/MEOPAR/SS36runs/calcFiles/comparePhytoN/ts_hindcast2015_CentralSJDF.nc\n",
"/data/eolson/MEOPAR/SS36runs/CedarRuns/rev_noMrubrum/ts_rev_noMrubrum_SentryShoal.nc\n",
"/data/eolson/MEOPAR/SS36runs/CedarRuns/rev_noMrubrum/ts_rev_noMrubrum_S3.nc\n",
"/data/eolson/MEOPAR/SS36runs/CedarRuns/rev_noMrubrum/ts_rev_noMrubrum_CentralNode.nc\n",
"/data/eolson/MEOPAR/SS36runs/CedarRuns/rev_noMrubrum/ts_rev_noMrubrum_CentralSJDF.nc\n"
]
}
],
"source": [
"saveloc='/data/eolson/MEOPAR/SS36runs/calcFiles/comparePhytoN/'\n",
"baseloc='/data/eolson/MEOPAR/SS36runs/CedarRuns/'\n",
"dirname1='hindcast2015'\n",
"dirname2='rev_noMrubrum'\n",
"dirnames=(dirname1,dirname2)\n",
"#dirnames=('spring2015_NewSink','spring2015_slowPP','spring2015_KhT','spring2015_diatHS')\n",
"varNameDict={'Sentry Shoal':'SentryShoal', 'S3':'S3', 'Central node':'CentralNode', 'Central SJDF':'CentralSJDF'}\n",
"with open('/data/eolson/results/MEOPAR/analysis-elise-2/notebooks/bioTuning/spathsMaster.txt') as f:\n",
" spaths = dict(x.strip().split() for x in f)\n",
"#spaths={'spring2015_NewSink':'/data/eolson/results/MEOPAR/SS36runs/CedarRuns/spring2015_NewSink/',\n",
"# 'spring2015_KhT':'/data/eolson/results/MEOPAR/SS36runs/CedarRuns/spring2015_KhT/',\n",
"# 'spring2015_uzoo':'/data/eolson/results/MEOPAR/SS36runs/CedarRuns/spring2015_uzoo/',\n",
"# 'spring2015_uzoo2':'/data/eolson/results/MEOPAR/SS36runs/CedarRuns/spring2015_uzoo2/',\n",
"# 'spring2015_slowPP':'/data/eolson/results/MEOPAR/SS36runs/OrcinusRuns/spring2015_slowPP/',\n",
"# 'spring2015_lowMuNano':'/data/eolson/results/MEOPAR/SS36runs/CedarRuns/spring2015_lowMuNano/',\n",
"# 'spring2015_SMELTBFastSink':'/data/eolson/results/MEOPAR/SS36runs/OrcinusRuns/spring2015_SMELTBFastSink/',\n",
"# 'spring2015_uzPref':'/data/eolson/results/MEOPAR/SS36runs/CedarRuns/spring2015_uzPref/',\n",
"# 'spring2015_hiNH':'/data/eolson/results/MEOPAR/SS36runs/CedarRuns/spring2015_hiNH/'}\n",
"ff=dict()\n",
"for idir in dirnames:\n",
" ff[idir]=dict()\n",
" for pl in plist:\n",
" print(spaths[idir]+'ts_'+idir+'_'+varNameDict[pl]+'.nc')\n",
" ff[idir][pl]=nc.Dataset(spaths[idir]+'ts_'+idir+'_'+varNameDict[pl]+'.nc')\n",
" #try:\n",
" # pl='Total'\n",
" # ff[idir][pl]=nc.Dataset(spaths[idir]+'ts_'+idir+'_'+varNameDict[pl]+'.nc')\n",
" #except:\n",
" # pass"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'/data/eolson/MEOPAR/SS36runs/CedarRuns/rev_noMrubrum/ts_rev_noMrubrum_CentralSJDF.nc'"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"spaths[idir]+'ts_'+idir+'_'+varNameDict[pl]+'.nc'"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1900-01-01 00:00:00\n",
"1900-01-01 00:00:00\n"
]
}
],
"source": [
"times=dict()\n",
"for idir in dirnames:\n",
" f=ff[idir]['S3']\n",
" torig=dt.datetime.strptime(f.variables['time_centered'].time_origin,'%Y-%m-%d %H:%M:%S')\n",
" print(torig)\n",
" times[idir]=np.array([torig + dt.timedelta(seconds=ii) for ii in f.variables['time_centered'][:]])"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"ename": "IndentationError",
"evalue": "unexpected indent (<ipython-input-19-74c615b5aa18>, line 2)",
"output_type": "error",
"traceback": [
"\u001b[0;36m File \u001b[0;32m\"<ipython-input-19-74c615b5aa18>\"\u001b[0;36m, line \u001b[0;32m2\u001b[0m\n\u001b[0;31m dirname2:{'diatoms':'g','flagellates':'blueviolet','ciliates':'orange'}}\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mIndentationError\u001b[0m\u001b[0;31m:\u001b[0m unexpected indent\n"
]
}
],
"source": [
"lcol={dirname1:{'diatoms':'darkgreen','flagellates':'mediumblue','ciliates':'maroon'},\n",
" dirname2:{'diatoms':'g','flagellates':'blueviolet','ciliates':'orange'}}\n",
"lsty={dirname1:'-',}\n",
" dirname2:'--'}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tmins=list()\n",
"tmaxs=list()\n",
"for idir in dirnames:\n",
" tmins.append(times[idir][0])\n",
" tmaxs.append(times[idir][-1])\n",
"xl=(np.min(np.array(tmins)),np.max(np.array(tmaxs)))\n",
"xl=(dt.datetime(2015, 5, 1, 0, 30), dt.datetime(2015, 9, 27, 23, 30))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"xl"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"sumVarTr={'diatoms':'diatSum','flagellates':'flagSum','ciliates':'myriSum'}"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yearsFmt = mdates.DateFormatter('%b %d')\n",
"fig,ax=plt.subplots(len(plist),1,figsize=(16,3*(len(plist))))\n",
"pp=dict()\n",
"for ii in range(0,len(plist)):\n",
" pl=plist[ii]\n",
" if pl=='S3':\n",
" ax[ii].plot(df.loc[iidfnd,['dts']],[float(ik[0]) for ik in df.loc[iidfnd,['Chlorophyll_ug']].values],'.',color='lightgray',alpha=.5)\n",
" pp[ii]=dict()\n",
" for idir in dirnames:\n",
" f=ff[idir][pl]\n",
" for var in ('diatoms','flagellates','ciliates'):\n",
" try:\n",
" pp[ii][var]=ax[ii].plot(times[idir],2.0*np.sum(f.variables[var][:,:3,0,0]*f.variables['e3t'][:,:3,0,0],1)/np.sum(f.variables['e3t'][:,:3,0,0],1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.5,label=idir+' '+var)\n",
" except:\n",
" pp[ii][var]=ax[ii].plot(times[idir],2.0*np.sum(f.variables[var][:,:3,0,0]*np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1)/np.sum(np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.5,label=idir+' '+var)\n",
" try:\n",
" pp[ii]['tot']=ax[ii].plot(times[idir],2.0*np.sum((f.variables['diatoms'][:,:3,0,0]+f.variables['flagellates'][:,:3,0,0]+f.variables['ciliates'][:,:3,0,0])*f.variables['e3t'][:,:3,0,0],1)/np.sum(f.variables['e3t'][:,:3,0,0],1),\n",
" linestyle=lsty[idir],color='k',alpha=.6,label=idir+' '+'total')\n",
" except:\n",
" pp[ii]['tot']=ax[ii].plot(times[idir],2.0*np.sum((f.variables['diatoms'][:,:3,0,0]+f.variables['flagellates'][:,:3,0,0]+f.variables['ciliates'][:,:3,0,0])*np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1)/np.sum(np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1),\n",
" linestyle=lsty[idir],color='k',alpha=.6,label=idir+' '+'total')\n",
" ax[ii].set_title(pl)\n",
" ax[ii].set_ylabel('chl ($\\mu$g/L)')\n",
" ax[ii].xaxis.set_major_formatter(yearsFmt)\n",
" ax[ii].set_xlim(xl)\n",
" ax[ii].set_ylim(0,40)\n",
" if ii==0:\n",
" ax[0].legend(bbox_to_anchor=(1,1),fontsize=12)\n",
"\n",
"plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yearsFmt = mdates.DateFormatter('%b %d')\n",
"fig,ax=plt.subplots(len(plist),1,figsize=(16,3*(len(plist))))\n",
"pp=dict()\n",
"for ii in range(0,len(plist)):\n",
" pl=plist[ii]\n",
" if pl=='S3':\n",
" ax[ii].plot(df.loc[iidfnd,['dts']],[float(ik[0]) for ik in df.loc[iidfnd,['Chlorophyll_ug']].values],'.',color='lightgray',alpha=.5)\n",
" pp[ii]=dict()\n",
" for idir in dirnames:\n",
" f=ff[idir][pl]\n",
" for var in ('diatoms',):#,'flagellates','ciliates'):\n",
" try:\n",
" pp[ii][var]=ax[ii].plot(times[idir],2.0*np.sum(f.variables[var][:,:3,0,0]*f.variables['e3t'][:,:3,0,0],1)/np.sum(f.variables['e3t'][:,:3,0,0],1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.5,label=idir+' '+var)\n",
" except:\n",
" pp[ii][var]=ax[ii].plot(times[idir],2.0*np.sum(f.variables[var][:,:3,0,0]*np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1)/np.sum(np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.5,label=idir+' '+var)\n",
" #try:\n",
" # pp[ii]['tot']=ax[ii].plot(times[idir],2.0*np.sum((f.variables['diatoms'][:,:3,0,0]+f.variables['flagellates'][:,:3,0,0]+f.variables['ciliates'][:,:3,0,0])*f.variables['e3t'][:,:3,0,0],1)/np.sum(f.variables['e3t'][:,:3,0,0],1),\n",
" # linestyle=lsty[idir],color='k',alpha=.6,label=idir+' '+'total')\n",
" #except:\n",
" # pp[ii]['tot']=ax[ii].plot(times[idir],2.0*np.sum((f.variables['diatoms'][:,:3,0,0]+f.variables['flagellates'][:,:3,0,0]+f.variables['ciliates'][:,:3,0,0])*np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1)/np.sum(np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1),\n",
" # linestyle=lsty[idir],color='k',alpha=.6,label=idir+' '+'total')\n",
" ax[ii].set_title(pl)\n",
" ax[ii].set_ylabel('chl ($\\mu$g/L)')\n",
" ax[ii].xaxis.set_major_formatter(yearsFmt)\n",
" ax[ii].set_xlim(xl)\n",
" ax[ii].set_ylim(0,40)\n",
" if ii==0:\n",
" ax[0].legend(bbox_to_anchor=(1,1),fontsize=12)\n",
"\n",
"plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yearsFmt = mdates.DateFormatter('%b %d')\n",
"fig,ax=plt.subplots(len(plist),1,figsize=(16,3*(len(plist))))\n",
"pp=dict()\n",
"for ii in range(0,len(plist)):\n",
" pl=plist[ii]\n",
" if pl=='S3':\n",
" ax[ii].plot(df.loc[iidfnd,['dts']],[float(ik[0]) for ik in df.loc[iidfnd,['Chlorophyll_ug']].values],'.',color='lightgray',alpha=.5)\n",
" pp[ii]=dict()\n",
" for idir in dirnames:\n",
" f=ff[idir][pl]\n",
" for var in ('flagellates',):#,'flagellates','ciliates'):\n",
" try:\n",
" pp[ii][var]=ax[ii].plot(times[idir],2.0*np.sum(f.variables[var][:,:3,0,0]*f.variables['e3t'][:,:3,0,0],1)/np.sum(f.variables['e3t'][:,:3,0,0],1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.5,label=idir+' '+var)\n",
" except:\n",
" pp[ii][var]=ax[ii].plot(times[idir],2.0*np.sum(f.variables[var][:,:3,0,0]*np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1)/np.sum(np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.5,label=idir+' '+var)\n",
" #try:\n",
" # pp[ii]['tot']=ax[ii].plot(times[idir],2.0*np.sum((f.variables['diatoms'][:,:3,0,0]+f.variables['flagellates'][:,:3,0,0]+f.variables['ciliates'][:,:3,0,0])*f.variables['e3t'][:,:3,0,0],1)/np.sum(f.variables['e3t'][:,:3,0,0],1),\n",
" # linestyle=lsty[idir],color='k',alpha=.6,label=idir+' '+'total')\n",
" #except:\n",
" # pp[ii]['tot']=ax[ii].plot(times[idir],2.0*np.sum((f.variables['diatoms'][:,:3,0,0]+f.variables['flagellates'][:,:3,0,0]+f.variables['ciliates'][:,:3,0,0])*np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1)/np.sum(np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1),\n",
" # linestyle=lsty[idir],color='k',alpha=.6,label=idir+' '+'total')\n",
" ax[ii].set_title(pl)\n",
" ax[ii].set_ylabel('chl ($\\mu$g/L)')\n",
" ax[ii].xaxis.set_major_formatter(yearsFmt)\n",
" ax[ii].set_xlim(xl)\n",
" ax[ii].set_ylim(0,40)\n",
" if ii==0:\n",
" ax[0].legend(bbox_to_anchor=(1,1),fontsize=12)\n",
"\n",
"plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"yearsFmt = mdates.DateFormatter('%b %d')\n",
"fig,ax=plt.subplots(len(plist),1,figsize=(16,3*(len(plist))))\n",
"pp=dict()\n",
"for ii in range(0,len(plist)):\n",
" pl=plist[ii]\n",
" if pl=='S3':\n",
" ax[ii].plot(df.loc[iidfnd,['dts']],[float(ik[0]) for ik in df.loc[iidfnd,['Chlorophyll_ug']].values],'.',color='lightgray',alpha=.5)\n",
" pp[ii]=dict()\n",
" for idir in dirnames:\n",
" f=ff[idir][pl]\n",
" for var in ('ciliates',):#,'flagellates','ciliates'):\n",
" try:\n",
" pp[ii][var]=ax[ii].plot(times[idir],2.0*np.sum(f.variables[var][:,:3,0,0]*f.variables['e3t'][:,:3,0,0],1)/np.sum(f.variables['e3t'][:,:3,0,0],1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.5,label=idir+' '+var)\n",
" except:\n",
" pp[ii][var]=ax[ii].plot(times[idir],2.0*np.sum(f.variables[var][:,:3,0,0]*np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1)/np.sum(np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.5,label=idir+' '+var)\n",
" #try:\n",
" # pp[ii]['tot']=ax[ii].plot(times[idir],2.0*np.sum((f.variables['diatoms'][:,:3,0,0]+f.variables['flagellates'][:,:3,0,0]+f.variables['ciliates'][:,:3,0,0])*f.variables['e3t'][:,:3,0,0],1)/np.sum(f.variables['e3t'][:,:3,0,0],1),\n",
" # linestyle=lsty[idir],color='k',alpha=.6,label=idir+' '+'total')\n",
" #except:\n",
" # pp[ii]['tot']=ax[ii].plot(times[idir],2.0*np.sum((f.variables['diatoms'][:,:3,0,0]+f.variables['flagellates'][:,:3,0,0]+f.variables['ciliates'][:,:3,0,0])*np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1)/np.sum(np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1),\n",
" # linestyle=lsty[idir],color='k',alpha=.6,label=idir+' '+'total')\n",
" ax[ii].set_title(pl)\n",
" ax[ii].set_ylabel('chl ($\\mu$g/L)')\n",
" ax[ii].xaxis.set_major_formatter(yearsFmt)\n",
" ax[ii].set_xlim(xl)\n",
" ax[ii].set_ylim(0,40)\n",
" if ii==0:\n",
" ax[0].legend(bbox_to_anchor=(1,1),fontsize=12)\n",
"\n",
"plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# one at a time\n",
"yearsFmt = mdates.DateFormatter('%b %d')\n",
"for idir in dirnames:\n",
" fig,ax=plt.subplots(len(plist),1,figsize=(16,3*(len(plist))))\n",
" pp=dict()\n",
" for ii in range(0,len(plist)):\n",
" pl=plist[ii]\n",
" if pl=='S3':\n",
" ax[ii].plot(df.loc[iidfnd,['dts']],[(1/2.0)*float(ik[0]) for ik in df.loc[iidfnd,['Chlorophyll_ug']].values],'.',color='lightgray',alpha=.5)\n",
" pp[ii]=dict()\n",
" f=ff[idir][pl]\n",
" for var in ('diatoms','flagellates','ciliates'):\n",
" try:\n",
" pp[ii][var]=ax[ii].plot(times[idir],np.sum(f.variables[var][:,:3,0,0]*f.variables['e3t'][:,:3,0,0],1)/np.sum(f.variables['e3t'][:,:3,0,0],1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.6,label=idir+' '+var)\n",
" except:\n",
" pp[ii][var]=ax[ii].plot(times[idir],np.sum(f.variables[var][:,:3,0,0]*np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1)/np.sum(np.tile(e3t_0[:,:3,0,0],(len(times[idir]),1)),1),\n",
" linestyle=lsty[idir],color=lcol[idir][var],alpha=.6,label=idir+' '+var)\n",
" ax[ii].set_title(pl)\n",
" ax[ii].set_ylabel('N (mmol m$^{-2}$)')\n",
" ax[ii].xaxis.set_major_formatter(yearsFmt)\n",
" ax[ii].set_xlim(xl)\n",
" if ii==0:\n",
" ax[0].legend(bbox_to_anchor=(1.25,1),fontsize=12)\n",
"\n",
" fig.suptitle(idir)\n",
" plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f.variables[var][:,:3,0,0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python (py38)",
"language": "python",
"name": "py38"
},
"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.8.1"
}
},
"nbformat": 4,
"nbformat_minor": 1
}