{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%% Project Setup\n"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import jupyter\n",
    "import sys\n",
    "import pandas as pd\n",
    "import geopandas as gp\n",
    "from scipy.signal import argrelextrema\n",
    "import numpy as np\n",
    "import math\n",
    "from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar, AnchoredDirectionArrows\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.font_manager as fm\n",
    "import matplotlib as mpl\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "import contextily as ctx\n",
    "import cmocean.cm as cmo\n",
    "from shapely.geometry import Point, LineString\n",
    "\n",
    "from xarray.backends import NetCDF4DataStore\n",
    "import xarray as xr\n",
    "\n",
    "from datetime import datetime, timedelta\n",
    "from netCDF4 import num2date\n",
    "from metpy.units import units\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "from metpy.plots import ctables\n",
    "\n",
    "from siphon.catalog import TDSCatalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [],
   "source": [
    "def find_time_var(var, time_basename='time'):\n",
    "    for coord_name in var.coords:\n",
    "        if coord_name.startswith(time_basename):\n",
    "            return var.coords[coord_name]\n",
    "    raise ValueError('No time variable found for ' + var.name)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%  Helper function for finding proper time variable\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\arey\\Anaconda3\\envs\\BairdBase\\lib\\site-packages\\openpyxl\\styles\\stylesheet.py:221: UserWarning: Workbook contains no default style, apply openpyxl's default\n",
      "  warn(\"Workbook contains no default style, apply openpyxl's default\")\n"
     ]
    }
   ],
   "source": [
    "importPaths = ['C:/Users/arey/files/Projects/West Coast/Observations/Great House/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Observations/Greensleeves/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Observations/Old Queens Fort/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Updated Observations/Great House/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Updated Observations/Greensleeves/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Updated Observations/Old Queens Fort/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Observations3/Great House/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Observations3/Greensleeves/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Observations3/Old Queens Fort/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Monitoring_Nov/Great House/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Monitoring_Nov/Greensleeves/',\n",
    "               'C:/Users/arey/files/Projects/West Coast/Monitoring_Nov/Old Queens Fort/']\n",
    "\n",
    "siteNames   = ['Great House',\n",
    "               'Greensleeves',\n",
    "               'Old Queens Fort',\n",
    "               'Great House',\n",
    "               'Greensleeves',\n",
    "               'Old Queens Fort',\n",
    "               'Great House',\n",
    "               'Greensleeves',\n",
    "               'Old Queens Fort',\n",
    "               'Great House',\n",
    "               'Greensleeves',\n",
    "               'Old Queens Fort']\n",
    "\n",
    "timeLabels=   ['Before Construction',\n",
    "               'Before Construction',\n",
    "               'Before Construction',\n",
    "               'During Construction',\n",
    "               'During Construction',\n",
    "               'During Construction',\n",
    "               'After Construction',\n",
    "               'After Construction',\n",
    "               'After Construction',\n",
    "               'November Monitoring',\n",
    "               'November Monitoring',\n",
    "               'November Monitoring']\n",
    "\n",
    "wave_bts_file = None\n",
    "\n",
    "for s in range(11, 12):\n",
    "    ## Define master import path\n",
    "    importPath  = importPaths[s]\n",
    "    siteName    = siteNames[s]\n",
    "    timeLabel   = timeLabels[s]\n",
    "    importFiles = os.listdir(importPath)\n",
    "\n",
    "    # Initialize import variables\n",
    "    RBR_File = None\n",
    "    JFE_File = None\n",
    "    GPS_File = None\n",
    "\n",
    "    for i in range(0, len(importFiles)):\n",
    "        if '.xlsx' in importFiles[i] and 'Summary' not in importFiles[i]:\n",
    "            RBR_File = importFiles[i]\n",
    "        elif '_A.csv' in importFiles[i] and 'Summary' not in importFiles[i]:\n",
    "            JFE_File = importFiles[i]\n",
    "        elif '.csv' in importFiles[i] and 'Summary' not in importFiles[i]:\n",
    "            GPS_File = importFiles[i]\n",
    "\n",
    "    #%% RBR Import Data\n",
    "    if RBR_File is not None:\n",
    "        RBR_Obs = pd.read_excel(importPath + RBR_File,\n",
    "                             sheet_name='Data', skiprows=0, header=1)\n",
    "\n",
    "    #%% JFE Import Data\n",
    "    if JFE_File is not None:\n",
    "        JFE_Obs = pd.read_csv(importPath + JFE_File, skiprows=30)\n",
    "\n",
    "    #%% GPS Import Data\n",
    "    GPS = pd.read_csv(importPath + GPS_File,\n",
    "                        header=None, names=['Index', 'Date1', 'Time1', 'Date2', 'Time2', 'Northing', 'North', 'Easting', 'East', 'Var1', 'Var2'])\n",
    "    #convert GPS data to geodataframe\n",
    "    GPS_gdf = gp.GeoDataFrame(GPS, geometry=gp.points_from_xy(-GPS.Easting, GPS.Northing, crs=\"EPSG:4326\"))\n",
    "\n",
    "    GPS_gdf['DateTime'] = pd.to_datetime(GPS_gdf['Date2'].astype(str) + ' ' + GPS_gdf['Time2'].astype(str))\n",
    "\n",
    "    GPS_gdf.set_index('DateTime', inplace=True)\n",
    "\n",
    "    # Convert to UTM\n",
    "    GPS_gdf.geometry = GPS_gdf.geometry.to_crs(\"EPSG:32621\")\n",
    "\n",
    "    # Read in site shapefile\n",
    "    siteShp = gp.read_file('C:/Users/arey/files/Projects/West Coast/SitePolygons.shp')\n",
    "    siteShp.geometry = siteShp.geometry.to_crs(\"EPSG:32621\")\n",
    "\n",
    "\n",
    "    #%% Merge GPS to RBR\n",
    "    # Process RBR into datetime\n",
    "    if RBR_File is not None:\n",
    "        RBR_Obs['DateTime'] = pd.to_datetime(RBR_Obs['Time'])\n",
    "        RBR_Obs.set_index('DateTime', inplace=True)\n",
    "\n",
    "        # Merge with GPS as dataframe\n",
    "        RBR_Obs_geo = pd.merge_asof(RBR_Obs, GPS_gdf,\n",
    "                                     left_index=True, right_index=True,  direction='nearest', tolerance=pd.Timedelta('300s'))\n",
    "        RBR_Obs_geo = gp.GeoDataFrame(RBR_Obs_geo, geometry=RBR_Obs_geo.geometry, crs=\"EPSG:32621\")\n",
    "\n",
    "\n",
    "\n",
    "    #%% Merge GPS to JFE\n",
    "    if JFE_File is not None:\n",
    "        # Process JFE into datetime\n",
    "        JFE_Obs['DateTime'] = pd.to_datetime(JFE_Obs['Date'])\n",
    "        JFE_Obs.set_index('DateTime', inplace=True)\n",
    "\n",
    "        # Merge with GPS as dataframe\n",
    "        JFE_Obs_geo = pd.merge_asof(JFE_Obs, GPS_gdf, left_index=True, right_index=True,  direction='nearest', tolerance=pd.Timedelta('60s'))\n",
    "        JFE_Obs_geo = gp.GeoDataFrame(JFE_Obs_geo, geometry=JFE_Obs_geo.geometry, crs=\"EPSG:32621\")\n",
    "\n",
    "    #%% Find and setup key plotting details\n",
    "\n",
    "    # Shaded Water\n",
    "    mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw'\n",
    "    # Sat water\n",
    "    # mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckekcw3pn08am19qmqbhtq8sb/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw'\n",
    "\n",
    "    if siteName == 'Great House':\n",
    "        axXlim  = (213210.7529575412, 213562.64172686986)\n",
    "        axYlim  = (1464769.2243017585, 1465135.2219089477)\n",
    "        GFS_Lon = -59.6441\n",
    "        GFS_Lat = 13.2372\n",
    "        RBR_Obs_geo['inArea'] = RBR_Obs_geo.within(siteShp.iloc[2,1])\n",
    "    elif siteName == 'Greensleeves':\n",
    "        axXlim  = (213269.99233348924, 213648.1643157148)\n",
    "        # axYlim  = (1463378.1020314451, 1463843.5442048472)\n",
    "        axYlim  = (1463378.1020314451, 1463950.5442048472)\n",
    "        GFS_Lon = -59.6428\n",
    "        GFS_Lat = 13.2289\n",
    "        RBR_Obs_geo['inArea'] = RBR_Obs_geo.within(siteShp.iloc[1,1])\n",
    "    elif siteName == 'Old Queens Fort':\n",
    "        axXlim  = (213368.59866770002, 213745.6997016811)\n",
    "        axYlim  = (1460192.707288096, 1460672.371780407)\n",
    "        GFS_Lon = -59.6419\n",
    "        GFS_Lat = 13.1960\n",
    "        RBR_Obs_geo['inArea'] = RBR_Obs_geo.within(siteShp.iloc[0,1])\n",
    "\n",
    "\n",
    "    # Set min and max times using conductivity\n",
    "    if JFE_File is None:\n",
    "        if RBR_Obs_geo['inArea'].any():\n",
    "            # First and last times from area in shapefile\n",
    "            minTime   = RBR_Obs_geo[RBR_Obs_geo['inArea']==True].iloc[20, 0]\n",
    "            maxTime   = RBR_Obs_geo[RBR_Obs_geo['inArea']==True].iloc[-20, 0]\n",
    "        else:\n",
    "            # First and last times if no GPS data\n",
    "            minTime   = RBR_Obs_geo.iloc[20, 0]\n",
    "            maxTime   = RBR_Obs_geo.iloc[-20, 0]\n",
    "    else:\n",
    "        minTime   = (RBR_Obs['Conductivity ']>5).idxmax()\n",
    "        minTime   = minTime + timedelta(seconds=30)\n",
    "        maxTime   = ((RBR_Obs['Conductivity ']<5) & (RBR_Obs['Time']>minTime)).idxmax()\n",
    "        maxTime   = maxTime - timedelta(seconds=30)\n",
    "        obsPeriod = maxTime-minTime\n",
    "\n",
    "        if (obsPeriod.seconds<180) | (maxTime<minTime):\n",
    "            minTime   = ((RBR_Obs['Conductivity ']>5) & (RBR_Obs['Time']>(minTime+timedelta(seconds=180)))).idxmax()\n",
    "            minTime   = minTime + timedelta(seconds=30)\n",
    "            maxTime   = ((RBR_Obs['Conductivity ']<5) & (RBR_Obs['Time']>minTime)).idxmax()\n",
    "            maxTime   = maxTime - timedelta(seconds=30)\n",
    "\n",
    "    metDate  = minTime - timedelta(\n",
    "                                hours = minTime.hour % 6,\n",
    "                                minutes=minTime.minute,\n",
    "                                seconds=minTime.second,\n",
    "                                microseconds=minTime.microsecond)\n",
    "\n",
    "\n",
    "    #%% GFS Met Import\n",
    "    var = ['Temperature_surface', 'Wind_speed_gust_surface',\n",
    "           'u-component_of_wind_height_above_ground', 'v-component_of_wind_height_above_ground']\n",
    "    var_precp = ['Total_precipitation_surface_6_Hour_Accumulation']\n",
    "\n",
    "    temp_1d = []\n",
    "    gust_1d = []\n",
    "    wndu_1d = []\n",
    "    wndv_1d = []\n",
    "    prep_1d = []\n",
    "    time_1d = []\n",
    "\n",
    "    # Set times to download\n",
    "    startdate  = metDate - timedelta(hours=18)\n",
    "    enddate    = metDate + timedelta(hours=6)\n",
    "    date_list = pd.date_range(startdate, enddate, freq='6H')\n",
    "\n",
    "    # Loop through dates\n",
    "    for date in date_list:\n",
    "        # Base URL for 0.5 degree GFS data\n",
    "        best_gfs = TDSCatalog('https://www.ncei.noaa.gov/thredds/catalog/model-gfs-g4-anl-files/' +\n",
    "                            date.strftime('%Y%m') + '/' + date.strftime('%Y%m%d') + '/' + 'catalog.xml')\n",
    "\n",
    "        # Generate URLs for specific grib file\n",
    "        best_ds         = best_gfs.datasets['gfs_4_'+date.strftime('%Y%m%d_%H%M')+'_000.grb2']\n",
    "        best_ds_precp   = best_gfs.datasets['gfs_4_'+date.strftime('%Y%m%d_%H%M')+'_006.grb2']\n",
    "\n",
    "        # Format the query parameters\n",
    "        ncss  = best_ds.subset()\n",
    "        query = ncss.query()\n",
    "\n",
    "        ncss_precp  = best_ds_precp.subset()\n",
    "        query_precp = ncss_precp.query()\n",
    "\n",
    "        # Extract data from specific point\n",
    "        query.lonlat_point(GFS_Lon, GFS_Lat).time(date)\n",
    "        query.accept('netcdf')\n",
    "        query.variables(var[0], var[1], var[2], var[3])\n",
    "        query.vertical_level(10)\n",
    "\n",
    "        data = ncss.get_data(query)\n",
    "        data = xr.open_dataset(NetCDF4DataStore(data), drop_variables='height_above_ground4')\n",
    "\n",
    "        query_precp.lonlat_point(GFS_Lon, GFS_Lat).time(date + timedelta(hours=6))\n",
    "        query_precp.accept('netcdf')\n",
    "        query_precp.variables(var_precp[0])\n",
    "\n",
    "\n",
    "        data_precp = ncss_precp.get_data(query_precp)\n",
    "        data_precp = xr.open_dataset(NetCDF4DataStore(data_precp))\n",
    "\n",
    "        temp_3d = data[var[0]]\n",
    "        gust_3d = data[var[1]]\n",
    "        wndu_3d = data[var[2]]\n",
    "        wndv_3d = data[var[3]]\n",
    "        prep_3d = data_precp[var_precp[0]]\n",
    "\n",
    "        # Read the individual point (with units) and append to the list\n",
    "        temp_1d.append(temp_3d.metpy.unit_array.squeeze())\n",
    "        gust_1d.append(gust_3d.metpy.unit_array.squeeze())\n",
    "        wndu_1d.append(wndu_3d.metpy.unit_array.squeeze())\n",
    "        wndv_1d.append(wndv_3d.metpy.unit_array.squeeze())\n",
    "        prep_1d.append(prep_3d.metpy.unit_array.squeeze())\n",
    "        time_1d.append(find_time_var(temp_3d))\n",
    "\n",
    "    #%% Process Met Data\n",
    "    # 24h Precipitation Total\n",
    "    # Time weighted average of last two points for everything else\n",
    "\n",
    "    met_prep = prep_1d[0] + prep_1d[1] + prep_1d[2] + prep_1d[3]\n",
    "\n",
    "    timeWeight1 = minTime-metDate\n",
    "    timeWeight2 = (metDate+timedelta(hours=6))-minTime\n",
    "\n",
    "    timeWeight1 = timeWeight1.seconds/21600\n",
    "    timeWeight2 = timeWeight2.seconds/21600\n",
    "\n",
    "    met_gust = gust_1d[3] * timeWeight1 + gust_1d[4] * timeWeight2\n",
    "    met_temp = temp_1d[3] * timeWeight1 + temp_1d[4] * timeWeight2\n",
    "    met_wind =  math.sqrt((wndv_1d[3].m.item(0) * timeWeight1 + wndv_1d[4].m.item(0)* timeWeight2) ** 2 +\n",
    "                          (wndu_1d[3].m.item(0) * timeWeight1 + wndu_1d[4].m.item(0)* timeWeight2) **2 )\n",
    "    met_wdir = math.degrees(math.atan2(\n",
    "                wndv_1d[3].m.item(0) * timeWeight1 + wndv_1d[4].m.item(0)* timeWeight2,\n",
    "                wndu_1d[3].m.item(0) * timeWeight1 + wndu_1d[4].m.item(0)* timeWeight2)) % 360\n",
    "\n",
    "    #%% Read in wave conditions from BTS file\n",
    "    if wave_bts_file is not None:\n",
    "        if siteName == 'Great House':\n",
    "            wave_bts = pd.read_csv(wave_bts_file,\n",
    "                        names=['date', 'time', 'HM0', 'TP', 'TM', 'MWD', 'DPK', 'HSWL', 'TSWL', 'DPSWL', 'HSEA', 'TSEA', 'DPSEA'],\n",
    "                        header=0, skiprows=22, delim_whitespace=True)\n",
    "            wave_bts['datetime'] = pd.to_datetime(wave_bts['date'] + ' ' + wave_bts['time'])\n",
    "            wave_bts.set_index('datetime', inplace=True)\n",
    "\n",
    "\n",
    "        met_hmo = wave_bts.iloc[wave_bts.index.get_loc(minTime, method='nearest'), :].HM0\n",
    "        met_tp  = wave_bts.iloc[wave_bts.index.get_loc(minTime, method='nearest'), :].TP\n",
    "        met_mwd = wave_bts.iloc[wave_bts.index.get_loc(minTime, method='nearest'), :].MWD\n",
    "    else:\n",
    "        met_hmo = -999\n",
    "        met_tp  = -999\n",
    "        met_mwd = -999\n",
    "\n",
    "\n",
    "    #%% Plot time series for Geo data\n",
    "    # RBRparam = ['Depth ', 'Salinity ', 'Dissolved O₂ saturation ', 'Temperature ']\n",
    "    # RBRparamName = ['Depth [m]', 'Salinity [PSU]', 'Dissolved O₂ saturation [%]', 'Temperature [degC]']\n",
    "    # RBRparamMin = [0.0, 34.5, 32.0, 29.0]\n",
    "    # RBRparamMax = [1.0, 36.0, 34.5, 31.0]\n",
    "    #\n",
    "    # JFEparam     = ['Turb. -M[FTU]', 'Chl-Flu.[ppb]', 'Chl-a[ug/l]']\n",
    "    # JFEparamName = ['Turbidity [FTU]', 'Chl-Flu. [ppb]', 'Chl-a [ug/l]']\n",
    "    # JFEparamMin  = [0.0, 0.0, 0.0]\n",
    "    # JFEparamMax  = [20.0, 1.0, 1.0]\n",
    "\n",
    "    if JFE_File is None:\n",
    "        fig, axesTMP = plt.subplots(nrows=1, ncols=1, figsize=(19, 5), constrained_layout=True)\n",
    "\n",
    "        RBRparam        = ['Turbidity ']\n",
    "        RBRparamName    = ['Turbidity [NTU]']\n",
    "        RBRparmCmap     = [cmo.turbid]\n",
    "        RBRparamMin     = [0.0]\n",
    "        RBRparamMax     = [60.0]\n",
    "\n",
    "        dataTable = np.zeros([9, 4])\n",
    "        axes = []\n",
    "        axes.append(axesTMP)\n",
    "    else:\n",
    "        fig, axes = plt.subplots(nrows=7, ncols=1, figsize=(19, 25), constrained_layout=True)\n",
    "        dataTable = np.zeros([15, 4])\n",
    "\n",
    "        RBRparam        = ['Depth ', 'Salinity ', 'Dissolved O₂ saturation ', 'Temperature ']\n",
    "        RBRparamName    = ['Depth [m]', 'Salinity [PSU]', 'Dissolved O₂ saturation [%]', 'Temperature [degC]']\n",
    "        RBRparmCmap     = [cmo.deep, 'cividis', cmo.dense, cmo.thermal]\n",
    "        RBRparamMin     = [0.0, 35.0, 32.5, 29.0]\n",
    "        RBRparamMax     = [1.0, 36.0, 33.5, 31.0]\n",
    "\n",
    "    fig.patch.set_facecolor('white')\n",
    "    # fig.tight_layout(pad=1.05)\n",
    "\n",
    "    fontprops = fm.FontProperties(size=25)\n",
    "\n",
    "    dataTable[0, 0] = met_temp.m.item(0)-272.15\n",
    "    dataTable[1, 0] = met_wind\n",
    "    dataTable[2, 0] = met_wdir\n",
    "    dataTable[3, 0] = met_prep.m.item(0)\n",
    "    dataTable[4, 0] = met_hmo\n",
    "    dataTable[5, 0] = met_tp\n",
    "    dataTable[6, 0] = met_mwd\n",
    "\n",
    "    for paramIDX, param in enumerate(RBRparam):\n",
    "        RBR_Obs_geo.loc[minTime:maxTime, param].plot(\n",
    "                        ax=axes[paramIDX], label='1 Second Observations', color='lightgrey') # Note the space in the col name\n",
    "\n",
    "        # Create mask for RBR data based on time and parameter minimum\n",
    "        RBR_mask = ((RBR_Obs_geo.index>minTime) &\n",
    "                         (RBR_Obs_geo.index<maxTime) &\n",
    "                         (RBR_Obs_geo[param]>RBRparamMin[paramIDX]))\n",
    "        RBR_smoothed = RBR_Obs_geo.loc[RBR_mask, param].rolling(\n",
    "                60, win_type='nuttall',center=True).mean()\n",
    "\n",
    "        RBR_smoothed.plot(\n",
    "                ax=axes[paramIDX], label='1 Minute Average', color='black', linewidth=3) # Note the space in the col name\n",
    "\n",
    "        # Find the local maximums\n",
    "        ilocs_max = argrelextrema(\n",
    "                            RBR_smoothed.values,\n",
    "                            np.greater_equal, order=60)[0]\n",
    "\n",
    "        # Add points at beginning and end\n",
    "        ilocs_max       = np.insert(ilocs_max, 0, 10)\n",
    "        ilocs_max_pts   = RBR_smoothed.iloc[ilocs_max].index.values\n",
    "        # ilocs_max[-1]   = len(RBR_smoothed.values)-10\n",
    "\n",
    "        plt.scatter(RBR_smoothed.iloc[\n",
    "                         ilocs_max].index, np.ones(len(ilocs_max)) * 30, 75,\n",
    "                         color='blue')\n",
    "\n",
    "        # Add labels\n",
    "        for a in range(0, len(ilocs_max)):\n",
    "            plt.annotate(str(a+1), (ilocs_max_pts[a], 30), fontsize=30)\n",
    "\n",
    "\n",
    "        dataTable[paramIDX+7, 0] = RBR_Obs_geo.loc[minTime:maxTime, param].mean(skipna=True)\n",
    "        dataTable[paramIDX+7, 1] = RBR_Obs_geo.loc[minTime:maxTime, param].std(skipna=True)\n",
    "        dataTable[paramIDX+7, 2] = max(RBR_Obs_geo.loc[minTime:maxTime, param].min(skipna=True), 0)\n",
    "        dataTable[paramIDX+7, 3] = RBR_Obs_geo.loc[minTime:maxTime, param].max(skipna=True)\n",
    "\n",
    "        axes[paramIDX].set_ylabel(RBRparamName[paramIDX])\n",
    "        axes[paramIDX].set_title('RBR: ' + RBRparamName[paramIDX])\n",
    "        axes[paramIDX].set_xlabel('')\n",
    "        axes[paramIDX].set_ylim(RBRparamMin[paramIDX], RBRparamMax[paramIDX])\n",
    "        axes[paramIDX].legend(loc='upper right')\n",
    "\n",
    "    if JFE_File is not None:\n",
    "        for paramIDX, param in enumerate(JFEparam):\n",
    "            JFE_Obs_geo.loc[minTime:maxTime, param].plot(\n",
    "                            ax=axes[paramIDX+4], label='15 Second Observations') # Note the space in the col name\n",
    "            # JFEDats[siteIDX].loc[minTime:maxTime, param].rolling(\n",
    "            #             300, win_type='gaussian',center=True).mean(std=5).plot(\n",
    "            #             ax=axes[paramIDX+4], label='15 Minute Average') # Note the space in the col name\n",
    "\n",
    "            dataTable[paramIDX+4+7, 0] = JFE_Obs_geo.loc[minTime:maxTime, param].mean(skipna=True)\n",
    "            dataTable[paramIDX+4+7, 1] = JFE_Obs_geo.loc[minTime:maxTime, param].std(skipna=True)\n",
    "            dataTable[paramIDX+4+7, 2] = max(JFE_Obs_geo.loc[minTime:maxTime, param].min(skipna=True), 0)\n",
    "            dataTable[paramIDX+4+7, 3] = JFE_Obs_geo.loc[minTime:maxTime, param].max(skipna=True)\n",
    "\n",
    "            axes[paramIDX+4].set_ylabel(JFEparamName[paramIDX])\n",
    "            axes[paramIDX+4].set_title('JFE: ' + JFEparamName[paramIDX])\n",
    "            axes[paramIDX+4].set_xlabel('')\n",
    "            # Set turbidity to log\n",
    "            # if paramIDX==2:\n",
    "            #     axes[6].set_yscale('log')\n",
    "            axes[paramIDX+4].set_ylim(JFEparamMin[paramIDX], JFEparamMax[paramIDX])\n",
    "            axes[paramIDX+4].legend(loc='upper right')\n",
    "\n",
    "        axes[paramIDX+4].set_xlabel(minTime.strftime('%B %#d, %Y'))\n",
    "    else:\n",
    "        axes[paramIDX].set_xlabel(minTime.strftime('%B %#d, %Y'))\n",
    "\n",
    "    dfOut = pd.DataFrame (dataTable[:, :])\n",
    "    # Change the column names\n",
    "    dfOut.columns =['Mean', 'Standard Deviation', 'Min', 'Max']\n",
    "    # Change the row indexes\n",
    "    if JFE_File is not None:\n",
    "        dfOut.iloc[14, 0] = minTime.strftime('%B %#d, %Y')\n",
    "        dfOut.index = ['Air Temperature [degC]', 'Wind Speed [m/s]', 'Wind Direction [deg]', '24h Precipitation [mm]',\n",
    "                   'Significant Wave Height Offshore [m]' ,'Peak Wave Period Offshore [s]',\n",
    "                   'Mean Wave Direction Offshore [deg]', 'Depth [m]', 'Salinity [PSU]',\n",
    "                   'Dissolved O₂ saturation [%]', 'Temperature [degC]',\n",
    "                   'Turbidity [FTU]', 'Chl-Flu. [ppb]', 'Chl-a [ug/l]', 'Observation Date']\n",
    "    else:\n",
    "        dfOut.iloc[8, 0] = minTime.strftime('%B %#d, %Y')\n",
    "        dfOut.index = ['Air Temperature [degC]', 'Wind Speed [m/s]', 'Wind Direction [deg]', '24h Precipitation [mm]',\n",
    "                   'Significant Wave Height Offshore [m]' ,'Peak Wave Period Offshore [s]',\n",
    "                   'Mean Wave Direction Offshore [deg]',\n",
    "                   'Turbidity [NTU]', 'Observation Date']\n",
    "\n",
    "    fig.suptitle(siteName + ', ' +minTime.strftime('%B %#d, %Y') + ' (' + timeLabel + ')', fontsize=30)\n",
    "\n",
    "    plt.show()\n",
    "    dfOut.to_excel(importPath + '/Summary_Stats_' + siteName + '.xlsx')\n",
    "    dfOut.to_csv(importPath + '/Summary_Stats_' + siteName + '.csv')\n",
    "\n",
    "    fig.savefig(importPath + '/Figures/SummaryTimeSeries_' + siteName + '.pdf',\n",
    "                 bbox_inches='tight')\n",
    "    fig.savefig(importPath + '/Figures/SummaryTimeSeries_' + siteName + '.png',\n",
    "                 bbox_inches='tight', dpi=500)\n",
    "\n",
    "\n",
    "    #%% Plot Points\n",
    "    if JFE_File is None:\n",
    "        # Only Turbidity Data\n",
    "        fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(9, 9), constrained_layout=True)\n",
    "        RBRparam        = ['Turbidity ']\n",
    "        RBRparamName    = ['Turbidity [NTU]']\n",
    "        RBRparmCmap     = [cmo.turbid]\n",
    "        RBRparamMin     = [0.0]\n",
    "        RBRparamMax     = [60.0]\n",
    "\n",
    "        ax = []\n",
    "        ax.append(axes)\n",
    "    else:\n",
    "        fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(19, 25), constrained_layout=True)\n",
    "        # List plotting parameters\n",
    "        RBRparam        = ['Depth ', 'Salinity ', 'Dissolved O₂ saturation ', 'Temperature ']\n",
    "        RBRparamName    = ['Depth [m]', 'Salinity [PSU]', 'Dissolved O₂ saturation [%]', 'Temperature [degC]']\n",
    "        RBRparmCmap     = [cmo.deep, 'cividis', cmo.dense, cmo.thermal]\n",
    "        RBRparamMin     = [0.0, 35.0, 32.5, 29.0]\n",
    "        RBRparamMax     = [1.0, 36.0, 33.5, 31.0]\n",
    "\n",
    "        ax = axes.flat\n",
    "\n",
    "    fig.patch.set_facecolor('white')\n",
    "    # fig.tight_layout(pad=1.05)\n",
    "\n",
    "    fontprops = fm.FontProperties(size=25)\n",
    "    x, y, arrow_length = 0.95, 0.93, 0.20\n",
    "    plt.rcParams.update({'font.size': 22})\n",
    "\n",
    "    axXlimTT = (RBR_Obs_geo.loc[minTime:maxTime].geometry.x.min()-100,\n",
    "                     RBR_Obs_geo.loc[minTime:maxTime].geometry.x.max()+100)\n",
    "    axYlimTT = (RBR_Obs_geo.loc[minTime:maxTime].geometry.y.min()-100,\n",
    "                     RBR_Obs_geo.loc[minTime:maxTime].geometry.y.max()+100)\n",
    "\n",
    "    plt.setp(axes, xlim=axXlim, ylim=axYlim)\n",
    "    # plt.setp(axes, xlim=axXlimTT, ylim=axYlimTT)\n",
    "\n",
    "    # Plot the RBR observations\n",
    "    # Salinity\n",
    "    for paramIDX, param in enumerate(RBRparam):\n",
    "        if RBR_File is not None:\n",
    "            RBR_Obs_geo.loc[minTime:maxTime].plot(\n",
    "                                column=param, ax=ax[paramIDX], vmin=RBRparamMin[paramIDX], vmax=RBRparamMax[paramIDX],\n",
    "                                legend=True, legend_kwds={'label': RBRparamName[paramIDX]},\n",
    "                                cmap=RBRparmCmap[paramIDX], markersize=20) # Note the space in the col name\n",
    "            ctx.add_basemap(ax[paramIDX], source=mapbox, crs='EPSG:32621')\n",
    "\n",
    "            ax[paramIDX].set_title(RBRparamName[paramIDX])\n",
    "            # ax[paramIDX].set_ylabel('UTM 21N [m]')\n",
    "            # ax[paramIDX].set_xlabel('UTM 21N [m]')\n",
    "            ax[paramIDX].locator_params(axis='y', nbins=3)\n",
    "            ax[paramIDX].ticklabel_format(useOffset=False, style='plain', axis='both')\n",
    "\n",
    "            ax[paramIDX].get_xaxis().set_ticks([])\n",
    "            ax[paramIDX].get_yaxis().set_ticks([])\n",
    "\n",
    "            #Add scale-bar\n",
    "            scalebar = AnchoredSizeBar(ax[paramIDX].transData,\n",
    "                                   100, '100 m', 'lower right', pad=0.5, size_vertical=10, fontproperties=fontprops)\n",
    "            ax[paramIDX].add_artist(scalebar)\n",
    "            ax[paramIDX].annotate('N', xy=(x, y), xytext=(x, y-arrow_length),\n",
    "                    arrowprops=dict(facecolor='black', width=6, headwidth=30),\n",
    "                    ha='center', va='center', fontsize=35,\n",
    "                    xycoords=ax[paramIDX].transAxes)\n",
    "\n",
    "    # Plot Plot JFE Points\n",
    "\n",
    "    # JFEparam     = ['Temp.[deg C]', 'Turb. -M[FTU]', 'Chl-a[ug/l]', 'Chl-Flu.[ppb]']\n",
    "    # JFEparamName = ['Temperature [degC]', 'Turbidity [FTU]', 'Chl-a [ug/l]', 'Chl-Flu. [ppb]']\n",
    "    # JFEparamCmp  = [cmo.thermal, cmo.turbid, cmo.algae, cmo.algae]\n",
    "    JFEparam     = ['Turb. -M[FTU]', 'Chl-Flu.[ppb]']\n",
    "    JFEparamName = ['Turbidity [FTU]', 'Chl-Flu. [ppb]']\n",
    "    JFEparamCmp  = [cmo.turbid,  cmo.algae]\n",
    "    JFEparamMin = [0.0, 0.0]\n",
    "    JFEparamMax = [1.0, 1.0]\n",
    "    if JFE_File is not None:\n",
    "        for paramIDX, param in enumerate(JFEparam):\n",
    "            if JFE_File is not None:\n",
    "                JFE_Obs_geo.loc[minTime:maxTime].plot(\n",
    "                                    column=param, ax=ax[paramIDX+4], vmin=JFEparamMin[paramIDX], vmax=JFEparamMax[paramIDX],\n",
    "                                    legend=True, legend_kwds={'label': JFEparamName[paramIDX]},\n",
    "                                    cmap=JFEparamCmp[paramIDX], markersize=20) # Note the space in the col name\n",
    "                ctx.add_basemap(ax[paramIDX+4], source=mapbox, crs='EPSG:32621')\n",
    "\n",
    "                ax[paramIDX+4].set_title(JFEparamName[paramIDX])\n",
    "                # ax[paramIDX+4].set_ylabel('UTM 21N [m]')\n",
    "                # ax[paramIDX+4].set_xlabel('UTM 21N [m]')\n",
    "                ax[paramIDX+4].locator_params(axis='y', nbins=3)\n",
    "                ax[paramIDX+4].ticklabel_format(useOffset=False, style='plain', axis='both')\n",
    "                ax[paramIDX+4].get_xaxis().set_ticks([])\n",
    "                ax[paramIDX+4].get_yaxis().set_ticks([])\n",
    "\n",
    "                #Add scale-bar\n",
    "                scalebar = AnchoredSizeBar(ax[paramIDX+4].transData,\n",
    "                                       100, '100 m', 'lower right', pad=0.5, size_vertical=10, fontproperties=fontprops)\n",
    "                ax[paramIDX+4].add_artist(scalebar)\n",
    "                ax[paramIDX+4].annotate('N', xy=(x, y), xytext=(x, y-arrow_length),\n",
    "                        arrowprops=dict(facecolor='black', width=6, headwidth=30),\n",
    "                        ha='center', va='center', fontsize=25,\n",
    "                        xycoords=ax[paramIDX+4].transAxes)\n",
    "\n",
    "    fig.suptitle(siteName + ', ' + minTime.strftime('%B %#d, %Y') + ' (' + timeLabel + ')', fontsize=30)\n",
    "    plt.show()\n",
    "\n",
    "    if not os.path.exists(importPath + '/Figures'):\n",
    "        os.mkdir(importPath + '/Figures')\n",
    "\n",
    "    fig.savefig(importPath + '/Figures/SummaryMap_' + siteName + '.pdf',\n",
    "                 bbox_inches='tight')\n",
    "\n",
    "    fig.savefig(importPath + '/Figures/SummaryMap_' + siteName + '.png',\n",
    "                 bbox_inches='tight', dpi=500)\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
     "is_executing": true
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 2736x3600 with 10 Axes>",
      "image/png": "[encoded data, 355276 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotvars = ['Air Temperature [degC]', 'Wind Speed [m/s]', 'Wind Direction [deg]', '24h Precipitation [mm]',\n",
    "           'Significant Wave Height [m]', 'Salinity [PSU]',\n",
    "           'Dissolved O₂ [%]', 'Temperature [degC]',\n",
    "           'Turbidity [FTU]', 'Chl-Flu. [ppb]']\n",
    "\n",
    "\n",
    "plotIDXs = [0, 3, 6]\n",
    "# plotIDXs = [1, 4, 7]\n",
    "# plotIDXs = [2, 5, 8]\n",
    "plotDates = []\n",
    "plotTable = np.zeros([10, len(plotIDXs)])\n",
    "\n",
    "for s, plotIDX in enumerate(plotIDXs):\n",
    "    ## Define master import path\n",
    "    importPath  = importPaths[plotIDX]\n",
    "    siteName    = siteNames[plotIDX]\n",
    "\n",
    "    obsStatsIN  = pd.read_excel(importPath + '/Summary_Stats.xlsx')\n",
    "\n",
    "    plotTable[:, s] = obsStatsIN.iloc[[0, 1, 2, 3, 4, 8, 9, 10, 11, 12], 1]\n",
    "    plotDates.append(datetime.strptime(obsStatsIN.iloc[14, 1], '%B %d, %Y'))\n",
    "\n",
    "\n",
    "fig, axes = plt.subplots(nrows=5, ncols=2, figsize=(19, 25))\n",
    "fig.patch.set_facecolor('white')\n",
    "fig.tight_layout(pad=3)\n",
    "ax = axes.flat\n",
    "\n",
    "for v in range(0, 10):\n",
    "    ax[v].scatter(plotDates, plotTable[v, :], 250)\n",
    "    ax[v].set_ylabel(plotvars[v])\n",
    "\n",
    "fig.suptitle(siteName, fontsize=35)\n",
    "plt.gcf().autofmt_xdate()\n",
    "plt.gcf().align_ylabels()\n",
    "plt.show()\n",
    "\n",
    "fig.savefig('C:/Users/arey/files/Projects/West Coast/Scatterplots/' + siteName + '.pdf',\n",
    "             bbox_inches='tight')\n",
    "fig.savefig('C:/Users/arey/files/Projects/West Coast/Scatterplots/' + siteName + '.png',\n",
    "             bbox_inches='tight', dpi=500)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Summary Plot\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PyCharm (Mustique)",
   "language": "python",
   "name": "pycharm-78fc648f"
  },
  "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}