some method changes

README.md

@@ -3,7 +3,7 @@
 <a href="http://creativecommons.org/licenses/by/4.0/"><img src="https://shields.io/github/license/romeroqe/mld-mtd" alt="License"></a>
 <a href="https://zenodo.org/badge/latestdoi/524132263"><img src="https://zenodo.org/badge/524132263.svg" alt="DOI"></a>
 
-This is a methodology to locate the minimum and maximum depth of the thermocline and its thickness by fitting the sigmoid function to the temperature profiles in the global ocean. This methodology can be applied to both global and local studies in those areas of the ocean where the water column can be divided into three layers according to its thermal structure.
+This is a methodology to locate the minimum and maximum depth of the thermocline, its thickness, and its strength by fitting the sigmoid function to the temperature profiles in the global ocean. This methodology can be applied to both global and local studies in those areas of the ocean where the water column can be divided into three layers according to its thermal structure.
 
 ## Installation
 To use this methodology, simply download the `mldmtd.py` file.
@@ -12,10 +12,10 @@ To use this methodology, simply download the `mldmtd.py` file.
 To locate the minimum and maximum depth of the thermocline of an Argo float profile::
 
 ```python
-from mldmtd import getProfileData, thermocline
+from mldmtd import getProfileDataFromArgoNc, thermocline
 
-df = getProfileData('.../aoml/4900432/profiles/D4900432_106.nc')
-pres_thermo, temp_thermo, pres_mld, temp_mld, r2, N2, pres_pred, temp_pred = thermocline(df)
+df = getProfileDataFromArgoNc('.../aoml/4900432/profiles/D4900432_106.nc')
+pres_mtd, temp_mtd, pres_mld, temp_mld, r2, N2T, pres_pred, temp_pred = thermocline(df)
 ```
 
 ## License

mldmtd.py

@@ -10,9 +10,24 @@
 from netCDF4 import Dataset
 
 ########################
-# Get data from netcdf
+# Get DataFrame from temperature, salinity and pressure data
 ########################
-def getProfileData(file):
+def getProfileDataFromTSP(temp, psal, pres, lon, lat):
+	n = len(pres)
+	df = pd.DataFrame({'latitude': np.repeat(lat, n), 'longitude': np.repeat(lon, n), 'pres': pres, 'psal': psal, 'temp': temp})
+	df = df.dropna().reset_index(drop=True)
+
+	########################
+	# Calculate absolute salinity and conservative temperature
+	########################
+	df['asal'] = gsw.SA_from_SP(df.psal, df.pres, df.longitude, df.latitude)
+	df['ctemp'] = gsw.CT_from_t(df.asal, df.temp, df.pres)
+	return df
+
+########################
+# Get DataFrame from Argo netcdf
+########################
+def getProfileDataFromArgoNc(file, QC_filter=True):
 	try:
 		nc = Dataset(file)
 		JULD = nc.variables['JULD']
@@ -22,13 +37,22 @@ def getProfileData(file):
 		PRES_ADJUSTED = nc.variables['PRES_ADJUSTED'][0]
 		TEMP_ADJUSTED = nc.variables['TEMP_ADJUSTED'][0]
 		PSAL_ADJUSTED = nc.variables['PSAL_ADJUSTED'][0]
+		PRES_ADJUSTED_QC = [QC.decode('UTF-8') for QC in nc.variables['PRES_ADJUSTED_QC'][0]]
+		TEMP_ADJUSTED_QC = [QC.decode('UTF-8') for QC in nc.variables['TEMP_ADJUSTED_QC'][0]]
+		PSAL_ADJUSTED_QC = [QC.decode('UTF-8') for QC in nc.variables['PSAL_ADJUSTED_QC'][0]]
 		nc.close()
 	except:
 		return pd.DataFrame([])
 
+	# QC flags filter
+	if QC_filter:
+		PRES_ADJUSTED[~np.isin(PRES_ADJUSTED_QC, ["1","2"])] = np.nan
+		TEMP_ADJUSTED[~np.isin(TEMP_ADJUSTED_QC, ["1","2"])] = np.nan
+		PSAL_ADJUSTED[~np.isin(PSAL_ADJUSTED_QC, ["1","2"])] = np.nan
+
 	n = len(PRES_ADJUSTED)
 	df = pd.DataFrame({'date': np.repeat(JULD, n), 'latitude': np.repeat(LATITUDE, n), 'longitude': np.repeat(LONGITUDE, n), 'pres': PRES_ADJUSTED, 'psal': PSAL_ADJUSTED, 'temp': TEMP_ADJUSTED})
-	df = df.dropna()
+	df = df.dropna().reset_index(drop=True)
 
 	########################
 	# Calculate absolute salinity and conservative temperature
@@ -55,6 +79,48 @@ def norm(y, y_min, y_max):
 def unnorm(y, y_min, y_max):
 	return y*(y_max-y_min)+y_min
 
+########################
+# NsquaredT
+########################
+def NsquaredT(SA, CT, p, lat=None, axis=0):
+	# Modified from gsw Nsquared function to get NsquaredT 
+	# (release https://github.com/TEOS-10/GSW-Python/releases/tag/v3.6.16.post1)
+	if lat is not None:
+		if np.any((lat < -90) | (lat > 90)):
+			raise ValueError('lat is out of range')
+		SA, CT, p, lat = np.broadcast_arrays(SA, CT, p, lat)
+		g = grav(lat, p)
+	else:
+		SA, CT, p = np.broadcast_arrays(SA, CT, p)
+		g = 9.7963
+
+	def axis_slicer(n, sl, axis):
+		itup = [slice(None)] * n
+		itup[axis] = sl
+		return tuple(itup)
+
+	db_to_pa = 1e4
+	shallow = axis_slicer(SA.ndim, slice(-1), axis)
+	deep = axis_slicer(SA.ndim, slice(1, None), axis)
+	if lat is not None:
+		g_local = 0.5 * (g[shallow] + g[deep])
+	else:
+		g_local = g
+
+	dSA = SA[deep] - SA[shallow]
+	dCT = CT[deep] - CT[shallow]
+	dp = p[deep] - p[shallow]
+	SA_mid = 0.5 * (SA[shallow] + SA[deep])
+	CT_mid = 0.5 * (CT[shallow] + CT[deep])
+	p_mid = 0.5 * (p[shallow] + p[deep])
+
+	specvol_mid, alpha_mid, beta_mid = gsw.specvol_alpha_beta(SA_mid, CT_mid, p_mid)
+
+	N2T = ((g_local**2) / (specvol_mid * db_to_pa * dp))
+	N2T = N2T * alpha_mid*dCT
+
+	return N2T, p_mid
+
 ########################
 # Calculate thermocline
 ########################
@@ -71,22 +137,22 @@ def thermocline(df, m_precision=0.01, threshold=0.2):
 	########################
 	# Calculates the buoyancy frequency squared (N2)
 	########################
-	N2, p_mid = gsw.Nsquared(df.asal, df.ctemp, df.pres)
-	N2 = np.vstack((N2, p_mid)).T
+	N2T, p_mid = NsquaredT(df.asal, df.ctemp, df.pres)
+	N2T = np.vstack((N2T, p_mid)).T
 
 	########################
-	# Obtain the pressure range where N2 is greater
+	# Obtain the pressure range where N2T is greater
 	########################
 	try:
-		index = np.argmax(N2[:,0])
+		index = np.argmax(np.abs(N2T[:,0]))
 	except:
 		return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan
 
-	max_pres = x_true[index]*2
+	max_pres = p_mid[index]*2
 	index = x_true < max_pres
 	x_true = x_true[index] 
 	y_true = y_true[index]
-	N2 = N2[index[:-1]]
+	N2T = N2T[index[:-1]]
 
 	########################
 	# Normalize data
@@ -133,7 +199,7 @@ def thermocline(df, m_precision=0.01, threshold=0.2):
 
 	pres_mld = np.nan
 	temp_mld = np.nan
-	for k in np.arange(pres_10m + m_precision, x_true[-1], m_precision): # check every 10 cm
+	for k in np.arange(pres_10m + m_precision, x_true[-1], m_precision): # check every cm
 		temp_mld = unnorm(fsigmoid(k, *popt), y_min, y_max)
 		if (temp_mld - temp_10m) >= threshold:
 			pres_mld = k
@@ -145,16 +211,16 @@ def thermocline(df, m_precision=0.01, threshold=0.2):
 	pres_deep = x_true[-1]
 	temp_deep = unnorm(fsigmoid(pres_deep, *popt), y_min, y_max)
 
-	pres_thermo = np.nan
-	temp_thermo = np.nan
+	pres_mtd = np.nan
+	temp_mtd = np.nan
 	if ~np.isnan(pres_mld):
-		for k in np.arange(pres_deep + m_precision, pres_mld, -m_precision): # verificación cada 10 cm
-			temp_thermo = unnorm(fsigmoid(k, *popt), y_min, y_max)
-			if (temp_deep - temp_thermo) >= threshold:
-				pres_thermo = k
+		for k in np.arange(pres_deep + m_precision, pres_mld, -m_precision): # check every cm
+			temp_mtd = unnorm(fsigmoid(k, *popt), y_min, y_max)
+			if (temp_deep - temp_mtd) >= threshold:
+				pres_mtd = k
 				break
 
-	if pres_mld > pres_thermo:
+	if pres_mld > pres_mtd:
 		return np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan
 
-	return pres_thermo, temp_thermo*sign, pres_mld, temp_mld*sign, r2, N2, x_pred, y_pred*sign
+	return pres_mtd, temp_mtd*sign, pres_mld, temp_mld*sign, r2, N2T, x_pred, y_pred*sign