{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% Setup\n"
    }
   },
   "outputs": [],
   "source": [
    "#import jupyter\n",
    "import pandas as pd\n",
    "import geopandas as gp\n",
    "import gpxpy\n",
    "import gpxpy.gpx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import contextily as ctx\n",
    "import datetime\n",
    "\n",
    "import scipy as sp\n",
    "\n",
    "from shapely import geometry, ops"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "pycharm": {
     "name": "#%% Load in Mill Site Data\n"
    }
   },
   "outputs": [],
   "source": [
    "# From CSV\n",
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/03_Hg fluxes and mass balances/Mill site data/All Samples FINAL.xlsx',\n",
    "                 sheet_name='Sheet1', skiprows=4)\n",
    "\n",
    "# Load Mill Locations\n",
    "gpx_file = open('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/03_Hg fluxes and mass balances/Mill site data/November fishing trip.gpx', 'r', encoding='utf-8-sig')\n",
    "gpx = gpxpy.parse(gpx_file)\n",
    "\n",
    "millLocs = []\n",
    "for waypoint in gpx.waypoints:\n",
    "    millLocs.append([waypoint.name, waypoint.latitude, waypoint.longitude, waypoint.time])\n",
    "\n",
    "# Convert to dataframe\n",
    "mill_df  = pd.DataFrame(millLocs, columns=['Name', 'Latitude', 'Longitude', 'Time'])\n",
    "# Change name strings to numbers\n",
    "mill_df['Name'] = mill_df['Name'].astype(float)\n",
    "df_in = df_in.set_index('Waypoint').join(mill_df.set_index('Name'))\n",
    "\n",
    "mill_gdf = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.Longitude, df_in.Latitude, crs=\"EPSG:4326\"))\n",
    "mill_gdf = mill_gdf.assign(media='SED')\n",
    "\n",
    "#Merge\n",
    "del df_in, gpx, mill_df, millLocs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "pycharm": {
     "name": "#%% Cores\n"
    }
   },
   "outputs": [],
   "source": [
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/06_Sediment Hg and characteristics/02_Cores/Core_Data_Summary_AJMR.xlsx',\n",
    "                 sheet_name='SUMMARY')\n",
    "core_gdf = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.Longitude, df_in.Latitude, crs=\"EPSG:4326\"))\n",
    "core_gdf = core_gdf.assign(media='SED')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "pycharm": {
     "name": "#%% Sediment\n"
    }
   },
   "outputs": [],
   "source": [
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/06_Sediment Hg and characteristics/04_Sediment/Wabigoon River_Spring 2017 MOECC Surface Sediment Sampling_Mercury TOC P_AJMR.xlsx',\n",
    "                       sheet_name='Data_Transpose')\n",
    "\n",
    "sed_gdf = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.Longitude, df_in.Latitude, crs=\"EPSG:4326\"))\n",
    "sed_gdf = sed_gdf.assign(media='SED')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "pycharm": {
     "name": "#%% Water Col\n"
    }
   },
   "outputs": [],
   "source": [
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/07_Water column Hg/Wabigoon River System Data  22Aug TS SE5ck Se 20 Nov 8 rch_AJMR.xlsx',\n",
    "                       sheet_name='Data Values & transposed_AJMR',skiprows=29)\n",
    "\n",
    "# Split date columns into separate rows\n",
    "\n",
    "# df_process = pd.DataFrame()\n",
    "i = 7\n",
    "s = 1\n",
    "df_process = df_in.iloc[0:13, [0, 1, 2, 3, 4, 5, 6, i]]\n",
    "df_process.loc[0:13, 'date'] = df_in.columns[i]\n",
    "# 'date' is col 8\n",
    "for i in range(8,20):\n",
    "    df_process = pd.concat([df_process, df_in.iloc[14*(s):14*(s)+13, [0, 1, 2, 3, 4, 5, 6, i]]])\n",
    "    df_process.iloc[-13:, 8] = df_in.columns[i]\n",
    "\n",
    "for s in range(2,7):\n",
    "    for i in range(7,20):\n",
    "        df_process = pd.concat([df_process, df_in.iloc[14*(s):14*(s)+13, [0, 1, 2, 3, 4, 5, 6, i]]])\n",
    "        df_process.iloc[-13:, 8] = df_in.columns[i]\n",
    "\n",
    "# Sum date cols to get one data col\n",
    "df_process['data'] = df_process.iloc[:, np.r_[7, 9:21]].sum(axis=1)\n",
    "\n",
    "# Drop other cols\n",
    "df_process.drop(columns=df_process.iloc[:, np.r_[7, 9:21]].columns.tolist(), inplace=True)\n",
    "\n",
    "waterCol_gdf = gp.GeoDataFrame(df_process, geometry=gp.points_from_xy(df_process.Lon, df_process.Lat, crs=\"EPSG:4326\"))\n",
    "\n",
    "waterCol_gdf = waterCol_gdf.assign(media='WATERCOL')\n",
    "del df_in, df_process, s, i\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "pycharm": {
     "name": "#%% Wetland\n"
    }
   },
   "outputs": [],
   "source": [
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/07_Water column Hg/Wabigoon River System Data  22Aug TS SE5ck Se 20 Nov 8 rch_AJMR.xlsx',\n",
    "                       sheet_name='Wetland sampling',skiprows=3)\n",
    "wetland_gdf = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.Long, df_in.Lat, crs=\"EPSG:4326\"))\n",
    "\n",
    "wetland_gdf = wetland_gdf.assign(media='WATERCOL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "pycharm": {
     "name": "#%% Parks River Data\n"
    }
   },
   "outputs": [],
   "source": [
    "parksPTS = gp.read_file('C:/Users/arey/files/Projects/Grassy Narrows/Hg Review/ParksPoints.shp')\n",
    "# Change to int to match excel data\n",
    "parksPTS = parksPTS.astype({'Name': 'int64'})\n",
    "parksPTS = parksPTS.to_crs(\"EPSG:4326\")\n",
    "\n",
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/07_Water column Hg/Parks () Wabigoon River field data J_PARKS_DEC_16_2010.xlsx',\n",
    "                       sheet_name='J_PARKS_DEC_16_2010')\n",
    "#Join locations\n",
    "df_in = df_in.set_index('STATION_NO').join(parksPTS.set_index('Name'))# wetland_gdf = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.Long, df_in.Lat, crs=\"EPSG:4326\"))\n",
    "\n",
    "parks_gdf = gp.GeoDataFrame(df_in, geometry=df_in.geometry)\n",
    "\n",
    "parks_gdf = parks_gdf.assign(media='WATERCOL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "pycharm": {
     "name": "#%% Temp DO data- River\n"
    }
   },
   "outputs": [],
   "source": [
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/07_Water column Hg/Temp  DO lake and river and MeHg profiles Grassy Project 2017 Updated Sept 7 TS Dec ck May CK_AJMR.xlsx',\n",
    "                        sheet_name='River Data_Transpose')\n",
    "\n",
    "riverTmpDO_gdf = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.Lon, df_in.Lat, crs=\"EPSG:4326\"))\n",
    "riverTmpDO_gdf = riverTmpDO_gdf.assign(media='WATERCOL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "pycharm": {
     "name": "#%% Temp DO data- Ball Lake\n"
    }
   },
   "outputs": [],
   "source": [
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/07_Water column Hg/Temp  DO lake and river and MeHg profiles Grassy Project 2017 Updated Sept 7 TS Dec ck May CK_AJMR.xlsx',\n",
    "                        sheet_name='Ball Lake')\n",
    "\n",
    "ballTmpDO_gdf = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.Lon, df_in.Lat, crs=\"EPSG:4326\"))\n",
    "ballTmpDO_gdf = ballTmpDO_gdf.assign(media='WATERCOL')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "pycharm": {
     "name": "#%% Temp DO data- Clay Lake\n"
    }
   },
   "outputs": [],
   "source": [
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/02_Field data - English Wabigoon/07_Water column Hg/Temp  DO lake and river and MeHg profiles Grassy Project 2017 Updated Sept 7 TS Dec ck May CK_AJMR.xlsx',\n",
    "                        sheet_name='Clay Lake')\n",
    "\n",
    "clayTmpDO_gdf = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.Lon, df_in.Lat, crs=\"EPSG:4326\"))\n",
    "clayTmpDO_gdf = clayTmpDO_gdf.assign(media='WATERCOL')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "pycharm": {
     "name": "#%% Consolidated Sheet\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\arey\\Anaconda3\\envs\\py39\\lib\\site-packages\\openpyxl\\worksheet\\header_footer.py:48: UserWarning: Cannot parse header or footer so it will be ignored\n",
      "  warn(\"\"\"Cannot parse header or footer so it will be ignored\"\"\")\n"
     ]
    }
   ],
   "source": [
    "df_in = pd.read_excel('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/25_Water Quality/EW_Sed_MeHg_More_AJMR.xlsx',\n",
    "                        sheet_name='EW_Sed_MeHg_More')\n",
    "\n",
    "# Remove points with no geodata\n",
    "df_in.dropna(subset='X_coord', inplace=True)\n",
    "\n",
    "consolidate_gdf = gp.GeoDataFrame(df_in, geometry=gp.points_from_xy(df_in.X_coord, df_in.Y_coord, crs=\"EPSG:32615\"))\n",
    "consolidate_gdf = consolidate_gdf.to_crs(\"EPSG:4326\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "pycharm": {
     "name": "#%% Merge\n"
    }
   },
   "outputs": [],
   "source": [
    "temp_gdf = consolidate_gdf.rename({'field_sample_date': 'date', 'mean_depth': 'depth', 'final_result':'data'}, axis=1)\n",
    "temp_gdf['date'] = pd.to_datetime(temp_gdf['date'], format='%Y-%m-%d %H:%M:%S')\n",
    "temp_gdf = temp_gdf.assign(Source='Consolidated')\n",
    "EWR_gdf = temp_gdf\n",
    "\n",
    "temp_gdf = mill_gdf.rename({'Sample Name': 'loc_name', 'Depth': 'depth', 'Concentration Dry Weight  (ng/g)':'data', 'Time':'date'}, axis=1)\n",
    "temp_gdf = temp_gdf.assign(param_name='Mercury')\n",
    "temp_gdf = temp_gdf.assign(result_uom='NG/G')\n",
    "temp_gdf = temp_gdf.assign(Source='Mill')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "temp_gdf = core_gdf.rename({'Site': 'loc_name', 'Mean Depth (m)': 'depth', 'Average ng/g DW':'data', 'Date2':'date'}, axis=1)\n",
    "temp_gdf['date']  = pd.to_datetime(temp_gdf['date'])\n",
    "temp_gdf['depth'] = temp_gdf['depth']*100 # convert to cm\n",
    "temp_gdf = temp_gdf.assign(param_name='Mercury')\n",
    "temp_gdf = temp_gdf.assign(result_uom='NG/G')\n",
    "temp_gdf = temp_gdf.assign(Source='Cores')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "temp_gdf = sed_gdf.rename({'ID': 'loc_name', 'Mercury (ng/g)':'data', 'Datetime':'date'}, axis=1)\n",
    "temp_gdf['date']  = pd.to_datetime(temp_gdf['date'])\n",
    "temp_gdf = temp_gdf.assign(depth=0)\n",
    "temp_gdf = temp_gdf.assign(param_name='Mercury')\n",
    "temp_gdf = temp_gdf.assign(result_uom='NG/G')\n",
    "temp_gdf = temp_gdf.assign(Source='Sediment')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "temp_gdf = waterCol_gdf.rename({'Site Name': 'loc_name', 'Datetime':'date', 'Unit':'result_uom', 'Variable':'param_name'}, axis=1)\n",
    "temp_gdf['date']  = pd.to_datetime(temp_gdf['date'])\n",
    "temp_gdf = temp_gdf.assign(depth=0)\n",
    "temp_gdf = temp_gdf.assign(Source='WaterCol')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "temp_gdf = wetland_gdf.rename({'Unnamed: 0': 'loc_name', 'Date':'date', 'Bulk THG, ng/L':'data'}, axis=1)\n",
    "temp_gdf['date']  = pd.to_datetime(temp_gdf['date'])\n",
    "temp_gdf = temp_gdf.assign(depth=0)\n",
    "temp_gdf = temp_gdf.assign(param_name='Mercury')\n",
    "temp_gdf = temp_gdf.assign(result_uom='NG/L')\n",
    "temp_gdf = temp_gdf.assign(Source='Wetland')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "temp_gdf = wetland_gdf.rename({'Unnamed: 0': 'loc_name', 'Date':'date', 'Bulk MeHg, ng/L':'data'}, axis=1)\n",
    "temp_gdf['date']  = pd.to_datetime(temp_gdf['date'])\n",
    "temp_gdf = temp_gdf.assign(depth=0)\n",
    "temp_gdf = temp_gdf.assign(param_name='MeHg')\n",
    "temp_gdf = temp_gdf.assign(result_uom='NG/L')\n",
    "temp_gdf = temp_gdf.assign(Source='Wetland')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "temp_gdf = wetland_gdf.rename({'Unnamed: 0': 'loc_name', 'Date':'date', '%MeHg':'data'}, axis=1)\n",
    "temp_gdf['date']  = pd.to_datetime(temp_gdf['date'])\n",
    "temp_gdf = temp_gdf.assign(depth=0)\n",
    "temp_gdf = temp_gdf.assign(param_name='%MeHg')\n",
    "temp_gdf = temp_gdf.assign(result_uom='%')\n",
    "temp_gdf = temp_gdf.assign(Source='Wetland')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "temp_gdf = wetland_gdf.rename({'Unnamed: 0': 'loc_name', 'Date':'date', 'ng/gdw':'data'}, axis=1)\n",
    "temp_gdf['date']  = pd.to_datetime(temp_gdf['date'])\n",
    "temp_gdf = temp_gdf.assign(depth=0)\n",
    "temp_gdf = temp_gdf.assign(param_name='MeHg')\n",
    "temp_gdf = temp_gdf.assign(result_uom='ng/gdw')\n",
    "temp_gdf = temp_gdf.assign(Source='Wetland')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "temp_gdf = parks_gdf.rename({'RESULT':'data', 'DATE':'date', 'UNITS':'result_uom', 'PARM':'param_name'}, axis=1)\n",
    "temp_gdf['date']  = pd.to_datetime(temp_gdf['date'])\n",
    "temp_gdf = temp_gdf.assign(depth=0)\n",
    "temp_gdf = temp_gdf.assign(loc_name=parks_gdf.index.values)\n",
    "temp_gdf = temp_gdf.assign(Source='Parks')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "temp_gdf = clayTmpDO_gdf.rename({'Unnamed: 1':'loc_name', 'Unnamed: 0':'date', 'm':'depth', 'MeHg ng/K':'data'}, axis=1)\n",
    "temp_gdf['date']  = pd.to_datetime(temp_gdf['date'])\n",
    "temp_gdf = temp_gdf.assign(param_name='MeHg')\n",
    "temp_gdf = temp_gdf.assign(result_uom='ng/L')\n",
    "temp_gdf = temp_gdf.assign(Source='ClayTempDO')\n",
    "EWR_gdf = pd.concat([EWR_gdf, temp_gdf])\n",
    "\n",
    "# Remove empty\n",
    "EWR_gdf = EWR_gdf.loc[EWR_gdf['geometry'].is_valid, :]\n",
    "EWR_gdf['date'] = pd.to_datetime(EWR_gdf['date'], utc=True)\n",
    "EWR_gdf['year'] = pd.DatetimeIndex(EWR_gdf['date']).year\n",
    "EWR_gdf['data'] = pd.to_numeric(EWR_gdf['data'], errors='coerce')\n",
    "\n",
    "EWR_gdf.to_csv('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Consolidated_AJMR_Dec1.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "pycharm": {
     "name": "#%% Filter and save merged\n"
    }
   },
   "outputs": [],
   "source": [
    "EWR_gdf_filt = EWR_gdf.drop(columns=['field_sample_id', 'tech_task_name', 'lab_id', 'analysis_method', 'final_qualifier',\n",
    "                                     'report_hit_flag', 'Unnamed: 0', 'Lab ID', 'Sample Name ', 'Sampler', 'Date', 'Site',\n",
    "                                     'Control?', 'Latitude', 'Longitude', 'River Kilometer', 'Location ', 'Site', 'GPS',\n",
    "                                     'Core', 'Depth  (cm)', 'Be-7 analysis done', 'Be-7 detected', 'Name', 'Unnamed: 2',\n",
    "                                     'Unnamed: 3', 'Site Number', 'Lat', 'Lon', 'River KM', 'Long', 'Unnamed: 10',\n",
    "                                     'Unnamed: 11', 'Unnamed: 12', 'PARM DESCRIPTION', 'TIME', 'FIELD_NO', 'LIMS_RCODE',\n",
    "                                     'VALQUALIFI', 'ANALMETH', 'Lat Lon Estimate'])\n",
    "\n",
    "\n",
    "# Drop rows with missing locations\n",
    "EWR_gdf_filt = EWR_gdf_filt.loc[~EWR_gdf_filt['geometry'].is_empty, :]\n",
    "\n",
    "EWR_gdf_filt.to_csv('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Consolidated_AJMR_April22.csv')\n",
    "# 'THg ng/g DW ist analysis', 'THg ng/g DW 2nd analysis', 'THg ng/g DW 3rd analysis',\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "pycharm": {
     "name": "#%% Read in centerline shapefile\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "<AxesSubplot:>"
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 1440x1008 with 1 Axes>",
      "image/png": "[encoded data, 1508428 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw'\n",
    "\n",
    "# river_centerline = gp.read_file('Y:/12828.101 English Wabigoon River/03_Data/02_Physical/16_Waterline/Centerline_for_Modelling_merge_UTMZ15.shp')\n",
    "river_centerline = gp.read_file('//srv-ott3.baird.com/Projects/12828.101 English Wabigoon River/03_Data/02_Physical/16_Waterline/Centerline_for_Modelling_UTMZ15.shp')\n",
    "# river_centerline = river_centerline.to_crs(\"EPSG:4326\")\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(20, 14))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "axes.set_xlim(350000, 520000)\n",
    "axes.set_ylim(5510000, 5580000)\n",
    "ctx.add_basemap(axes, source=mapbox, crs='EPSG:32615')\n",
    "\n",
    "river_centerlineExploded = river_centerline.explode(ignore_index=True)\n",
    "river_centerlineExploded.reset_index(inplace=True)\n",
    "\n",
    "river_centerlineExploded.loc[[0], 'geometry'].plot(ax=axes, color='tab:blue')\n",
    "river_centerlineExploded.loc[[1], 'geometry'].plot(ax=axes, color='tab:orange')\n",
    "river_centerlineExploded.loc[[2], 'geometry'].plot(ax=axes, color='tab:green')\n",
    "river_centerlineExploded.loc[[3], 'geometry'].plot(ax=axes, color='tab:red')\n",
    "river_centerlineExploded.loc[[4], 'geometry'].plot(ax=axes, color='tab:purple')\n",
    "river_centerlineExploded.loc[[5], 'geometry'].plot(ax=axes, color='tab:brown')\n",
    "river_centerlineExploded.loc[[6], 'geometry'].plot(ax=axes, color='tab:pink')\n",
    "river_centerlineExploded.loc[[7], 'geometry'].plot(ax=axes, color='tab:cyan')\n",
    "river_centerlineExploded.loc[[8], 'geometry'].plot(ax=axes, color='tab:grey')\n",
    "river_centerlineExploded.loc[[9], 'geometry'].plot(ax=axes, color='tab:olive')\n",
    "\n",
    "tempMulti = river_centerlineExploded.iloc[[5,0,1,2,3,4,6,7,9], 4]\n",
    "# Put the sub-line coordinates into a list of sublists\n",
    "outcoords = [list(i.coords) for i in tempMulti]\n",
    "# Flatten the list of sublists and use it to make a new line\n",
    "river_centerline_merge = geometry.LineString([i for sublist in outcoords for i in sublist])\n",
    "river_centerline_merge_gpd = gp.GeoSeries(river_centerline_merge)\n",
    "\n",
    "river_centerline_merge_gpd.plot(ax=axes, color='k', linestyle='--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "pycharm": {
     "name": "#%% Find distane of each merged datapoint along the river length\n"
    }
   },
   "outputs": [],
   "source": [
    "EWR_gdf_filt_UTM =  EWR_gdf_filt.to_crs(\"EPSG:26915\")\n",
    "EWR_gdf_filt_UTM.reset_index(inplace=True)\n",
    "\n",
    "for i in range(0, EWR_gdf_filt_UTM.shape[0]):\n",
    "    EWR_gdf_filt_UTM.loc[i, 'RiverKM'] = river_centerline_merge.project(\n",
    "                                            EWR_gdf_filt_UTM.loc[i, 'geometry'])  # 201"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "pycharm": {
     "name": "#%% Find and Plot Centerline Distance\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 1440x504 with 2 Axes>",
      "image/png": "[encoded data, 984820 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "river_centerline_merge_gpd2 =\\\n",
    "    gp.GeoDataFrame(geometry=gp.points_from_xy(\n",
    "        river_centerline_merge.xy[0], river_centerline_merge.xy[1], crs=\"EPSG:32615\"))\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(20, 7))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "axes.set_xlim(350000, 520000)\n",
    "axes.set_ylim(5510000, 5580000)\n",
    "ctx.add_basemap(axes, source=mapbox, crs='EPSG:32615')\n",
    "\n",
    "for i in range(0, river_centerline_merge_gpd2.shape[0]):\n",
    "    river_centerline_merge_gpd2.loc[i, 'RiverKM'] = river_centerline_merge.project(\n",
    "                                            river_centerline_merge_gpd2.loc[i, 'geometry'])  # 201\n",
    "\n",
    "rc = river_centerline_merge_gpd2.plot(column='RiverKM', ax=axes, legend=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "pycharm": {
     "name": "#%% Convert Units\n"
    }
   },
   "outputs": [],
   "source": [
    "EWR_gdf_filt_UTM.loc[EWR_gdf_filt_UTM['result_uom']=='UG/G','data']=EWR_gdf_filt_UTM.loc[\n",
    "    EWR_gdf_filt_UTM['result_uom']=='UG/G','data']*1000\n",
    "EWR_gdf_filt_UTM.loc[EWR_gdf_filt_UTM['result_uom']=='UG/G','result_uom']='NG/G'\n",
    "\n",
    "EWR_gdf_filt_UTM.loc[EWR_gdf_filt_UTM['result_uom']=='MG/KG','data']=EWR_gdf_filt_UTM.loc[\n",
    "    EWR_gdf_filt_UTM['result_uom']=='MG/KG','data']*1000\n",
    "EWR_gdf_filt_UTM.loc[EWR_gdf_filt_UTM['result_uom']=='MG/KG','result_uom']='NG/G'\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "pycharm": {
     "name": "#%% Profile Plotting\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 936x504 with 1 Axes>",
      "image/png": "[encoded data, 37340 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(13, 7))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "# plotvars = ['MeHg', 'Bulk MeHg']\n",
    "# plotvars = ['Mercury']\n",
    "# plotvars = ['Methyl mercury']\n",
    "plotvars = ['Mercury', 'Bulk THG', 'Filtered THg', 'Part THg']\n",
    "\n",
    "# plotSources = ['Consolidated', 'Mill', 'Sediment', 'Cores']\n",
    "plotSources = ['WaterCol', 'Parks', 'Wetland']\n",
    "\n",
    "# sourceLabels = ['Consolidated Reed', 'Mill', 'MOECC Surface Sediment', 'Cores']\n",
    "sourceLabels = ['WaterCol', 'Parks', 'Wetland']\n",
    "\n",
    "sourceMarkers = ['o', 's', '^', '*']\n",
    "\n",
    "for pIDX, source in enumerate(plotSources):\n",
    "    sc  = axes.scatter(EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                (EWR_gdf_filt_UTM['Source']==source), 'RiverKM'],\n",
    "                   EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                (EWR_gdf_filt_UTM['Source']==source), 'data'], 35,\n",
    "                                marker=sourceMarkers[pIDX], label=sourceLabels[pIDX])\n",
    "\n",
    "axes.set_xlabel('Distance along river [m]')\n",
    "# axes.set_ylabel('Surface Sediment MeHg Concentration [Ng/g]')\n",
    "# axes.set_ylabel('Surface Sediment Hg Concentration [Ng/g]')\n",
    "axes.set_ylabel('Water Hg Concentration [Ng/L]')\n",
    "\n",
    "# axes.set_ylim(0, 15000)\n",
    "axes.set_ylim(0, 80)\n",
    "plt.legend()\n",
    "\n",
    "plt.rcParams.update({'font.size': 15})\n",
    "plt.show()\n",
    "fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/WaterHg_Prof.png',\n",
    "             bbox_inches='tight', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\arey\\Anaconda3\\envs\\py39\\lib\\site-packages\\geopandas\\plotting.py:637: FutureWarning: 'axes' is deprecated, please use 'ax' instead (for consistency with pandas)\n",
      "  warnings.warn(\n",
      "C:\\Users\\arey\\Anaconda3\\envs\\py39\\lib\\site-packages\\geopandas\\plotting.py:637: FutureWarning: 'axes' is deprecated, please use 'ax' instead (for consistency with pandas)\n",
      "  warnings.warn(\n",
      "C:\\Users\\arey\\Anaconda3\\envs\\py39\\lib\\site-packages\\geopandas\\plotting.py:681: UserWarning: The GeoDataFrame you are attempting to plot is empty. Nothing has been displayed.\n",
      "  warnings.warn(\n",
      "C:\\Users\\arey\\Anaconda3\\envs\\py39\\lib\\site-packages\\geopandas\\plotting.py:637: FutureWarning: 'axes' is deprecated, please use 'ax' instead (for consistency with pandas)\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 1440x504 with 2 Axes>",
      "image/png": "[encoded data, 994688 characters]\n"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(20, 7))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "axes.set_xlim(350000, 520000)\n",
    "axes.set_ylim(5510000, 5580000)\n",
    "ctx.add_basemap(axes, source=mapbox, crs='EPSG:32615')\n",
    "\n",
    "plotSources = ['Consolidated', 'Mill', 'Sediment', 'Cores']\n",
    "sourceLabels = ['Consolidated Reed', 'Mill', 'MOECC Surface Sediment', 'Cores']\n",
    "\n",
    "plotvars = ['Mercury', 'Bulk THG', 'Filtered THg', 'Part THg']\n",
    "plotSources = ['WaterCol', 'Parks', 'Wetland']\n",
    "sourceLabels = ['WaterCol', 'Parks', 'Wetland']\n",
    "\n",
    "sourceMarkers = ['o', 's', '^', '*']\n",
    "# plotvars = ['Mercury']\n",
    "# plotvars = ['Methyl mercury']\n",
    "legendPlot = [True, False, False, False]\n",
    "\n",
    "for pIDX, source in enumerate(plotSources):\n",
    "    sc  = EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                (EWR_gdf_filt_UTM['Source']==source), :].plot(axes=axes,\n",
    "                                marker=sourceMarkers[pIDX], label=sourceLabels[pIDX],\n",
    "                                vmin=0, vmax=30, column='data', legend=legendPlot[pIDX])\n",
    "\n",
    "plt.show()\n",
    "# fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/HgSed_Map.png',\n",
    "#              bbox_inches='tight', dpi=300)\n",
    "# fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/MeHgSed_Map.png',\n",
    "#              bbox_inches='tight', dpi=300)\n",
    "fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/HgWater_Map.png',\n",
    "             bbox_inches='tight', dpi=300)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% Hg Map Plotting\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# sc = axes.scatter(EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                   (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                   (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                   (EWR_gdf_filt_UTM['media'] == 'SED'),'RiverKM'],\n",
    "#             EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                   (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                   (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                   (EWR_gdf_filt_UTM['media'] == 'SED'), 'data'], 35)\n",
    "\n",
    "\n",
    "# sc  = axes.scatter(EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 (EWR_gdf_filt_UTM['Source']=='Consolidated'), 'RiverKM'],\n",
    "#                    EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 (EWR_gdf_filt_UTM['Source']=='Consolidated'), 'data'], 35,\n",
    "#                                 marker='^', label='Consolidated Reed')\n",
    "#\n",
    "# sc  = axes.scatter(EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 (EWR_gdf_filt_UTM['Source']=='Mill'), 'RiverKM'],\n",
    "#                    EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 (EWR_gdf_filt_UTM['Source']=='Mill'), 'data'], 35,\n",
    "#                                 marker='o', label='Mill')\n",
    "#\n",
    "# sc  = axes.scatter(EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 (EWR_gdf_filt_UTM['Source']=='Sediment'), 'RiverKM'],\n",
    "#                    EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 (EWR_gdf_filt_UTM['Source']=='Sediment'), 'data'], 35,\n",
    "#                                 marker='*', label='MOECC Surface Sediment')\n",
    "#\n",
    "# sc  = axes.scatter(EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 (EWR_gdf_filt_UTM['Source']=='Cores'), 'RiverKM'],\n",
    "#                    EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 (EWR_gdf_filt_UTM['Source']=='Cores'), 'data'], 35,\n",
    "#                                 marker='s', label='Cores')\n",
    "\n",
    "# sc  = axes.scatter(EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 ((EWR_gdf_filt_UTM['Source']=='Sediment') |\n",
    "#                                 ((EWR_gdf_filt_UTM['Source']=='Consolidated') & (EWR_gdf_filt_UTM['RiverKM']>100000))), 'RiverKM'],\n",
    "#                    EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "#                                 (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "#                                 (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "#                                 ((EWR_gdf_filt_UTM['Source']=='Sediment') |\n",
    "#                                 ((EWR_gdf_filt_UTM['Source']=='Consolidated') & (EWR_gdf_filt_UTM['RiverKM']>100000))), 'data'], 35,\n",
    "#                                 marker='^', label='Filtered')\n",
    "\n",
    "\n",
    "# axes.set_ylim(0, 1)\n",
    "# cbar = plt.colorbar(sc)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%% Mapping!\n"
    }
   },
   "outputs": [],
   "source": [
    "mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw'\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(27, 14))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "plt.rcParams.update({'font.size': 25})\n",
    "# plotvars = ['Mercury', 'Bulk THG', 'Filtered THg', 'Part THg']\n",
    "plotvars = ['Mercury']\n",
    "\n",
    "pltDat  = EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                (EWR_gdf_filt_UTM['Source']=='Consolidated'), :].plot(column='data',\n",
    "                ax=axes, markersize=30, legend=False, vmin=0, vmax=2000, marker='^')\n",
    "pltDat  = EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                (EWR_gdf_filt_UTM['Source']=='Mill'), :].plot(column='data',\n",
    "                ax=axes, markersize=30, legend=False, vmin=0, vmax=2000, marker='o')\n",
    "pltDat  = EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                (EWR_gdf_filt_UTM['Source']=='Sediment'), :].plot(column='data',\n",
    "                ax=axes, markersize=30, legend=False, vmin=0, vmax=2000, marker='*')\n",
    "pltDat  = EWR_gdf_filt_UTM.loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                (EWR_gdf_filt_UTM['Source']=='Cores'), :].plot(column='data',\n",
    "                ax=axes, markersize=30, legend=True, vmin=0, vmax=2000, marker='s')\n",
    "\n",
    "ctx.add_basemap(axes, source=mapbox, crs='EPSG:26915')\n",
    "axes.set_xlabel('UTM 15N Easting [m]')\n",
    "axes.set_ylabel('UTM 15N Northing [m]')\n",
    "\n",
    "fig.savefig('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/HgSed.png',\n",
    "             bbox_inches='tight', dpi=250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%% Save Surface as XYZ\n"
    }
   },
   "outputs": [],
   "source": [
    "HgSurfaceFilt = np.vstack((np.array(EWR_gdf_filt_UTM.dropna(subset=['data']).loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                ((EWR_gdf_filt_UTM['Source']=='Sediment') |\n",
    "                                ((EWR_gdf_filt_UTM['Source']=='Consolidated') & (EWR_gdf_filt_UTM['RiverKM']>100000))), 'geometry'].x),\n",
    "                    np.array(EWR_gdf_filt_UTM.dropna(subset=['data']).loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                ((EWR_gdf_filt_UTM['Source']=='Sediment') |\n",
    "                                ((EWR_gdf_filt_UTM['Source']=='Consolidated') & (EWR_gdf_filt_UTM['RiverKM']>100000))), 'geometry'].y),\n",
    "                    np.array(EWR_gdf_filt_UTM.dropna(subset=['data']).loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10) &\n",
    "                                ((EWR_gdf_filt_UTM['Source']=='Sediment') |\n",
    "                                ((EWR_gdf_filt_UTM['Source']=='Consolidated') & (EWR_gdf_filt_UTM['RiverKM']>100000))), 'data'])))\n",
    "HgSurfaceFilt = np.transpose(HgSurfaceFilt)\n",
    "\n",
    "np.savetxt('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/SurfaceHg_Filt.xyz',\n",
    "           HgSurfaceFilt, delimiter=\" \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%% Save Filtered Surface as XYZ\n"
    }
   },
   "outputs": [],
   "source": [
    "HgSurface = np.vstack((np.array(EWR_gdf_filt_UTM.dropna(subset=['data']).loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10), 'geometry'].x),\n",
    "                    np.array(EWR_gdf_filt_UTM.dropna(subset=['data']).loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10), 'geometry'].y),\n",
    "                    np.array(EWR_gdf_filt_UTM.dropna(subset=['data']).loc[(EWR_gdf_filt_UTM['param_name'].isin(plotvars)) &\n",
    "                                (EWR_gdf_filt_UTM['date']>'2000-01-01') &\n",
    "                                (abs(EWR_gdf_filt_UTM['depth'])<10), 'data'])))\n",
    "HgSurface = np.transpose(HgSurface)\n",
    "np.savetxt('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/SurfaceHg2.xyz',\n",
    "           HgSurface, delimiter=\" \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%% Interpolate Filtered Surface\n"
    }
   },
   "outputs": [],
   "source": [
    "# Read in bath xyz\n",
    "bathy_xyz = pd.read_csv('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Bed_Level.xyz',\n",
    "                        names=['x', 'y', 'z'], header=0, delim_whitespace=True)\n",
    "\n",
    "rough_xyz = pd.read_csv('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/Roughness.xyz',\n",
    "                        names=['x', 'y', 'z'], header=0, delim_whitespace=True)\n",
    "\n",
    "HgInterp  = sp.interpolate.griddata(HgSurfaceFilt[:,0:1], HgSurfaceFilt[:,2], [bathy_xyz.x, bathy_xyz.y],\n",
    "                                    method='nearest')\n",
    "\n",
    "HgInterpOut = np.vstack((bathy_xyz.x, bathy_xyz.y, HgInterp[0, :]))\n",
    "HgInterpOut = np.transpose(HgInterpOut)\n",
    "\n",
    "A = np.transpose(np.array([rough_xyz.x, rough_xyz.y]))\n",
    "B = np.transpose(np.array([bathy_xyz.x, bathy_xyz.y]))\n",
    "\n",
    "# Set Wetlands to zero\n",
    "# Interpolate Roughness\n",
    "RoughInterp = sp.interpolate.griddata(np.transpose(np.array([rough_xyz.x, rough_xyz.y])), rough_xyz.z, np.transpose(np.array([bathy_xyz.x, bathy_xyz.y])),\n",
    "                                    method='nearest')\n",
    "\n",
    "HgInterpOut[RoughInterp==0.05, 2] = 0\n",
    "\n",
    "np.savetxt('C:/Users/arey/files/Projects/Grassy Narrows/LocalData/SurfaceInterFilt.xyz',\n",
    "           HgInterpOut, delimiter=\" \")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%% Mapping!\n"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw'\n",
    "\n",
    "fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(10, 14))\n",
    "fig.patch.set_facecolor('white')\n",
    "\n",
    "consolidate_gdf.dropna(subset=['X_coord']).plot(ax=axes[00], markersize=12, label='Consolidated Points')\n",
    "mill_gdf.plot(ax=axes[0], markersize=12, label='Mill')\n",
    "core_gdf.plot(ax=axes[0], markersize=12, label='Cores')\n",
    "sed_gdf.plot(ax=axes[0], markersize=12, label='Sediment')\n",
    "waterCol_gdf.plot(ax=axes[0], markersize=12, label='Water Column')\n",
    "wetland_gdf.plot(ax=axes[0], markersize=12, label='Wetland')\n",
    "parks_gdf.plot(ax=axes[0], markersize=12, label='Parks')\n",
    "riverTmpDO_gdf.plot(ax=axes[0], markersize=12, label='River DO')\n",
    "ballTmpDO_gdf.plot(ax=axes[0], markersize=12, label='Ball Lake DO')\n",
    "clayTmpDO_gdf.plot(ax=axes[0], markersize=12, label='Clay Lake DO')\n",
    "\n",
    "ctx.add_basemap(axes[0], source=mapbox, crs='EPSG:4326')\n",
    "axes[0].set_xlabel('UTM 15N Easting [m]')\n",
    "axes[0].set_ylabel('UTM 15N Northing [m]')\n",
    "axes[0].legend()\n",
    "\n",
    "\n",
    "pltDat  = EWR_gdf.plot(column='year', ax=axes[1], markersize=12, legend=True)\n",
    "\n",
    "ctx.add_basemap(axes[1], source=mapbox, crs='EPSG:4326')\n",
    "axes[1].set_xlabel('UTM 15N Easting [m]')\n",
    "axes[1].set_ylabel('UTM 15N Northing [m]')\n",
    "\n",
    "\n",
    "# pltDat  = EWR_gdf.loc[EWR_gdf.param_name=='Mercury', :].plot(column='data',\n",
    "#                 ax=axes[2], markersize=12, legend=True, vmin=0, vmax=500)\n",
    "# pltDat  = consolidate_gdf.dropna(subset=['X_coord']).loc[consolidate_gdf.param_name=='Methyl mercury', :].plot(column='final_result',\n",
    "#                 ax=axes[2], markersize=12, legend=True, vmin=0, vmax=10)\n",
    "pltDat  = EWR_gdf.loc[EWR_gdf.param_name=='DO', :].plot(column='data',\n",
    "                ax=axes[2], markersize=12, legend=True, vmin=0, vmax=20)\n",
    "\n",
    "ctx.add_basemap(axes[2], source=mapbox, crs='EPSG:4326')\n",
    "axes[2].set_xlabel('UTM 15N Easting [m]')\n",
    "axes[2].set_ylabel('UTM 15N Northing [m]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%  Save as pickle\n"
    }
   },
   "outputs": [],
   "source": [
    "filename='C:/Users/arey/files/Projects/Grassy Narrows/HgDat.out'\n",
    "my_shelf = shelve.open(filename,'n') # 'n' for new\n",
    "\n",
    "for key in dir():\n",
    "    try:\n",
    "        my_shelf[key] = globals()[key]\n",
    "    except TypeError:\n",
    "        #\n",
    "        # __builtins__, my_shelf, and imported modules can not be shelved.\n",
    "        #\n",
    "        print('ERROR shelving: {0}'.format(key))\n",
    "my_shelf.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "pycharm-48af278f",
   "language": "python",
   "display_name": "PyCharm (EWR_Data)"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}