diff --git a/README.md b/README.md index 6f85d80..6ad16b3 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ License DOI -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 diff --git a/mldmtd.py b/mldmtd.py index f277d58..26d325a 100644 --- a/mldmtd.py +++ b/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