{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\"\"\"\n",
    "@author: AJMR\n",
    "\n",
    "Script to compare EFDC results against observations\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%% Import\n"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "import numpy as np\n",
    "import geopandas as gp\n",
    "gp.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'\n",
    "from shapely.ops import nearest_points\n",
    "from shapely.geometry import LineString\n",
    "\n",
    "from shapely.geometry import Point, MultiPoint\n",
    "from shapely.ops import nearest_points\n",
    "\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "\n",
    "from dfm_tools.get_nc import get_netdata, get_ncmodeldata, plot_netmapdata\n",
    "from dfm_tools.get_nc_helpers import get_timesfromnc\n",
    "from dfm_tools.regulargrid import scatter_to_regulargrid\n",
    "from pathlib import Path\n",
    "import netCDF4 as nc\n",
    "\n",
    "import xarray as xr\n",
    "import pytz\n",
    "import pathlib as pl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%% Function to find nearest point\n"
    }
   },
   "outputs": [],
   "source": [
    "def get_nearest_values(row, other_gdf, point_column='geometry', value_column=\"geometry\"):\n",
    "    \"\"\"Find the nearest point and return the corresponding value from specified value column.\"\"\"\n",
    "\n",
    "    # Create an union of the other GeoDataFrame's geometries:\n",
    "    other_points = other_gdf[\"geometry\"].unary_union\n",
    "\n",
    "    # Find the nearest points\n",
    "    nearest_geoms = nearest_points(row[point_column], other_points)\n",
    "\n",
    "    # Get corresponding values from the other df\n",
    "    nearest_data = other_gdf.loc[other_gdf[\"geometry\"] == nearest_geoms[1]]\n",
    "\n",
    "    nearest_value = nearest_data[value_column].get_values()[0]\n",
    "\n",
    "    return nearest_value\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%  Load in Water Level Data\n"
    }
   },
   "outputs": [],
   "source": [
    "# Gauge Locations from map\n",
    "# gdf_gaugeLoc = gp.read_file('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/NTC6.kml')\n",
    "gdf_gaugeLoc = gp.read_file('//srv-ott3.baird.com/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/NTC6.kml')\n",
    "gdf_gaugeLocUSSP = gdf_gaugeLoc.to_crs(\"EPSG:32118\")\n",
    "\n",
    "# All in Feet NAVD88\n",
    "# df_wl_FFG  = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/04 Gauge Data/NCP1_Gauge_Elev_Data_20130521/NC_Gauge_Elev_Data_20130521.xlsx',\n",
    "df_wl_FFG  = pd.read_excel('//srv-ott3.baird.com/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/04 Gauge Data/NCP1_Gauge_Elev_Data_20130521/NC_Gauge_Elev_Data_20130521.xlsx',\n",
    "                               sheet_name='Field_Facility_Gauge')\n",
    "df_wl_FFG.set_index('Date_time', inplace=True)\n",
    "df_wl_FFG.index = df_wl_FFG.index.tz_localize(\n",
    "    pytz.timezone('America/New_York')).tz_convert(pytz.utc) #Convert to UTC\n",
    "\n",
    "gdf_wl_FFG = gp.GeoDataFrame(df_wl_FFG,\n",
    "                             geometry=gp.points_from_xy(np.ones([len(df_wl_FFG), 1]) * gdf_gaugeLoc['geometry'].x[2],\n",
    "                             np.ones([len(df_wl_FFG), 1]) * gdf_gaugeLoc['geometry'].y[2]), crs=\"EPSG:4326\")\n",
    "gdf_wl_FFG.geometry = gdf_wl_FFG.geometry.to_crs(\"EPSG:32118\")\n",
    "\n",
    "\n",
    "df_wl_NGG1 = pd.read_excel('//srv-ott3.baird.com/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/04 Gauge Data/NCP1_Gauge_Elev_Data_20130521/NC_Gauge_Elev_Data_20130521.xlsx',\n",
    "                               sheet_name='National_Grid_Gauge')\n",
    "df_wl_NGG1.set_index('Date_time', inplace=True)\n",
    "df_wl_NGG1.index    = df_wl_NGG1.index.tz_localize(\n",
    "    pytz.timezone('America/New_York'), ambiguous=True).tz_convert(pytz.utc) #Convert to UTC\n",
    "\n",
    "gdf_wl_NGG1 = gp.GeoDataFrame(df_wl_NGG1,\n",
    "                             geometry=gp.points_from_xy(np.ones([len(df_wl_NGG1), 1]) * gdf_gaugeLoc['geometry'].x[3],\n",
    "                             np.ones([len(df_wl_NGG1), 1]) * gdf_gaugeLoc['geometry'].y[3]), crs=\"EPSG:4326\")\n",
    "gdf_wl_NGG1.geometry = gdf_wl_NGG1.geometry.to_crs(\"EPSG:32118\")\n",
    "\n",
    "\n",
    "# Also includes temperature\n",
    "df_wl_NGG2 = pd.read_excel('//srv-ott3.baird.com/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/04 Gauge Data/NCP2_NG_TideGauge_20160113/NCP2_NG_TideGauge_20160113.xlsx',\n",
    "                               sheet_name='Sheet1')\n",
    "gdf_wl_NGG2 = gp.GeoDataFrame(df_wl_NGG2,\n",
    "                             geometry=gp.points_from_xy(np.ones([len(df_wl_NGG2), 1]) * gdf_gaugeLoc['geometry'].x[3],\n",
    "                             np.ones([len(df_wl_NGG2), 1]) * gdf_gaugeLoc['geometry'].y[3]), crs=\"EPSG:4326\")\n",
    "gdf_wl_NGG2.geometry = gdf_wl_NGG2.geometry.to_crs(\"EPSG:32118\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "pycharm": {
     "name": "#%%  Read in EFDC Data to dataframe\n"
    }
   },
   "outputs": [],
   "source": [
    "# efdc6_df = pd.read_csv('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/2012_Mon6/READ_FROM_EFDC_GRAPHICS_OUT_BIN.CSV')\n",
    "# efdc7_df = pd.read_csv('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/2012_Mon7/READ_FROM_EFDC_GRAPHICS_OUT_BIN.CSV')\n",
    "# efdc8_df = pd.read_csv('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/2012_Mon8/READ_FROM_EFDC_GRAPHICS_OUT_BIN.CSV')\n",
    "# efdc9_df = pd.read_csv('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/2012_Mon9/READ_FROM_EFDC_GRAPHICS_OUT_BIN.CSV')\n",
    "\n",
    "efdc12_df = pd.read_csv('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/2012_Mon12/READ_FROM_EFDC_GRAPHICS_OUT_BIN.CSV')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "pycharm": {
     "name": "#%%  Read in EFDC grid\n"
    }
   },
   "outputs": [],
   "source": [
    "efdc_grid_df = pd.read_csv('//Ott-flavius/j/INPUT_FILES/GRID_FILES/NC_ER_lxly_20140129.inp',\n",
    "                           delim_whitespace=True, skiprows=4,\n",
    "                           names=['I', 'J', 'X', 'Y', 'CUE', 'CVE', 'CUN', 'CVN'])\n",
    "\n",
    "efdc_grid_gdf = gp.GeoDataFrame(\n",
    "    efdc_grid_df, geometry=gp.points_from_xy(efdc_grid_df.X, efdc_grid_df.Y), crs=\"EPSG:32118\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "pycharm": {
     "name": "#%%  Merge EFDC outputs with grid locations and add datetimes\n"
    }
   },
   "outputs": [],
   "source": [
    "efdc6_merged = pd.merge(efdc6_df, efdc_grid_df,  how='left', left_on=['I_MOD','J_MOD'],\n",
    "                    right_on = ['I','J']).drop(columns=[\n",
    "                    ' DUMPID', 'END_hr', 'I_MOD', 'J_MOD', 'I_MOD', 'X', 'Y', 'CUE', 'CVE', 'CUN', 'CVN'])\n",
    "efdc6_merged['date'] =  pd.to_datetime([datetime(2012, 6, 1, 00, 00, 00) +\n",
    "                timedelta(hours=h) for h in efdc6_merged['ST_hr']])\n",
    "efdc6_merged.set_index('date', inplace=True)\n",
    "efdc6_merged.index = efdc6_merged.index.tz_localize(\n",
    "    pytz.timezone('EST')).tz_convert(pytz.utc) #Convert to UTC\n",
    "del efdc6_df\n",
    "\n",
    "\n",
    "efdc7_merged = pd.merge(efdc7_df, efdc_grid_df,  how='left', left_on=['I_MOD','J_MOD'],\n",
    "                    right_on = ['I','J']).drop(columns=[\n",
    "                    ' DUMPID', 'END_hr', 'I_MOD', 'J_MOD', 'I_MOD', 'X', 'Y', 'CUE', 'CVE', 'CUN', 'CVN'])\n",
    "efdc7_merged['date'] =  pd.to_datetime([datetime(2012, 7, 1, 00, 00, 00) +\n",
    "                timedelta(hours=h) for h in efdc7_merged['ST_hr']])\n",
    "efdc7_merged.set_index('date', inplace=True)\n",
    "efdc7_merged.index = efdc7_merged.index.tz_localize(\n",
    "    pytz.timezone('EST')).tz_convert(pytz.utc) #Convert to UTC\n",
    "del efdc7_df\n",
    "\n",
    "\n",
    "efdc8_merged = pd.merge(efdc8_df, efdc_grid_df,  how='left', left_on=['I_MOD','J_MOD'],\n",
    "                    right_on = ['I','J']).drop(columns=[\n",
    "                    ' DUMPID', 'END_hr', 'I_MOD', 'J_MOD', 'I_MOD', 'X', 'Y', 'CUE', 'CVE', 'CUN', 'CVN'])\n",
    "efdc8_merged['date'] =  pd.to_datetime([datetime(2012, 8, 1, 00, 00, 00) +\n",
    "                timedelta(hours=h) for h in efdc8_merged['ST_hr']])\n",
    "efdc8_merged.set_index('date', inplace=True)\n",
    "efdc8_merged.index = efdc8_merged.index.tz_localize(\n",
    "    pytz.timezone('EST')).tz_convert(pytz.utc) #Convert to UTC\n",
    "del efdc8_df\n",
    "\n",
    "\n",
    "efdc9_merged = pd.merge(efdc9_df, efdc_grid_df,  how='left', left_on=['I_MOD','J_MOD'],\n",
    "                    right_on = ['I','J']).drop(columns=[\n",
    "                    ' DUMPID', 'END_hr', 'I_MOD', 'J_MOD', 'I_MOD', 'X', 'Y', 'CUE', 'CVE', 'CUN', 'CVN'])\n",
    "efdc9_merged['date'] =  pd.to_datetime([datetime(2012, 9, 1, 00, 00, 00) +\n",
    "                timedelta(hours=h) for h in efdc9_merged['ST_hr']])\n",
    "efdc9_merged.set_index('date', inplace=True)\n",
    "efdc9_merged.index = efdc9_merged.index.tz_localize(\n",
    "    pytz.timezone('EST')).tz_convert(pytz.utc) #Convert to UTC\n",
    "del efdc9_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "outputs": [],
   "source": [
    "efdc12_merged = pd.merge(efdc12_df, efdc_grid_df,  how='left', left_on=['I_MOD','J_MOD'],\n",
    "                    right_on = ['I','J']).drop(columns=[\n",
    "                    ' DUMPID', 'END_hr', 'I_MOD', 'J_MOD', 'I_MOD', 'X', 'Y', 'CUE', 'CVE', 'CUN', 'CVN'])\n",
    "efdc12_merged['date'] =  pd.to_datetime([datetime(2012, 12, 1, 00, 00, 00) +\n",
    "                timedelta(hours=h) for h in efdc12_merged['ST_hr']])\n",
    "efdc12_merged.set_index('date', inplace=True)\n",
    "efdc12_merged.index = efdc12_merged.index.tz_localize(\n",
    "    pytz.timezone('EST')).tz_convert(pytz.utc) #Convert to UTC\n",
    "del efdc12_df"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [],
   "source": [
    "efdc_merged  = pd.concat([efdc6_merged, efdc7_merged, efdc8_merged, efdc9_merged])\n",
    "del efdc6_merged\n",
    "del efdc7_merged\n",
    "del efdc8_merged\n",
    "del efdc9_merged\n",
    "\n",
    "efdc_merged_gdf = gp.GeoDataFrame(\n",
    "    efdc_merged, geometry=efdc_merged['geometry'], crs=\"EPSG:32118\")\n",
    "del efdc_merged"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%  Merge EFDC dataframes and convert to geodataframe\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "outputs": [],
   "source": [
    "efdc_merged = efdc12_merged\n",
    "\n",
    "efdc_merged_gdf = gp.GeoDataFrame(\n",
    "    efdc_merged, geometry=efdc_merged['geometry'], crs=\"EPSG:32118\")\n",
    "del efdc_merged\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "pycharm": {
     "name": "#%% Find nearest grid point to station FFG\n"
    }
   },
   "outputs": [],
   "source": [
    "nearest_geoms   = nearest_points(Point(gdf_wl_FFG.iloc[0,:].geometry), MultiPoint(efdc_grid_gdf.geometry))\n",
    "station_gdf     = efdc_merged_gdf.loc[efdc_merged_gdf['geometry']==nearest_geoms[1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "outputs": [],
   "source": [
    "nearest_geoms   = nearest_points(Point(gdf_wl_NGG1.iloc[0,:].geometry), MultiPoint(efdc_grid_gdf.geometry))\n",
    "stationNGG_gdf  = efdc_merged_gdf.loc[efdc_merged_gdf['geometry']==nearest_geoms[1]]\n",
    "\n",
    "# efdc_grid_gdf.loc[efdc_grid_gdf['geometry']== nearest_geoms[1]].value"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Find nearest grid point to station BGG\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\arey\\AppData\\Local\\Temp/ipykernel_34632/1501513982.py:1: UserWarning: this is an initial implementation of Parquet/Feather file support and associated metadata.  This is tracking version 0.1.0 of the metadata specification at https://github.com/geopandas/geo-arrow-spec\n",
      "\n",
      "This metadata specification does not yet make stability promises.  We do not yet recommend using this in a production setting unless you are able to rewrite your Parquet/Feather files.\n",
      "\n",
      "To further ignore this warning, you can do: \n",
      "import warnings; warnings.filterwarnings('ignore', message='.*initial implementation of Parquet.*')\n",
      "  efdc_merged_gdf.to_parquet('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/efdc_merged_Dec.parquet',\n",
      "C:\\Users\\arey\\AppData\\Local\\Temp/ipykernel_34632/1501513982.py:4: UserWarning: this is an initial implementation of Parquet/Feather file support and associated metadata.  This is tracking version 0.1.0 of the metadata specification at https://github.com/geopandas/geo-arrow-spec\n",
      "\n",
      "This metadata specification does not yet make stability promises.  We do not yet recommend using this in a production setting unless you are able to rewrite your Parquet/Feather files.\n",
      "\n",
      "To further ignore this warning, you can do: \n",
      "import warnings; warnings.filterwarnings('ignore', message='.*initial implementation of Parquet.*')\n",
      "  station_gdf.to_parquet('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/station_gdf_Dec.parquet',\n",
      "C:\\Users\\arey\\AppData\\Local\\Temp/ipykernel_34632/1501513982.py:7: UserWarning: this is an initial implementation of Parquet/Feather file support and associated metadata.  This is tracking version 0.1.0 of the metadata specification at https://github.com/geopandas/geo-arrow-spec\n",
      "\n",
      "This metadata specification does not yet make stability promises.  We do not yet recommend using this in a production setting unless you are able to rewrite your Parquet/Feather files.\n",
      "\n",
      "To further ignore this warning, you can do: \n",
      "import warnings; warnings.filterwarnings('ignore', message='.*initial implementation of Parquet.*')\n",
      "  stationNGG_gdf.to_parquet('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/stationNGG_gdf_Dec.parquet',\n"
     ]
    }
   ],
   "source": [
    "efdc_merged_gdf.to_parquet('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/efdc_merged_Dec.parquet',\n",
    "                           index=True, compression='gzip')\n",
    "\n",
    "station_gdf.to_parquet('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/station_gdf_Dec.parquet',\n",
    "                      index=True, compression='gzip')\n",
    "\n",
    "stationNGG_gdf.to_parquet('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/stationNGG_gdf_Dec.parquet',\n",
    "                      index=True, compression='gzip')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Save EFDC as Parquet for faster importing\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Exception ignored in: <function WKBReader.__del__ at 0x000002A4D5B2A3A0>\n",
      "Traceback (most recent call last):\n",
      "  File \"C:\\Users\\arey\\Anaconda3\\envs\\D3DFM\\lib\\site-packages\\shapely\\geos.py\", line 421, in __del__\n",
      "    self._lgeos.GEOSWKBReader_destroy(self._reader)\n",
      "KeyboardInterrupt: \n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m",
      "\u001B[1;31mKeyboardInterrupt\u001B[0m                         Traceback (most recent call last)",
      "\u001B[1;32m~\\AppData\\Local\\Temp/ipykernel_15068/921050448.py\u001B[0m in \u001B[0;36m<module>\u001B[1;34m\u001B[0m\n\u001B[0;32m      2\u001B[0m \u001B[0mstation_gdf\u001B[0m    \u001B[1;33m=\u001B[0m \u001B[0mgp\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mread_parquet\u001B[0m\u001B[1;33m(\u001B[0m\u001B[1;34m'C:/Users/arey/files/Projects/Newtown/RevisedEFDC/station_gdf.parquet'\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m      3\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m----> 4\u001B[1;33m \u001B[0mefdc_merged_gdf\u001B[0m \u001B[1;33m=\u001B[0m \u001B[0mgp\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mread_parquet\u001B[0m\u001B[1;33m(\u001B[0m\u001B[1;34m'C:/Users/arey/files/Projects/Newtown/RevisedEFDC/efdc_merged.parquet'\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m      5\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;32m~\\Anaconda3\\envs\\D3DFM\\lib\\site-packages\\geopandas\\io\\arrow.py\u001B[0m in \u001B[0;36m_read_parquet\u001B[1;34m(path, columns, **kwargs)\u001B[0m\n\u001B[0;32m    393\u001B[0m     \u001B[0mtable\u001B[0m \u001B[1;33m=\u001B[0m \u001B[0mparquet\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mread_table\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mpath\u001B[0m\u001B[1;33m,\u001B[0m \u001B[0mcolumns\u001B[0m\u001B[1;33m=\u001B[0m\u001B[0mcolumns\u001B[0m\u001B[1;33m,\u001B[0m \u001B[1;33m**\u001B[0m\u001B[0mkwargs\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    394\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m--> 395\u001B[1;33m     \u001B[1;32mreturn\u001B[0m \u001B[0m_arrow_to_geopandas\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mtable\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m    396\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    397\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;32m~\\Anaconda3\\envs\\D3DFM\\lib\\site-packages\\geopandas\\io\\arrow.py\u001B[0m in \u001B[0;36m_arrow_to_geopandas\u001B[1;34m(table)\u001B[0m\n\u001B[0;32m    335\u001B[0m     \u001B[1;31m# Convert the WKB columns that are present back to geometry.\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    336\u001B[0m     \u001B[1;32mfor\u001B[0m \u001B[0mcol\u001B[0m \u001B[1;32min\u001B[0m \u001B[0mgeometry_columns\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m--> 337\u001B[1;33m         \u001B[0mdf\u001B[0m\u001B[1;33m[\u001B[0m\u001B[0mcol\u001B[0m\u001B[1;33m]\u001B[0m \u001B[1;33m=\u001B[0m \u001B[0mfrom_wkb\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mdf\u001B[0m\u001B[1;33m[\u001B[0m\u001B[0mcol\u001B[0m\u001B[1;33m]\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mvalues\u001B[0m\u001B[1;33m,\u001B[0m \u001B[0mcrs\u001B[0m\u001B[1;33m=\u001B[0m\u001B[0mmetadata\u001B[0m\u001B[1;33m[\u001B[0m\u001B[1;34m\"columns\"\u001B[0m\u001B[1;33m]\u001B[0m\u001B[1;33m[\u001B[0m\u001B[0mcol\u001B[0m\u001B[1;33m]\u001B[0m\u001B[1;33m[\u001B[0m\u001B[1;34m\"crs\"\u001B[0m\u001B[1;33m]\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m    338\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    339\u001B[0m     \u001B[1;32mreturn\u001B[0m \u001B[0mGeoDataFrame\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mdf\u001B[0m\u001B[1;33m,\u001B[0m \u001B[0mgeometry\u001B[0m\u001B[1;33m=\u001B[0m\u001B[0mgeometry\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;32m~\\Anaconda3\\envs\\D3DFM\\lib\\site-packages\\geopandas\\array.py\u001B[0m in \u001B[0;36mfrom_wkb\u001B[1;34m(data, crs)\u001B[0m\n\u001B[0;32m    192\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    193\u001B[0m     \"\"\"\n\u001B[1;32m--> 194\u001B[1;33m     \u001B[1;32mreturn\u001B[0m \u001B[0mGeometryArray\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mvectorized\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mfrom_wkb\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mdata\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m,\u001B[0m \u001B[0mcrs\u001B[0m\u001B[1;33m=\u001B[0m\u001B[0mcrs\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m    195\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    196\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;32m~\\Anaconda3\\envs\\D3DFM\\lib\\site-packages\\geopandas\\_vectorized.py\u001B[0m in \u001B[0;36mfrom_wkb\u001B[1;34m(data)\u001B[0m\n\u001B[0;32m    167\u001B[0m     \u001B[1;32mfor\u001B[0m \u001B[0mgeom\u001B[0m \u001B[1;32min\u001B[0m \u001B[0mdata\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    168\u001B[0m         \u001B[1;32mif\u001B[0m \u001B[0mgeom\u001B[0m \u001B[1;32mis\u001B[0m \u001B[1;32mnot\u001B[0m \u001B[1;32mNone\u001B[0m \u001B[1;32mand\u001B[0m \u001B[0mlen\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mgeom\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m--> 169\u001B[1;33m             \u001B[0mgeom\u001B[0m \u001B[1;33m=\u001B[0m \u001B[0mshapely\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mwkb\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mloads\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mgeom\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m    170\u001B[0m         \u001B[1;32melse\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    171\u001B[0m             \u001B[0mgeom\u001B[0m \u001B[1;33m=\u001B[0m \u001B[1;32mNone\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;32m~\\Anaconda3\\envs\\D3DFM\\lib\\site-packages\\shapely\\wkb.py\u001B[0m in \u001B[0;36mloads\u001B[1;34m(data, hex)\u001B[0m\n\u001B[0;32m     15\u001B[0m         \u001B[1;32mreturn\u001B[0m \u001B[0mreader\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mread_hex\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mdata\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m     16\u001B[0m     \u001B[1;32melse\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m---> 17\u001B[1;33m         \u001B[1;32mreturn\u001B[0m \u001B[0mreader\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mread\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mdata\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m     18\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m     19\u001B[0m \u001B[1;32mdef\u001B[0m \u001B[0mload\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mfp\u001B[0m\u001B[1;33m,\u001B[0m \u001B[0mhex\u001B[0m\u001B[1;33m=\u001B[0m\u001B[1;32mFalse\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;32m~\\Anaconda3\\envs\\D3DFM\\lib\\site-packages\\shapely\\geos.py\u001B[0m in \u001B[0;36mread\u001B[1;34m(self, data)\u001B[0m\n\u001B[0;32m    433\u001B[0m         \u001B[1;31m# avoid circular import dependency\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    434\u001B[0m         \u001B[1;32mfrom\u001B[0m \u001B[0mshapely\u001B[0m \u001B[1;32mimport\u001B[0m \u001B[0mgeometry\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m--> 435\u001B[1;33m         \u001B[1;32mreturn\u001B[0m \u001B[0mgeometry\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mbase\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mgeom_factory\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mgeom\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m    436\u001B[0m \u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m    437\u001B[0m     \u001B[1;32mdef\u001B[0m \u001B[0mread_hex\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mself\u001B[0m\u001B[1;33m,\u001B[0m \u001B[0mdata\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;32m~\\Anaconda3\\envs\\D3DFM\\lib\\site-packages\\shapely\\geometry\\base.py\u001B[0m in \u001B[0;36mgeom_factory\u001B[1;34m(g, parent)\u001B[0m\n\u001B[0;32m     88\u001B[0m     \u001B[0mob\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0m__class__\u001B[0m \u001B[1;33m=\u001B[0m \u001B[0mgetattr\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mmod\u001B[0m\u001B[1;33m,\u001B[0m \u001B[0mgeom_type\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m     89\u001B[0m     \u001B[0mob\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0m_geom\u001B[0m \u001B[1;33m=\u001B[0m \u001B[0mg\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[1;32m---> 90\u001B[1;33m     \u001B[0mob\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0m__p__\u001B[0m \u001B[1;33m=\u001B[0m \u001B[0mparent\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0m\u001B[0;32m     91\u001B[0m     \u001B[1;32mif\u001B[0m \u001B[0mlgeos\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0mmethods\u001B[0m\u001B[1;33m[\u001B[0m\u001B[1;34m'has_z'\u001B[0m\u001B[1;33m]\u001B[0m\u001B[1;33m(\u001B[0m\u001B[0mg\u001B[0m\u001B[1;33m)\u001B[0m\u001B[1;33m:\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n\u001B[0;32m     92\u001B[0m         \u001B[0mob\u001B[0m\u001B[1;33m.\u001B[0m\u001B[0m_ndim\u001B[0m \u001B[1;33m=\u001B[0m \u001B[1;36m3\u001B[0m\u001B[1;33m\u001B[0m\u001B[1;33m\u001B[0m\u001B[0m\n",
      "\u001B[1;31mKeyboardInterrupt\u001B[0m: "
     ]
    }
   ],
   "source": [
    "stationNGG_gdf  = gp.read_parquet('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/stationNGG_gdf.parquet')\n",
    "station_gdf     = gp.read_parquet('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/station_gdf.parquet')\n",
    "\n",
    "efdc_merged_gdf = gp.read_parquet('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/efdc_merged.parquet')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Read in EFDC as Paraquet\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [],
   "source": [
    "runLog = pd.read_excel('C:/Users/arey/files\\Projects/Newtown/Model Log NTC.xlsx', 'ModelLog')\n",
    "\n",
    "dataPath = \"FlowFM_map.nc\""
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Read Model Log\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reading time dimension: reconstruct array\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "analyzing netdata from domain 0000 of 0000 (counting max number of facenodes)\n",
      "processing netdata from domain 0000 of 0000\n",
      "reading time dimension: only requested indices\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "reading time dimension: only requested indices\n",
      "retrieving original timezone succeeded, no conversion to UTC/GMT applied\n",
      "varname: mesh2d_s1  (337, 11096)  ('time', 'mesh2d_nFaces'), coordinates=(mesh2d_face_x mesh2d_face_y)\n",
      "varname: mesh2d_face_x  (11096,)  ('mesh2d_nFaces',), coordinates=(None)\n",
      "varname: mesh2d_face_y  (11096,)  ('mesh2d_nFaces',), coordinates=(None)\n"
     ]
    }
   ],
   "source": [
    "modelPlot  = [31]\n",
    "\n",
    "ugrid_all           = [None] * (max(modelPlot)+1)\n",
    "data_frommap_wl     = [None] * (max(modelPlot)+1)\n",
    "facex               = [None] * (max(modelPlot)+1)\n",
    "facey               = [None] * (max(modelPlot)+1)\n",
    "\n",
    "nearest_IDX     = list(range(100))\n",
    "nearest_IDX_NGG = list(range(100))\n",
    "dsloc           = list(range(100))\n",
    "dsloc_NGG       = list(range(100))\n",
    "\n",
    "# Find nearest point\n",
    "for i in modelPlot:\n",
    "    dataPath            = pl.Path(runLog['Run Location'][i],\n",
    "                                  'FlowFM', 'dflowfm', 'output', 'FlowFM_map.nc')\n",
    "\n",
    "    file_nc_map         = dataPath.as_posix()\n",
    "\n",
    "    tSteps              = get_timesfromnc(file_nc=file_nc_map, varname='time')\n",
    "    tStep               = len(tSteps)-1\n",
    "\n",
    "    ugrid_all[i]        = get_netdata(file_nc=file_nc_map)#,multipart=False)\n",
    "    data_frommap_wl[i]  = get_ncmodeldata(file_nc=file_nc_map, varname='mesh2d_s1', timestep=tStep)\n",
    "    facex[i]            = get_ncmodeldata(file_nc=file_nc_map, varname='mesh2d_face_x')\n",
    "    facey[i]            = get_ncmodeldata(file_nc=file_nc_map, varname='mesh2d_face_y')\n",
    "\n",
    "    # Find matching Delft3D index\n",
    "    s = gp.GeoSeries(map(Point, zip(facex[i], facey[i])))\n",
    "    nearest_geoms       = nearest_points(Point(gdf_wl_FFG.iloc[0,:].geometry), MultiPoint(s))\n",
    "    nearest_IDX[i]      = np.where(s==nearest_geoms[1])\n",
    "\n",
    "    nearest_geoms       = nearest_points(Point(gdf_wl_NGG1.iloc[0,:].geometry), MultiPoint(s))\n",
    "    nearest_IDX[i]      = np.where(s==nearest_geoms[1])\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Import Delft3D results\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [],
   "source": [
    "modelPlot  = [31]\n",
    "\n",
    "ugrid_all           = [None] * (max(modelPlot)+1)\n",
    "data_frommap_wl     = [None] * (max(modelPlot)+1)\n",
    "facex               = [None] * (max(modelPlot)+1)\n",
    "facey               = [None] * (max(modelPlot)+1)\n",
    "\n",
    "nearest_IDX     = list(range(100))\n",
    "nearest_IDX_NGG = list(range(100))\n",
    "dsloc           = list(range(100))\n",
    "dsloc_NGG       = list(range(100))\n",
    "\n",
    "# Note, this uses the standard netcdf lib + xarray\n",
    "for i in modelPlot:\n",
    "    dataPath            = pl.Path(runLog['Run Location'][i],\n",
    "                                  'FlowFM', 'dflowfm', 'output', 'FlowFM_map.nc')\n",
    "\n",
    "    file_nc_map         = dataPath.as_posix()\n",
    "\n",
    "    ds = xr.open_dataset(file_nc_map)\n",
    "\n",
    "    # Find nearest point\n",
    "    ds_x = ds['mesh2d_face_x']\n",
    "    ds_y = ds['mesh2d_face_y']\n",
    "\n",
    "    s = gp.GeoSeries(map(Point, zip(ds_x, ds_y)))\n",
    "    nearest_geoms       = nearest_points(Point(gdf_wl_FFG.iloc[0,:].geometry), MultiPoint(s))\n",
    "    nearest_IDX[i]      = np.where(s==nearest_geoms[1])\n",
    "\n",
    "    nearest_geoms_NGG   = nearest_points(Point(gdf_wl_NGG1.iloc[0,:].geometry), MultiPoint(s))\n",
    "    nearest_IDX_NGG[i]  = np.where(s==nearest_geoms_NGG[1])\n",
    "\n",
    "    dsloc[i]       = ds['mesh2d_s1'].sel(mesh2d_nFaces=nearest_IDX[i][0])\n",
    "    dsloc_NGG[i]   = ds['mesh2d_s1'].sel(mesh2d_nFaces=nearest_IDX_NGG[i][0])"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Read in time series\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "outputs": [],
   "source": [
    "modelPlot  = [31]\n",
    "\n",
    "df_wl_D3D       = list(range(100))\n",
    "df_wl_D3D_NGG   = list(range(100))\n",
    "\n",
    "d3dTzones       = list(range(100))\n",
    "d3dTzones[25]   = 'EST'\n",
    "d3dTzones[26]   = 'EST'\n",
    "d3dTzones[31]   = 'EST'\n",
    "\n",
    "for i in modelPlot:\n",
    "    df_wl_D3D[i] = dsloc[i].to_dataframe()\n",
    "    df_wl_D3D[i].index = df_wl_D3D[i].index.droplevel(-1)\n",
    "\n",
    "    df_wl_D3D[i].index = df_wl_D3D[i].index.tz_localize(\n",
    "        pytz.timezone(d3dTzones[i])).tz_convert(pytz.utc) #Convert to UTC\n",
    "\n",
    "    df_wl_D3D_NGG[i] = dsloc_NGG[i].to_dataframe()\n",
    "    df_wl_D3D_NGG[i].index = df_wl_D3D_NGG[i].index.droplevel(-1)\n",
    "\n",
    "    df_wl_D3D_NGG[i].index = df_wl_D3D_NGG[i].index.tz_localize(\n",
    "        pytz.timezone(d3dTzones[i])).tz_convert(pytz.utc) #Convert to UTC"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Convert D3D Timezone\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "pycharm": {
     "name": "#%% Plot water levels at nearest point\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 648x360 with 2 Axes>",
      "image/png": "[encoded data, 97688 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 5))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "for week in range(0, 2):\n",
    "    weekOffset = 0\n",
    "    axes[week+weekOffset].plot(station_gdf.index, station_gdf['WSE_m'], label='EFDC Model')\n",
    "\n",
    "    axes[week+weekOffset].plot(gdf_wl_FFG.index+timedelta(hours=0),\n",
    "              gdf_wl_FFG['Water Surface Elevation (ft)']*0.3048,\n",
    "              label='Field Facility Gauge', color='k')\n",
    "\n",
    "    for i in modelPlot:\n",
    "        axes[week+weekOffset].plot(df_wl_D3D[i].index, df_wl_D3D[i].mesh2d_s1, label=runLog['Run'][i])\n",
    "\n",
    "    axes[week+weekOffset].set_xlim(datetime(2012, 12,  15, 0, 0, 0) + timedelta(days=week*7),\n",
    "                        datetime(2012, 12, 15+7, 0, 0, 0) + timedelta(days=week*7))\n",
    "    axes[week+weekOffset].set_ylim(-1, 2.5)\n",
    "\n",
    "\n",
    "    axes[week+weekOffset].set_ylabel('Water Surface Elevation [m]')\n",
    "\n",
    "axes[0].legend(loc='upper left')\n",
    "axes[0].title.set_text('Field Facility Gauge')\n",
    "\n",
    "plt.show()\n",
    "fig.savefig('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/Figures/DecTS_FFG.png',\n",
    "             bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 648x360 with 2 Axes>",
      "image/png": "[encoded data, 100556 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 5))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "for week in range(0, 2):\n",
    "    weekOffset = 0\n",
    "    axes[week+weekOffset].plot(stationNGG_gdf.index, stationNGG_gdf['WSE_m'], label='EFDC Model')\n",
    "\n",
    "    axes[week+weekOffset].plot(gdf_wl_NGG1.index+timedelta(hours=0),\n",
    "              gdf_wl_NGG1['Water Surface Elevation (ft)']*0.3048,\n",
    "              label='National Grid Gauge', color='k')\n",
    "\n",
    "    for i in modelPlot:\n",
    "        axes[week+weekOffset].plot(df_wl_D3D_NGG[i].index, df_wl_D3D_NGG[i].mesh2d_s1,\n",
    "                                   label=runLog['Run'][i])\n",
    "\n",
    "    axes[week+weekOffset].set_xlim(datetime(2012, 12,  15, 0, 0, 0) + timedelta(days=week*7),\n",
    "                        datetime(2012, 12, 15+7, 0, 0, 0) + timedelta(days=week*7))\n",
    "    axes[week+weekOffset].set_ylim(-1, 2.5)\n",
    "\n",
    "    axes[week+weekOffset].set_ylabel('Water Surface Elevation [m]')\n",
    "\n",
    "axes[0].legend(loc='upper left')\n",
    "axes[0].title.set_text('National Grid Gauge')\n",
    "\n",
    "plt.show()\n",
    "fig.savefig('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/Figures/DecTS_NGG.png',\n",
    "             bbox_inches='tight')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Plot water levels at National Grid\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 720x720 with 1 Axes>",
      "image/png": "[encoded data, 79380 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 10))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "interpTimes = pd.date_range(\"2012-06-02T00:00:00\", \"2012-10-01T00:00:00\", freq=\"60min\")\n",
    "interpTimes = pd.date_range(\"2012-07-07T00:00:00\", \"2012-07-27T00:00:00\", freq=\"60min\")\n",
    "\n",
    "for i in modelPlot:\n",
    "    axes.scatter(np.interp(interpTimes, gdf_wl_FFG.index,\n",
    "             gdf_wl_FFG['Water Surface Elevation (ft)']*0.3048),\n",
    "             np.interp(interpTimes, df_wl_D3D[i].index,\n",
    "             df_wl_D3D[i].mesh2d_s1),\n",
    "             label=runLog['Run'][i])\n",
    "\n",
    "axes.scatter(np.interp(interpTimes, gdf_wl_FFG.index,\n",
    "             gdf_wl_FFG['Water Surface Elevation (ft)']*0.3048),\n",
    "             np.interp(interpTimes, station_gdf.index,\n",
    "             station_gdf['WSE_m']),\n",
    "             label='EFDC')\n",
    "\n",
    "linepts = np.array([-1.0, 1.5])\n",
    "\n",
    "plt.plot(linepts, linepts, color='k')\n",
    "\n",
    "axes.set_xlim(-1.0, 1.5)\n",
    "axes.set_ylim(-1.0, 1.5)\n",
    "axes.set_aspect('equal')\n",
    "axes.set_title('Field Facility Gauge')\n",
    "axes.legend(loc='upper left')\n",
    "axes.set_ylabel('Model Water Surface Elevation [m]')\n",
    "axes.set_xlabel('Observed Water Surface Elevation [m]')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "fig.savefig('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/Figures/FFG_ScatterUpdated.pdf',\n",
    "             bbox_inches='tight')\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Plot scatters\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 720x720 with 1 Axes>",
      "image/png": "[encoded data, 80164 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 10))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "interpTimes = pd.date_range(\"2012-06-02T00:00:00\", \"2012-10-01T00:00:00\", freq=\"60min\")\n",
    "interpTimes = pd.date_range(\"2012-07-07T00:00:00\", \"2012-07-27T00:00:00\", freq=\"60min\")\n",
    "\n",
    "for i in modelPlot:\n",
    "    axes.scatter(np.interp(interpTimes, gdf_wl_NGG1.index,\n",
    "             gdf_wl_NGG1['Water Surface Elevation (ft)']*0.3048),\n",
    "             np.interp(interpTimes, df_wl_D3D_NGG[i].index,\n",
    "             df_wl_D3D_NGG[i].mesh2d_s1),\n",
    "             label=runLog['Run'][i])\n",
    "\n",
    "axes.scatter(np.interp(interpTimes, gdf_wl_NGG1.index,\n",
    "             gdf_wl_NGG1['Water Surface Elevation (ft)']*0.3048),\n",
    "             np.interp(interpTimes, stationNGG_gdf.index,\n",
    "             stationNGG_gdf['WSE_m']),\n",
    "             label='EFDC')\n",
    "\n",
    "linepts = np.array([-1.0, 1.5])\n",
    "\n",
    "plt.plot(linepts, linepts, color='k')\n",
    "\n",
    "axes.set_xlim(-1.0, 1.5)\n",
    "axes.set_ylim(-1.0, 1.5)\n",
    "axes.set_aspect('equal')\n",
    "axes.set_title('National Grid Gauge')\n",
    "axes.legend(loc='upper left')\n",
    "axes.set_ylabel('Model Water Surface Elevation [m]')\n",
    "axes.set_xlabel('Observed Water Surface Elevation [m]')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "fig.savefig('C:/Users/arey/files/Projects/Newtown/RevisedEFDC/Figures/NGG_ScatterUpdated.pdf',\n",
    "             bbox_inches='tight')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Plot scatters NGG\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PyCharm (NTC_DFM)",
   "language": "python",
   "name": "pycharm-b2185061"
  },
  "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
}