Add data management scripts
This commit is contained in:
Alexander Rey
parent
785f39de3e
commit
86f881f09f
|
|
@ -0,0 +1,91 @@
|
||||||
|
import pandas as pd
|
||||||
|
import matplotlib as mpl
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
from math import sin, cos, sqrt, atan2, radians
|
||||||
|
import numpy as np
|
||||||
|
import geopandas as gp
|
||||||
|
from matplotlib.tri import Triangulation
|
||||||
|
|
||||||
|
#read in data
|
||||||
|
df = pd.read_csv('NC_CurrentMeter_All_Phase1_all_mobile_2013_05_02.csv')
|
||||||
|
#convert data to CRS in delft model
|
||||||
|
gdf = gp.GeoDataFrame(df, geometry = gp.points_from_xy(df.Longitude,df.Latitude))
|
||||||
|
gdf.crs = {'init':'epsg:4326'}
|
||||||
|
gdf.geometry = gdf.geometry.to_crs({'init':'epsg:32118'})
|
||||||
|
gdf.geometry.crs = (4326)
|
||||||
|
gdf.geometry = gdf.geometry.to_crs(32118)
|
||||||
|
gdf['UTM_X'] = gdf.geometry.x
|
||||||
|
gdf['UTM_Y'] = gdf.geometry.y
|
||||||
|
gdf.to_csv('NC_CurrentMeter_All_Phase1_all_mobile_2013_05_02_EPSG_32118_V2.csv',index=False)
|
||||||
|
#Define variables and columns for calculations
|
||||||
|
df['latR'] = ''
|
||||||
|
df['lonR'] = ''
|
||||||
|
df['dlat']= ''
|
||||||
|
df['dlon']= ''
|
||||||
|
df['a']= ''
|
||||||
|
df['c']= ''
|
||||||
|
df['d']= ''
|
||||||
|
R = 6373 #Radius of Earth
|
||||||
|
#loop through all transects
|
||||||
|
for transect_id in df['transect'].unique():
|
||||||
|
T1 = df.loc[(df['transect']==transect_id)]
|
||||||
|
for i in range(1,len(T1)):#loop to calculate distance between points
|
||||||
|
T1.iloc[0,T1.columns.get_loc('latR')] = radians(T1.iloc[0,9])
|
||||||
|
T1.iloc[0,T1.columns.get_loc('lonR')] = radians(T1.iloc[0,10])
|
||||||
|
T1.iloc[i,T1.columns.get_loc('latR')] = radians(T1.iloc[i,9])#convert lat to radians
|
||||||
|
T1.iloc[i,T1.columns.get_loc('lonR')] = radians(T1.iloc[i,10])#convert lon to radians
|
||||||
|
T1.iloc[i,T1.columns.get_loc('dlat')] = T1.iloc[i,T1.columns.get_loc('latR')]-T1.iloc[i-1,T1.columns.get_loc('latR')]#get dif between lat
|
||||||
|
T1.iloc[i,T1.columns.get_loc('dlon')] = T1.iloc[i,T1.columns.get_loc('lonR')]-T1.iloc[i-1,T1.columns.get_loc('lonR')]#get dif between lon
|
||||||
|
T1.iloc[i,T1.columns.get_loc('a')] = (sin(T1.iloc[i,T1.columns.get_loc('dlat')]/2))**2 + cos(T1.iloc[i,T1.columns.get_loc('latR')]) * cos(T1.iloc[i-1,T1.columns.get_loc('latR')]) * (sin(T1.iloc[i,T1.columns.get_loc('dlon')]/2))**2 #calc a
|
||||||
|
T1.iloc[i,T1.columns.get_loc('c')] = 2 * atan2(sqrt(T1.iloc[i,T1.columns.get_loc('a')]), sqrt(1 - T1.iloc[i,T1.columns.get_loc('a')])) #calc c
|
||||||
|
T1.iloc[i,T1.columns.get_loc('d')] = 1000*R*T1.iloc[i,T1.columns.get_loc('c')] #calc dist between points
|
||||||
|
T1['d']=pd.to_numeric(T1['d'])
|
||||||
|
T1['dist']=T1['d'].cumsum() #get cumulative sum of distances
|
||||||
|
#Column names to distances data for plotting
|
||||||
|
depths1 = T1.columns[11:43] #ADCP sensor depth reading column names to list
|
||||||
|
depths = [float(d.replace('m','')) for d in depths1]
|
||||||
|
V = T1.values[:,11:43].astype(float) #get velocity data all rows columns 11-43
|
||||||
|
T1['date'] = T1['year'].astype(str) + '-' + T1['mon'].astype(str).str.zfill(2) + '-' + T1['day'].astype(str).str.zfill(2)
|
||||||
|
# V = V.flatten() #flatten matrix into list
|
||||||
|
# V = pd.to_numeric(V)
|
||||||
|
# dist = np.repeat(T1['dist'],len(depths1))
|
||||||
|
#Plotting
|
||||||
|
# print(f'Plotting Transect {transect_id} Location {T1.iloc[i,0]}')
|
||||||
|
# plt.scatter(dist,depths,c=V)
|
||||||
|
# plt.title(f"Location {T1.iloc[i,0]} Transect {transect_id}")
|
||||||
|
# plt.xlabel("Distance Along Transect (m)")
|
||||||
|
# plt.ylabel("Height Above Bed (m)")
|
||||||
|
# plt.colorbar()
|
||||||
|
# plt.savefig(f"Location {T1.iloc[i,0]} Transect {transect_id}")
|
||||||
|
# #__________________________________________________________________________________________________________________________
|
||||||
|
#quad contourf
|
||||||
|
from math import ceil
|
||||||
|
|
||||||
|
#x = np.reshape(dist.values, (len(T1), len(depths1)))
|
||||||
|
#y = np.reshape(depths, (len(T1), len(depths1)))
|
||||||
|
#c = np.reshape(V, (len(T1), len(depths1)))
|
||||||
|
|
||||||
|
fig, ax = plt.subplots(figsize=(16,4))
|
||||||
|
ax.set_xlabel('Distance along Transect (m)')
|
||||||
|
ax.set_ylabel('Depth below WSL (m)')
|
||||||
|
ax.grid(linestyle=':')
|
||||||
|
colormap = 'jet'
|
||||||
|
vmin=0
|
||||||
|
vmax=0.5
|
||||||
|
#vmax=ceil(np.nanmax(V))
|
||||||
|
cnt = ax.imshow(V.T, extent=[T1['dist'].min(), T1['dist'].max(), min(depths), max(depths)],
|
||||||
|
aspect='auto', origin='lower', cmap=colormap, vmin=vmin, vmax=vmax)
|
||||||
|
# cnt = ax.contourf(x,y,c, cmap='jet', levels=100, vmin=0, vmax=ceil(np.nanmax(c)))
|
||||||
|
ax.set_xlim(0, ceil(np.nanmax(T1['dist'])/10)*10)
|
||||||
|
ax.set_ylim(ceil(np.nanmax(depths)),0)
|
||||||
|
cm = plt.cm.ScalarMappable(cmap=colormap)
|
||||||
|
cm.set_array(V.T)
|
||||||
|
cm.set_clim(vmin, vmax)
|
||||||
|
cb = fig.colorbar(cm)
|
||||||
|
cb.set_label('Velocity Magnitude (m/s)')
|
||||||
|
cb.ax.tick_params(labelsize=8)
|
||||||
|
cb.set_ticks(np.arange(0, ceil(np.nanmax(V))+0.1, 0.1))
|
||||||
|
ax.set_title(f"Transect {transect_id} : {T1['date'].iloc[0]}")
|
||||||
|
fig.savefig(f"transect_{transect_id}.png", dpi=400, bbox_to_inches='tight', pad_inches=0)
|
||||||
|
plt.close(fig)
|
||||||
|
|
||||||
|
|
@ -0,0 +1,546 @@
|
||||||
|
# -*- coding: utf-8 -*-
|
||||||
|
"""
|
||||||
|
@author: aforsythe
|
||||||
|
|
||||||
|
Copied from "P:/11934.201 Newtown Creek TPP – Privileged and Confidential/05_Analyses/07 ADCP/NC_CurrentMeter_All_Phase1_all_data_2012_05_20/ADCP_Plot_v4.py"
|
||||||
|
_AJMR
|
||||||
|
"""
|
||||||
|
|
||||||
|
|
||||||
|
import pandas as pd
|
||||||
|
import matplotlib as mpl
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
import matplotlib.dates as mdates
|
||||||
|
from math import ceil, isnan
|
||||||
|
import numpy as np
|
||||||
|
import geopandas as gp
|
||||||
|
gp.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'
|
||||||
|
|
||||||
|
import scipy as sp
|
||||||
|
import scipy.ndimage
|
||||||
|
from astropy.convolution import convolve, Gaussian2DKernel, interpolate_replace_nans
|
||||||
|
import cartopy.crs as ccrs
|
||||||
|
import contextily as ctx
|
||||||
|
import os
|
||||||
|
import datetime
|
||||||
|
from shapely.geometry import Point
|
||||||
|
|
||||||
|
#%% read in data
|
||||||
|
dataPath = '//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/05_Analyses/07 ADCP/NC_CurrentMeter_All_Phase1_all_data_2012_05_20/'
|
||||||
|
|
||||||
|
df = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/05 Currents/NC_CurrentMeter_All_Phase1_all_data_2013_07_16/NC_CurrentMeter_All_Phase1_all_mobile_2013_05_02.xlsx',
|
||||||
|
sheet_name='mag_all_mobile')
|
||||||
|
|
||||||
|
#convert data to geodataframe
|
||||||
|
gdf = gp.GeoDataFrame(df, geometry=gp.points_from_xy(df.Longitude, df.Latitude, crs="EPSG:4326"))
|
||||||
|
|
||||||
|
#convert data to CRS of D3D
|
||||||
|
gdf.geometry = gdf.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
# Add new cols to geodataframe
|
||||||
|
gdf['Distance'] = np.zeros([len(gdf), 1])
|
||||||
|
gdf['DistanceSmth'] = np.zeros([len(gdf), 1])
|
||||||
|
gdf['dist'] = np.zeros([len(gdf), 1])
|
||||||
|
gdf['date'] = np.zeros([len(gdf), 1])
|
||||||
|
|
||||||
|
# Set Date
|
||||||
|
gdf['date'] = gdf['year'].astype(str) + '-' + gdf['mon'].astype(str).str.zfill(2) + '-' + gdf['day'].astype(
|
||||||
|
str).str.zfill(2)
|
||||||
|
|
||||||
|
# Column names to distances data for plotting
|
||||||
|
depths1 = gdf.columns[11:43] # ADCP sensor depth reading column names to list
|
||||||
|
depths = np.array([float(d.replace('m', '')) for d in depths1]) * -1
|
||||||
|
|
||||||
|
#loop through all transects
|
||||||
|
|
||||||
|
transectCount = 0
|
||||||
|
transects = gdf['transect'].unique()
|
||||||
|
|
||||||
|
kernel = Gaussian2DKernel(x_stddev=0.25)
|
||||||
|
|
||||||
|
for transect_id in transects[range(221, 241)]:
|
||||||
|
# Select a given trasect
|
||||||
|
tMask = (gdf['transect']==transect_id)
|
||||||
|
# Remove rows without locations
|
||||||
|
tMask[pd.isna(gdf['Latitude'])] = False
|
||||||
|
|
||||||
|
# Find distance betwen points in m
|
||||||
|
gdf.loc[tMask, 'Distance'] = gdf[tMask].distance(gdf[tMask].shift())
|
||||||
|
|
||||||
|
# Many of the points are recorded as being at the same location...
|
||||||
|
# Where the distance is zero, set it to half of the following distance
|
||||||
|
|
||||||
|
# Find zeros and iloc after zeros
|
||||||
|
zeroDist = (gdf['Distance'] == 0) & tMask
|
||||||
|
zeroDistShift = zeroDist.shift()
|
||||||
|
zeroDistShift.iloc[0] = False
|
||||||
|
distShift = gdf['Distance'].shift(-1)
|
||||||
|
distShift.iloc[-1] = False
|
||||||
|
gdf.loc[tMask, 'DistanceSmth'] = gdf['Distance'][tMask]
|
||||||
|
|
||||||
|
# Set zeros to half of following
|
||||||
|
gdf.loc[zeroDist, 'DistanceSmth'] = distShift[zeroDist]/2
|
||||||
|
# Set following to half
|
||||||
|
gdf.loc[zeroDistShift, 'DistanceSmth'] = gdf['Distance']/2
|
||||||
|
|
||||||
|
# Set initial to zero
|
||||||
|
gdf.loc[gdf[tMask].index[0], 'DistanceSmth'] = 0
|
||||||
|
|
||||||
|
# If final location is duplicate, set to previous value
|
||||||
|
if gdf.loc[gdf[tMask].index[-1], 'Distance']==0:
|
||||||
|
gdf.loc[gdf[tMask].index[-1], 'DistanceSmth'] = gdf.loc[gdf[tMask].index[-2], 'DistanceSmth']
|
||||||
|
|
||||||
|
# Get cumulative sum of distances
|
||||||
|
gdf.loc[tMask, 'dist'] = gdf.loc[tMask, 'DistanceSmth'].cumsum()
|
||||||
|
|
||||||
|
# get velocity data all rows columns 11-43
|
||||||
|
V = gdf.values[tMask, 11:43].astype(float)
|
||||||
|
# #__________________________________________________________________________________________________________________________
|
||||||
|
# %% Plotting
|
||||||
|
if transectCount == 0:
|
||||||
|
fig, axes = plt.subplots(nrows=5, ncols=4, figsize=(16, 8))
|
||||||
|
fig.tight_layout(pad=2.5)
|
||||||
|
ax = axes.flat
|
||||||
|
|
||||||
|
colormap = 'jet'
|
||||||
|
vmin=0
|
||||||
|
vmax=0.5
|
||||||
|
|
||||||
|
VS = sp.ndimage.filters.gaussian_filter(V, [1, 1], mode='nearest')
|
||||||
|
|
||||||
|
# vDatEnd = (~np.isnan(V)).cumsum(1).argmax(1).astype(int) + 1
|
||||||
|
# VS = convolve(V, kernel)
|
||||||
|
# for i in range(0, len(VS)):
|
||||||
|
# VS[i, vDatEnd[i]:-1] = np.nan
|
||||||
|
|
||||||
|
pltDat = ax[transectCount].pcolormesh(gdf.loc[tMask, 'dist'], depths, np.transpose(VS),
|
||||||
|
shading='auto', vmin=vmin, vmax=vmax, cmap=colormap)
|
||||||
|
|
||||||
|
cbar = fig.colorbar(pltDat, ax=ax[transectCount],
|
||||||
|
shrink=0.95,aspect=5)
|
||||||
|
cbar.set_label('Magnitude [m/s]')
|
||||||
|
ax[transectCount].set_xlabel('Distance Along Transect [m]')
|
||||||
|
ax[transectCount].set_ylabel('Depth below WSL [m]')
|
||||||
|
stationStart = next((i for i, j in enumerate(tMask) if j), None)
|
||||||
|
ax[transectCount].set_title(gdf.loc[stationStart, 'station'])
|
||||||
|
ax[transectCount].set_ylim(-6, 0)
|
||||||
|
|
||||||
|
transectCount = transectCount + 1
|
||||||
|
|
||||||
|
plt.show()
|
||||||
|
fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/Transects221-241.png',
|
||||||
|
bbox_inches='tight', dpi=300)
|
||||||
|
|
||||||
|
|
||||||
|
# %% Load in moored data
|
||||||
|
df_moored_data = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/05 Currents/NC_CurrentMeter_All_Phase1_all_data_2013_07_16/NC_CurrentMeter_All_Phase1_all_moored_2013_05_20.xlsx',
|
||||||
|
sheet_name='mag_all_moored')
|
||||||
|
|
||||||
|
# Shift col names to put second back in
|
||||||
|
colIN = list(df_moored_data.columns)
|
||||||
|
colOUT = colIN[0:6] + ['second'] + colIN[6:-1]
|
||||||
|
df_moored_data.columns = colOUT
|
||||||
|
|
||||||
|
### Moored current meter locations from here:
|
||||||
|
### file://srv-ott3/Projects/11934.201%20Newtown%20Creek%20TPP%20%E2%80%93%20Privileged%20and%20Confidential/03_Data/01_BkgrdReports/NC_DRAFT_DSR_Submittal_No_3_2013-07-03.pdf
|
||||||
|
### Table 3-1, PDF page 65
|
||||||
|
### Locations change between deployment 1-9 and 10/11
|
||||||
|
|
||||||
|
df_moored = pd.DataFrame(
|
||||||
|
{'depOrig': ['NC083CM', 'NC081CM', 'NC082CM', 'NC086CM', 'EK023CM'],
|
||||||
|
'depCurr': ['NC086CM', 'NC087CM', 'NC088CM', 'NC089CM', 'EK023CM'],
|
||||||
|
'Northing': [208519.95, 206083.24, 203835.10, 201381.55, 200827.33],
|
||||||
|
'Easting': [996198.97, 1000959.45, 1004387.42, 1005283.07, 1004644.90],
|
||||||
|
'bedElev': [-16, -17, -20, -18, -20]})
|
||||||
|
gdf_moored_loc = gp.GeoDataFrame(
|
||||||
|
df_moored, geometry=gp.points_from_xy(df_moored.Easting, df_moored.Northing), crs="EPSG:2263")
|
||||||
|
gdf_moored_loc.geometry = gdf_moored_loc.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
# Loop through Station IDs for deployments 1-9 and assign locations
|
||||||
|
# for d in range(1,10):
|
||||||
|
# depMask = df_moored_data['deployment'] == ('dep' + str(d))
|
||||||
|
# for stationIDX, station in enumerate(df_moored['depOrig']):
|
||||||
|
# stationMask = (df_moored_data['station'] == station) & depMask
|
||||||
|
#
|
||||||
|
# # Assign geography to station plus deployment
|
||||||
|
# df_moored_data.loc[stationMask, 'Northing'] = df_moored.loc[stationIDX, 'Northing']
|
||||||
|
# df_moored_data.loc[stationMask, 'Easting'] = df_moored.loc[stationIDX, 'Easting']
|
||||||
|
|
||||||
|
# Loop through Station IDs for deployments 10-11 and assign locations
|
||||||
|
for d in range(1,12):
|
||||||
|
depMask = df_moored_data['deployment'] == 'dep' + str(d)
|
||||||
|
for stationIDX, station in enumerate(df_moored['depCurr']):
|
||||||
|
stationMask = (df_moored_data['station'] == station) & depMask
|
||||||
|
|
||||||
|
# Assign geography to station plus deployment
|
||||||
|
df_moored_data.loc[stationMask, 'Northing'] = df_moored.loc[stationIDX, 'Northing']
|
||||||
|
df_moored_data.loc[stationMask, 'Easting'] = df_moored.loc[stationIDX, 'Easting']
|
||||||
|
|
||||||
|
# Create geodataframe and convert to USSP
|
||||||
|
gdf_moored = gp.GeoDataFrame(
|
||||||
|
df_moored_data, geometry=gp.points_from_xy(df_moored_data.Easting, df_moored_data.Northing), crs="EPSG:2263")
|
||||||
|
|
||||||
|
#convert data to CRS of D3D
|
||||||
|
gdf_moored.geometry = gdf_moored.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
gdf_moored['date'] = pd.to_datetime(gdf_moored['year'].astype(str) + '-' +
|
||||||
|
gdf_moored['month'].astype(str).str.zfill(2) + '-' +
|
||||||
|
gdf_moored['day'].astype(str).str.zfill(2) + ' ' +
|
||||||
|
gdf_moored['hour'].astype(str).str.zfill(2) + ':' +
|
||||||
|
gdf_moored['minute'].astype(str).str.zfill(2))
|
||||||
|
|
||||||
|
|
||||||
|
# %% Plot moored data
|
||||||
|
# Column names to distances data for plotting
|
||||||
|
depths1 = gdf_moored.columns[8:51] # ADCP sensor depth reading column names to list
|
||||||
|
depths_moored = np.array([float(d.replace('m', '')) for d in depths1])
|
||||||
|
|
||||||
|
colormap = 'jet'
|
||||||
|
vmin = 0
|
||||||
|
vmax = 0.2
|
||||||
|
|
||||||
|
fig, axes = plt.subplots(nrows=5, ncols=1, figsize=(16, 8))
|
||||||
|
fig.tight_layout(pad=2)
|
||||||
|
|
||||||
|
# Loop through Station IDs for deployments 1-9 and assign locations
|
||||||
|
for d in range(1,12):
|
||||||
|
depMask = gdf_moored['deployment'] == ('dep' + str(d))
|
||||||
|
for stationIDX, station in enumerate(df_moored['depCurr']):
|
||||||
|
stationMask = (gdf_moored['station'] == station) & depMask
|
||||||
|
|
||||||
|
V = gdf_moored.values[stationMask, 8:51].astype(float)
|
||||||
|
|
||||||
|
if len(gdf_moored.loc[stationMask, 'date']) != 0:
|
||||||
|
pltDat = axes[stationIDX].pcolormesh(gdf_moored.loc[stationMask, 'date'],
|
||||||
|
depths_moored, np.transpose(V),
|
||||||
|
shading='auto', vmin=vmin, vmax=vmax, cmap=colormap)
|
||||||
|
|
||||||
|
axes[stationIDX].set_ylim(0, 7)
|
||||||
|
fmt_half_year = mdates.MonthLocator(interval=1)
|
||||||
|
axes[stationIDX].xaxis.set_major_locator(fmt_half_year)
|
||||||
|
axes[stationIDX].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
|
||||||
|
|
||||||
|
if d == 1:
|
||||||
|
axes[stationIDX].set_ylabel('Elevation [m]')
|
||||||
|
axes[stationIDX].set_title(station)
|
||||||
|
cbar = fig.colorbar(pltDat, ax=axes[stationIDX],
|
||||||
|
shrink=1.1, aspect=3)
|
||||||
|
cbar.set_label('Magnitude [m/s]')
|
||||||
|
pltDat.set_clim([0, 0.2])
|
||||||
|
|
||||||
|
plt.show()
|
||||||
|
|
||||||
|
fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/ADCP_2012.png',
|
||||||
|
bbox_inches='tight', dpi=300)
|
||||||
|
|
||||||
|
|
||||||
|
# %% Load in ADCP2 data
|
||||||
|
|
||||||
|
# Read in locations
|
||||||
|
df_adcp2_locs = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/05 Currents/NCP2_ADCP_D1-D3/AQ_LocationsSO_ADCP_ADV_overview_Coords_20150126_AJMR.xlsx',
|
||||||
|
sheet_name='ADCP')
|
||||||
|
gdf_adcp2_locs = gp.GeoDataFrame(df_adcp2_locs,
|
||||||
|
geometry=gp.points_from_xy(df_adcp2_locs.X_NYSPLI, df_adcp2_locs.Y_NYSPLI), crs="EPSG:2263")
|
||||||
|
gdf_adcp2_locs = gdf_adcp2_locs.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
|
||||||
|
adcp2_data_path = '//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/05 Currents/'
|
||||||
|
|
||||||
|
adcp2_paths = ['NCP2_ADCP_D1-D3/ADCP_D1_070314_090814', 'NCP2_ADCP_D1-D3/ADCP_D2_090914_110214',
|
||||||
|
'NCP2_ADCP_D1-D3/ADCP_D3_110414_010615', 'NCP2_ADCP_D4-D5/D4_010715_030315',
|
||||||
|
'NCP2_ADCP_D4-D5/D5_030415_050515']
|
||||||
|
|
||||||
|
adcp2_gdfs = dict()
|
||||||
|
|
||||||
|
for depIDX, adcp2_path in enumerate(adcp2_paths):
|
||||||
|
adcp2_files = os.listdir(adcp2_data_path + adcp2_path) # returns list of files in adv folder
|
||||||
|
|
||||||
|
for adcp2_file in adcp2_files:
|
||||||
|
if '.xls' in adcp2_file or '.csv' in adcp2_file:
|
||||||
|
if '.xls' in adcp2_file:
|
||||||
|
df_in = pd.read_excel(adcp2_data_path + adcp2_path + '/' + adcp2_file)
|
||||||
|
else:
|
||||||
|
df_in = pd.read_csv(adcp2_data_path + adcp2_path + '/' + adcp2_file)
|
||||||
|
|
||||||
|
for stationIDX, station in enumerate(df_adcp2_locs['parent_loc_code']):
|
||||||
|
advStrIDX = adcp2_file.find('_')+1
|
||||||
|
if adcp2_file[advStrIDX:advStrIDX+3] in station:
|
||||||
|
adcp2_geo_x = np.ones([len(df_in), 1]) * df_adcp2_locs.X_NYSPLI[stationIDX]
|
||||||
|
adcp2_geo_y = np.ones([len(df_in), 1]) * df_adcp2_locs.Y_NYSPLI[stationIDX]
|
||||||
|
|
||||||
|
|
||||||
|
df_in['date'] = pd.to_datetime(df_in['Year'].astype(str) + '-' + df_in['Month'].astype(
|
||||||
|
str).str.zfill(2) + '-' + df_in['Day'].astype(
|
||||||
|
str).str.zfill(2) + ' ' + df_in['Hour'].astype(
|
||||||
|
str).str.zfill(2) + ':' + df_in['Minute'].astype(
|
||||||
|
str).str.zfill(2) + ':' + df_in['Second'].astype(str).str.zfill(2))
|
||||||
|
|
||||||
|
gdf_in = gp.GeoDataFrame(
|
||||||
|
df_in, geometry=gp.points_from_xy(adcp2_geo_x, adcp2_geo_y), crs="EPSG:2263")
|
||||||
|
|
||||||
|
if adcp2_file[advStrIDX:advStrIDX + 3] not in adcp2_gdfs:
|
||||||
|
adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]] = dict()
|
||||||
|
|
||||||
|
if adcp2_file[0:advStrIDX-1] not in adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]]:
|
||||||
|
adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]][adcp2_file[0:advStrIDX - 1]] = dict()
|
||||||
|
|
||||||
|
if depIDX+1 not in adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]][adcp2_file[0:advStrIDX - 1]]:
|
||||||
|
adcp2_gdfs[adcp2_file[advStrIDX:advStrIDX + 3]][adcp2_file[0:advStrIDX - 1]][depIDX+1] = gdf_in
|
||||||
|
|
||||||
|
# %% Plot ADCP2 Data
|
||||||
|
fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(9, 8))
|
||||||
|
fig.tight_layout(pad=2)
|
||||||
|
|
||||||
|
colormap = 'jet'
|
||||||
|
vmin = 0
|
||||||
|
vmax = 0.2
|
||||||
|
|
||||||
|
for stationIDX, stat in enumerate(adcp2_gdfs):
|
||||||
|
for depIDX, depdat in enumerate(adcp2_gdfs[stat]['mag']):
|
||||||
|
depths1 = adcp2_gdfs[stat]['mag'][depdat].columns[6:-2] # ADCP sensor depth reading column names to list
|
||||||
|
depths_moored = np.array([float(d.replace('m', '')) for d in depths1])
|
||||||
|
|
||||||
|
V = adcp2_gdfs[stat]['mag'][depdat].values[:, 6:-2].astype(float)
|
||||||
|
|
||||||
|
pltDat = axes[stationIDX].pcolormesh(adcp2_gdfs[stat]['mag'][depdat].loc[:, 'date'],
|
||||||
|
depths_moored, np.transpose(V),
|
||||||
|
shading='auto', vmin=vmin, vmax=vmax, cmap=colormap)
|
||||||
|
|
||||||
|
axes[stationIDX].set_ylim(0, 7)
|
||||||
|
axes[stationIDX].set_xlim(pd.to_datetime("2014-07-01"), pd.to_datetime('2015-05-15'))
|
||||||
|
#axes[stationIDX, depIDX].format_xdata = mdates.DateFormatter('%Y-%m')
|
||||||
|
fmt_half_year = mdates.MonthLocator(interval=1)
|
||||||
|
axes[stationIDX].xaxis.set_major_locator(fmt_half_year)
|
||||||
|
axes[stationIDX].xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
|
||||||
|
|
||||||
|
axes[stationIDX].set_title(str(stat))
|
||||||
|
|
||||||
|
axes[stationIDX].set_ylabel('Elevation [m]')
|
||||||
|
|
||||||
|
cbar = fig.colorbar(pltDat, ax=axes[stationIDX],
|
||||||
|
shrink=1.1, aspect=3)
|
||||||
|
cbar.set_label('Magnitude [m/s]')
|
||||||
|
pltDat.set_clim([0, 0.2])
|
||||||
|
|
||||||
|
plt.show()
|
||||||
|
|
||||||
|
|
||||||
|
fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/ADCP_2014.png',
|
||||||
|
bbox_inches='tight', dpi=300)
|
||||||
|
|
||||||
|
# %% Load in Water Level Data
|
||||||
|
# Gauge Locations from map
|
||||||
|
gdf_gaugeLoc = gp.read_file('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/NTC6.kml')
|
||||||
|
gdf_gaugeLocUSSP = gdf_gaugeLoc.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
# All in Feet NAVD88
|
||||||
|
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',
|
||||||
|
sheet_name='Field_Facility_Gauge')
|
||||||
|
gdf_wl_FFG = gp.GeoDataFrame(df_wl_FFG,
|
||||||
|
geometry=gp.points_from_xy(np.ones([len(df_wl_FFG), 1]) * gdf_gaugeLoc['geometry'].x[2],
|
||||||
|
np.ones([len(df_wl_FFG), 1]) * gdf_gaugeLoc['geometry'].y[2]), crs="EPSG:4326")
|
||||||
|
gdf_wl_FFG.geometry = gdf_wl_FFG.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
|
||||||
|
df_wl_NGG1 = 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',
|
||||||
|
sheet_name='National_Grid_Gauge')
|
||||||
|
gdf_wl_NGG1 = gp.GeoDataFrame(df_wl_NGG1,
|
||||||
|
geometry=gp.points_from_xy(np.ones([len(df_wl_NGG1), 1]) * gdf_gaugeLoc['geometry'].x[3],
|
||||||
|
np.ones([len(df_wl_NGG1), 1]) * gdf_gaugeLoc['geometry'].y[3]), crs="EPSG:4326")
|
||||||
|
gdf_wl_NGG1.geometry = gdf_wl_NGG1.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
# Also includes temperature
|
||||||
|
df_wl_NGG2 = pd.read_excel('//srv-ott3/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',
|
||||||
|
sheet_name='Sheet1')
|
||||||
|
gdf_wl_NGG2 = gp.GeoDataFrame(df_wl_NGG2,
|
||||||
|
geometry=gp.points_from_xy(np.ones([len(df_wl_NGG2), 1]) * gdf_gaugeLoc['geometry'].x[3],
|
||||||
|
np.ones([len(df_wl_NGG2), 1]) * gdf_gaugeLoc['geometry'].y[3]), crs="EPSG:4326")
|
||||||
|
gdf_wl_NGG2.geometry = gdf_wl_NGG2.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
# %% Plot Water Level Data
|
||||||
|
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 4))
|
||||||
|
|
||||||
|
axes.plot(gdf_wl_FFG.Date_time,
|
||||||
|
gdf_wl_FFG['Water Surface Elevation (ft)']*0.3048,
|
||||||
|
label='Field Facility Gauge')
|
||||||
|
|
||||||
|
axes.plot(gdf_wl_NGG1.Date_time,
|
||||||
|
gdf_wl_NGG1['Water Surface Elevation (ft)']*0.3048,
|
||||||
|
label='National Grid Gauge Deployment 1')
|
||||||
|
|
||||||
|
axes.plot(gdf_wl_NGG2.datetime,
|
||||||
|
gdf_wl_NGG2['water_surface_elevation']*0.3048,
|
||||||
|
label='National Grid Gauge Deployment 2')
|
||||||
|
|
||||||
|
fmt_half_year = mdates.MonthLocator(interval=6)
|
||||||
|
axes.xaxis.set_major_locator(fmt_half_year)
|
||||||
|
axes.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
|
||||||
|
axes.set_ylabel('Water Surface Elevation [mNAVD88]')
|
||||||
|
axes.set_title('Water Surface Elevation')
|
||||||
|
|
||||||
|
axes.legend()
|
||||||
|
|
||||||
|
fig.show()
|
||||||
|
fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/WaterLevel.png',
|
||||||
|
bbox_inches='tight', dpi=300)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
# %% Load temperature data
|
||||||
|
df_tdat = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/10_Salinity/tData.xlsx',
|
||||||
|
skiprows=[1])
|
||||||
|
|
||||||
|
# Day Month order is reversed
|
||||||
|
df_tdat['date'] = pd.to_datetime(df_tdat['CollectionDate'])
|
||||||
|
|
||||||
|
df_tsample = pd.read_csv('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/10_Salinity/tSampleLocation.csv')
|
||||||
|
|
||||||
|
# Create geodataframe and convert to USSP
|
||||||
|
gdf_tsample = gp.GeoDataFrame(
|
||||||
|
df_tsample, geometry=gp.points_from_xy(df_tsample.EastCoordinate, df_tsample.NorthCoordinate), crs="EPSG:2263")
|
||||||
|
|
||||||
|
#convert data to CRS of D3D
|
||||||
|
gdf_tsample.geometry = gdf_tsample.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
tdat_geo = list()
|
||||||
|
tdat_facility = list()
|
||||||
|
|
||||||
|
for i in range(0, len(df_tdat)):
|
||||||
|
#sampleID = df_tdat['LocationID'][i][0:df_tdat['LocationID'][i].find('_')]
|
||||||
|
#geoFind = df_tsample.loc[df_tsample['ParentLocationID'] == sampleID]
|
||||||
|
geoFind = df_tsample.loc[df_tsample['LocationID'] == df_tdat['LocationID'][i]].geometry.values
|
||||||
|
|
||||||
|
if len(geoFind) !=0:
|
||||||
|
tdat_geo.append(df_tsample.loc[df_tsample['LocationID'] == df_tdat['LocationID'][i]].geometry.values[0])
|
||||||
|
tdat_facility.append(df_tsample.loc[df_tsample['LocationID'] == df_tdat['LocationID'][i]].SourceArea.values[0])
|
||||||
|
else:
|
||||||
|
tdat_geo.append(Point())
|
||||||
|
tdat_facility.append('')
|
||||||
|
|
||||||
|
# Create geodataframe
|
||||||
|
gdf_tdat = gp.GeoDataFrame(df_tdat, geometry=tdat_geo, crs="EPSG:32118")
|
||||||
|
gdf_tdat.loc[:, 'SourceArea'] = tdat_facility
|
||||||
|
|
||||||
|
gdf_tdat.loc[gdf_tdat.loc[:, 'DepthUnit']=='ft', 'DepthM'] = gdf_tdat.loc[gdf_tdat.loc[:, 'DepthUnit']=='ft', 'BeginDepth']*0.3048
|
||||||
|
gdf_tdat.loc[gdf_tdat.loc[:, 'DepthUnit']=='cm', 'DepthM'] = gdf_tdat.loc[gdf_tdat.loc[:, 'DepthUnit']=='cm', 'BeginDepth']*0.01
|
||||||
|
|
||||||
|
|
||||||
|
# Import spring observations
|
||||||
|
df_springDat = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/11_Temperature/NC_Spr2012_Benthic_Water_Quality_FieldData&Observations_20121022.xlsx')
|
||||||
|
|
||||||
|
# Create geodataframe and convert to USSP
|
||||||
|
gdf_springDat = gp.GeoDataFrame(
|
||||||
|
df_springDat, geometry=gp.points_from_xy(df_springDat.x_coord_as_numeric, df_springDat.y_coord_as_numeric), crs="EPSG:2263")
|
||||||
|
gdf_springDat.geometry = gdf_springDat.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
|
||||||
|
# Import Summer observations
|
||||||
|
df_summerDat = pd.read_excel('//srv-ott3/Projects/11934.201 Newtown Creek TPP – Privileged and Confidential/03_Data/02_Physical/11_Temperature/NC_Sum2012_Benthic_Water_Quality_FieldData&Observations_20130205.xlsx')
|
||||||
|
|
||||||
|
# Create geodataframe and convert to USSP
|
||||||
|
gdf_summerDat = gp.GeoDataFrame(
|
||||||
|
df_summerDat, geometry=gp.points_from_xy(df_summerDat.x_coord_as_numeric, df_summerDat.y_coord_as_numeric), crs="EPSG:2263")
|
||||||
|
gdf_summerDat.geometry = gdf_summerDat.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
# Merged
|
||||||
|
df_SpringSummerDat = pd.concat([df_springDat, gdf_summerDat], ignore_index=True)
|
||||||
|
gdf_SpringSummerDat = gp.GeoDataFrame(
|
||||||
|
df_SpringSummerDat, geometry=gp.points_from_xy(df_SpringSummerDat.x_coord_as_numeric, df_SpringSummerDat.y_coord_as_numeric), crs="EPSG:2263")
|
||||||
|
gdf_SpringSummerDat.geometry = gdf_SpringSummerDat.geometry.to_crs("EPSG:32118")
|
||||||
|
|
||||||
|
# %% Plot Salinity Time series
|
||||||
|
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 8))
|
||||||
|
plotMaskStation = (gdf_tdat.SourceArea == 'Newtown Creek') | (gdf_tdat.SourceArea == 'East Branch of Newtown Creek')
|
||||||
|
plotMaskWater = (gdf_tdat.SampleMedium == 'Surface Water')
|
||||||
|
plotMask = plotMaskStation & plotMaskWater
|
||||||
|
|
||||||
|
pltDat = axes.scatter(gdf_tdat.loc[plotMask, 'date'],
|
||||||
|
gdf_tdat.loc[plotMask, 'DepthM'], 10,
|
||||||
|
gdf_tdat.loc[plotMask, 'NumericResult'])
|
||||||
|
axes.set_ylim(10, 0)
|
||||||
|
axes.set_ylabel('Depth below water surface [m]')
|
||||||
|
axes.set_title('Surface Water Salinity Samples')
|
||||||
|
axes.set_xlabel('Date')
|
||||||
|
|
||||||
|
fmt_half_year = mdates.MonthLocator(interval=12)
|
||||||
|
axes.xaxis.set_major_locator(fmt_half_year)
|
||||||
|
axes.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
|
||||||
|
|
||||||
|
cbar = fig.colorbar(pltDat, ax=axes, shrink=0.95)
|
||||||
|
cbar.set_label('Salinity [PSU]')
|
||||||
|
pltDat.set_clim([5, 40])
|
||||||
|
|
||||||
|
fig.show()
|
||||||
|
fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/Salinity.png',
|
||||||
|
bbox_inches='tight', dpi=300)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
# %% Plot Temperature Time series
|
||||||
|
## Additional salinity obs are here!
|
||||||
|
fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(6, 8))
|
||||||
|
fig.tight_layout(pad=2)
|
||||||
|
|
||||||
|
for i, miles in enumerate(np.arange(0, 3, 0.5)):
|
||||||
|
plotMaskDistance = (gdf_SpringSummerDat.miles_from_NC_mouth > miles) & (
|
||||||
|
gdf_SpringSummerDat.miles_from_NC_mouth < miles + 0.5)
|
||||||
|
plotMaskStation = (gdf_SpringSummerDat.subfacility_code == 'Newtown Creek')
|
||||||
|
plotMaskVar = (gdf_SpringSummerDat.chemical_name == 'Temperature (field)') # Temperature (field)'
|
||||||
|
|
||||||
|
plotMask = plotMaskDistance & plotMaskStation & plotMaskVar
|
||||||
|
|
||||||
|
|
||||||
|
pltDat = axes[i].scatter(gdf_SpringSummerDat.loc[plotMask, 'location_start_date'],
|
||||||
|
gdf_SpringSummerDat.loc[plotMask, 'elev']*0.3048, 10,
|
||||||
|
gdf_SpringSummerDat.loc[plotMask, 'result_value'])
|
||||||
|
|
||||||
|
axes[i].set_ylim(-10, 0)
|
||||||
|
axes[i].set_ylabel('Elevation [m]')
|
||||||
|
|
||||||
|
cbar = fig.colorbar(pltDat, ax=axes[i], shrink=0.95)
|
||||||
|
cbar.set_label('Temp [degC]')
|
||||||
|
pltDat.set_clim([5, 30])
|
||||||
|
|
||||||
|
fmt_half_year = mdates.MonthLocator(interval=1)
|
||||||
|
axes[i].xaxis.set_major_locator(fmt_half_year)
|
||||||
|
axes[i].xaxis.set_major_formatter(mdates.DateFormatter('%m-%Y'))
|
||||||
|
|
||||||
|
axes[0].set_title('Surface Water Temperature Samples')
|
||||||
|
axes[i].set_xlabel('Date')
|
||||||
|
|
||||||
|
fig.show()
|
||||||
|
|
||||||
|
|
||||||
|
fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/Temperature.png',
|
||||||
|
bbox_inches='tight', dpi=300)
|
||||||
|
|
||||||
|
# %% Plot Map
|
||||||
|
mapbox = 'https://api.mapbox.com/styles/v1/alexander0042/ckemxgtk51fgp19nybfmdcb1e/tiles/256/{z}/{x}/{y}@2x?access_token=pk.eyJ1IjoiYWxleGFuZGVyMDA0MiIsImEiOiJjazVmdG4zbncwMHY4M2VrcThwZGUzZDFhIn0.w6oDHoo1eCeRlSBpwzwVtw'
|
||||||
|
|
||||||
|
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(8, 8))
|
||||||
|
|
||||||
|
axes.set_xlim(303500, 306500)
|
||||||
|
axes.set_ylim(61000, 63750)
|
||||||
|
gdf.plot(ax=axes, markersize=10, color='blue', label='Mobile ADCP')
|
||||||
|
gdf_SpringSummerDat.plot(ax=axes, markersize=12, color='magenta', label='Temperature & Salinity')
|
||||||
|
gdf_moored_loc.plot(ax=axes, markersize=20, color='red', label='Moored ADCP 2012')
|
||||||
|
gdf_adcp2_locs.plot(ax=axes, markersize=20, color='orange', label='Moored ADCP 2014')
|
||||||
|
gdf_gaugeLocUSSP.loc[2:3, 'geometry'].plot(ax=axes, markersize=20, color='green', label='Water Level Gauge')
|
||||||
|
|
||||||
|
|
||||||
|
ctx.add_basemap(axes, source=mapbox, crs='EPSG:32118')
|
||||||
|
axes.set_xlabel('New York State Plane Easting [m]')
|
||||||
|
axes.set_ylabel('New York State Plane Northing [m]')
|
||||||
|
axes.legend()
|
||||||
|
|
||||||
|
# axes[1].set_xlim(303500, 306500)
|
||||||
|
# axes[1].set_ylim(61000, 63750)
|
||||||
|
# gdf_SpringSummerDat.plot(ax=axes[1], markersize=12, color='magenta', label='Temperature & Salinity')
|
||||||
|
# axes[1].set_xlabel('New York State Plane Easting [m]')
|
||||||
|
# axes[1].legend()
|
||||||
|
# ctx.add_basemap(axes[1], source=mapbox, crs='EPSG:32118')
|
||||||
|
|
||||||
|
fig.show()
|
||||||
|
fig.savefig('C:/Users/arey/files/Projects/Newtown/DataFigs/DataMap.png',
|
||||||
|
bbox_inches='tight', dpi=300)
|
||||||
Loading…
Reference in New Issue