##################################
###### diag_ts_python.py
###### https://oceanpython.org/2013/02/17/t-s-diagram/
######  APetrenko  Jan. 2023
###### to prepare PHYBIO
##################################

import numpy as np
import gsw
import matplotlib.pyplot as plt
import pandas as pd
 
# Extract data from file *********************************
# f = open('CTD_profile.csv', 'r')
f=open('dPHYBIO_2022_ST_K_CTD_12.cnv','rb')

mydata = pd.read_csv(f,
                     header=260,                     
                     sep="\s+",
                     decimal=".",
                     usecols=[13, 16],
                     names=["salt" , "temp"],
                     encoding="ISO-8859-1")

# Figure out boudaries (mins and maxs)
smin = mydata.salt.min() - (0.01 * mydata.salt.min())
smax = mydata.salt.max() + (0.01 * mydata.salt.max())
tmin = mydata.temp.min() - (0.1 * mydata.temp.max())
tmax = mydata.temp.max() + (0.1 * mydata.temp.max())
 
# Calculate how many gridcells we need in the x and y dimensions
xdim = int(round((smax-smin)/0.1+1,0))
ydim = int(round((tmax-tmin)+1,0))
 
# Create empty grid of zeros
dens = np.zeros((ydim,xdim))
 
# Create temp and salt vectors of appropriate dimensions
ti = np.linspace(1,ydim-1,ydim)+tmin
si = np.linspace(1,xdim-1,xdim)*0.1+smin
 
# Loop to fill in grid with densities
for j in range(0,int(ydim)):
    for i in range(0, int(xdim)):
        dens[j,i]=gsw.rho(si[i],ti[j],0)
 
# Substract 1000 to convert to sigma-t
dens = dens - 1000
 
# Plot data ***********************************************
fig1 = plt.figure()
ax1 = fig1.add_subplot(111)
CS = plt.contour(si,ti,dens, linestyles='dashed', colors='k')
plt.clabel(CS, fontsize=12, inline=1, fmt='%0.1f') # Label every second level
 
ax1.plot(mydata.salt,mydata.temp,'-r',markersize=9)
plt.title("PHYBIO_2022_ST_K_CTD_12")
ax1.set_xlabel('Salinity (PSU)')
ax1.set_ylabel('Temperature (°C)')

plt.show()