import numpy as np                         # NumPy is an array and math library
import matplotlib.pyplot as plt            # Matplotlib is a visualization (plotting) library
import pandas as pd                        # Pandas lets us work with spreadsheet (.csv) data
import xarray as xr
from datetime import datetime, timedelta   # Datetime helps us work with dates and times

# May 2023 R/V Carson cruise data
# filepaths = ['/content/2023051001001_Carkeek.cnv','/content/2023051101001_Carkeek.cnv']
# export_filepaths = ['/content/2023051001001_Carkeek.csv','/content/2023051101001_Carkeek.csv']

# August 2023 R/V Carson MATE cruise data
filepaths = ['/content/20230825ctd01_Carkeek.cnv','/content/20230825ctd02_Carkeek.cnv']
export_filepaths = ['/content/20230825ctd01_Carkeek.csv','/content/20230825ctd02_Carkeek.csv']

casts = []
for idx, filepath in enumerate(filepaths):
  # Identify header line number (= line_idx - 1) and extract column names
  file_object = open(filepath,'r')
  all_lines = file_object.readlines()
  file_object.close()
  header_names = []
  preamble_to_save = []
  for line_idx, line in enumerate(all_lines):
    if line[0] == '*':
      preamble_to_save.append('# ' + line)
    if line[0] == '#':
      preamble_to_save.append(line)
    if '*END*' in line:
      break
    elif ' name' in line:
      header_names.append(line.split('= ')[1].split(':')[0])
  print('First line of data after header:',line_idx, line)
  print('Header column names:',header_names,'\n')

  # Load data
  cnv_data = pd.read_csv(filepath,header=None,names=header_names,
                        skiprows=line_idx+1,delim_whitespace=True)
  # data = pd.read_csv(filepath,comment='$',delim_whitespace=True,na_values='*',
  #                    header=None,names=col_names[2:])

  # Extract approximate cast (note: this only works for a CNV file with a single cast)
  first_sample_idx = np.argmax(np.diff(cnv_data['depSM'].rolling(window=17,center=True).mean()) > 0.03)
  last_sample_idx = np.where(np.diff(cnv_data['depSM'].rolling(window=17,center=True).mean()) > 0.03)[0][-1]
  cast = cnv_data.loc[first_sample_idx:last_sample_idx]
  cast.reset_index(inplace=True,drop=True)

  # De-spike cast data
  neg_spike_idx = cast['depSM'].index[cast['depSM'].diff() < -10].values
  pos_spike_idx = cast['depSM'].index[cast['depSM'].diff() > 10].values
  oxy_spike_idx = cast['sbeox0ML/L'].index[cast['sbeox0ML/L'] < 0.0].values
  all_spike_idx = np.sort(np.concatenate((neg_spike_idx,pos_spike_idx,oxy_spike_idx)))
  if len(all_spike_idx) > 0:
    cast = cast.drop(index=all_spike_idx)

  # Display and save cast data
  display(cast)
  casts.append(cast)

  # Export data and pre-pend with original header preamble
  cast.to_csv(export_filepaths[idx])
  with open(export_filepaths[idx], 'r+') as f:
    content = f.read()
    f.seek(0,0)
    f.write("".join(preamble_to_save) + content)

plt.figure(dpi=400)
plt.plot(cnv_data['depSM'])
plt.xlabel('Sample number')
plt.ylabel('Depth (m)')
plt.title('R/V Carson MATE cruise (8/25/23) - CTD #1')
plt.grid();

plt.plot(np.diff(cnv_data['depSM']),lw=0.5,c='k')

plt.plot(np.diff(cnv_data['depSM'].rolling(window=19,center=True).mean()),lw=0.5,c='k')

plt.plot(casts[1]['depSM'])
plt.xlabel('Sample number')
plt.ylabel('Depth (m)');

plt.plot(casts[0]['t090C'],casts[0]['depSM'],label='Cast #1')
plt.plot(casts[1]['t090C'],casts[1]['depSM'],label='Cast #2')
plt.legend()
plt.gca().invert_yaxis()
plt.grid(alpha=0.5)
plt.xlabel('Temperature (°C)')
plt.ylabel('Depth (m)');

casts[1]['t090C'].mean()

plt.figure(dpi=400)
plt.plot(casts[0]['t090C'],casts[0]['depSM'],label='Cast #1 (morning)')
plt.plot(casts[1]['t090C'],casts[1]['depSM'],label='Cast #2 (afternoon)')
plt.legend()
plt.gca().invert_yaxis()
plt.grid(alpha=0.5)
plt.xlabel('Temperature (°C)')
plt.ylabel('Depth (m)')
plt.title('R/V Carson MATE cruise (8/25/23)');

oxy_bot_dep = [185,130,90,30,20,0]
oxy_bot_val = [6.23,5.56,5.54,6.57,7.09,8.06]

plt.figure(dpi=400)
plt.plot(casts[1]['sbeox0Mg/L'],casts[1]['depSM'],label='Cast #2 (afternoon) - SBE 43 sensor')
plt.scatter(oxy_bot_val,oxy_bot_dep,c='purple',s=30,zorder=2,label='Cast #2 (afternoon) - bottle titrations')
plt.legend()
plt.gca().invert_yaxis()
plt.grid(alpha=0.5)
plt.xlabel('Dissolved oxygen (mg/L)')
plt.ylabel('Depth (m)')
plt.title('R/V Carson MATE cruise (8/25/23)');

chl_bot_dep = [[185,185,185],150,120,90,60,40,30,20,10,0]
chl_bot_val = [[0.612,0.813,0.799],0.468,0.460,0.274,0.345,0.426,1.111,1.791,4.589,4.165]

plt.figure(dpi=400)
plt.plot(casts[1]['flECO-AFL'],casts[1]['depSM'],label='Cast #2 (afternoon) - ECO-AFL/FL sensor')
for idx in range(len(chl_bot_dep)):
  if type(chl_bot_val[idx]) is not list:
    if idx == len(chl_bot_dep)-1: label = 'Cast #2 (afternoon) - bottle sample analyses'
    else:                         label = None
    plt.scatter(chl_bot_val[idx],chl_bot_dep[idx],c='purple',s=30,zorder=2,label=label)
  else:
    # plt.scatter(chl_bot_val[idx],chl_bot_dep[idx],c='purple',s=10,marker='x',zorder=2,label=label)
    plt.errorbar(np.mean(chl_bot_val[idx]),np.mean(chl_bot_dep[idx]),xerr=np.std(chl_bot_val[idx]),
                 marker='o',ms=3,capsize=3,c='purple')
plt.legend()
plt.gca().invert_yaxis()
plt.grid(alpha=0.5)
plt.xlabel('Chlorophyll fluorescence (mg/m$^3$)')
plt.ylabel('Depth (m)')
plt.title('R/V Carson MATE cruise (8/25/23)');

chl_bot_dep[idx]

filepath = '/content/p1750005.nc'

data = xr.open_dataset(filepath)
display(data)

data_to_export = data[['time','pressure','depth','temperature','conductivity','salinity',
                       'sigma_theta','density','latitude','longitude','buoyancy','speed',
                       'glide_angle','horz_speed','vert_speed','polar_heading']].drop(['ctd_time','ctd_depth'])

export_filepath = '/content/20230512_glider.csv'
data_to_export.to_pandas().to_csv(export_filepath)

plt.plot(data['speed'],data['depth'])

plt.scatter(data['ctd_time'],data['ctd_depth'],c=data['buoyancy'])
plt.colorbar()

plt.plot(casts[0]['longitude'],casts[0]['latitude'],label='R/V Rachel Carson cast #1')
plt.plot(casts[1]['longitude'],casts[1]['latitude'],label='R/V Rachel Carson cast #2')
plt.plot(data['longitude'],data['latitude'],label='Glider')

plt.legend()
plt.xlabel('Longitude (°E)')
plt.ylabel('Latitude (°N)');

plt.plot(casts[0]['t090C'],casts[0]['depSM'],label='R/V Carson cast #1')
plt.plot(casts[1]['t090C'],casts[1]['depSM'],label='R/V Carson cast #2')
plt.plot(data['temperature'],data['depth'],label='Glider')
plt.legend()
plt.gca().invert_yaxis()
plt.grid(alpha=0.5)
plt.xlabel('Temperature (°C)')
plt.ylabel('Depth (m)');
plt.ylim([20,0])

plt.plot(casts[0]['sal00'],casts[0]['depSM'],label='R/V Carson cast #1')
plt.plot(casts[1]['sal00'],casts[1]['depSM'],label='R/V Carson cast #2')
plt.plot(data['salinity'],data['depth'],label='Glider')
plt.legend()
plt.gca().invert_yaxis()
plt.grid(alpha=0.5)
plt.xlabel('Salinity (PSU)')
plt.ylabel('Depth (m)');
plt.ylim([175,130])

plt.plot(casts[0]['sal00'],casts[0]['depSM'],label='R/V Carson cast #1')
plt.plot(casts[1]['sal00'],casts[1]['depSM'],label='R/V Carson cast #2')
plt.plot(data['salinity'],data['depth'],label='Glider')
plt.legend()
plt.gca().invert_yaxis()
plt.grid(alpha=0.5)
plt.xlabel('Salinity (PSU)')
plt.ylabel('Depth (m)');
plt.ylim([20,0])

casts[0]['t090C'].mean()

casts[1]['t090C'].mean()

data['temperature'].mean()