diff --git a/.gitignore b/.gitignore index 571fe028..18ef39c1 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,6 @@ makefile *.grib /.pytest_cache/ .hypothesis/* +bitsea-copilot-coder-bundle.zip +docs +*pyc diff --git a/src/bitsea/Float/CORIOLIS_checks.py b/src/bitsea/Float/CORIOLIS_checks.py index fd8b822a..0cb6c066 100644 --- a/src/bitsea/Float/CORIOLIS_checks.py +++ b/src/bitsea/Float/CORIOLIS_checks.py @@ -30,16 +30,26 @@ def check_lenght_timeseries(df): days = duration.days return(days>365) -# check over nans +# check over nans def nans_check(df, namecol): - """ input : dataframe & nameof column selected for the check + """ output : Boolean method : fill nans inside the dataset to remove the starting part of time series: create list_index nan are counted in the core of time series df_ix[list_index] - """ + """ df[namecol] = pd.to_numeric(df[namecol]) - list_index = list((df.interpolate()).dropna().index) - return(df.iloc[list_index].VAR.isna().sum() < df.shape[0]/2) + list_index = list((df.VAR.interpolate()).dropna().index) + sliced_time_df = df.iloc[list_index] + if sliced_time_df.empty: + return False + if len(sliced_time_df) <5: + return False + sliced_time_df["time"] = pd.to_datetime(sliced_time_df["time"]) + duration_days = (sliced_time_df["time"].iloc[-1] - sliced_time_df["time"].iloc[0]).days + + if duration_days > 365 : return True + else: return False + def lenght_timeseries(df, namecol): df1=df.copy() diff --git a/src/bitsea/Float/Float_opt/superfloat_cdom.py b/src/bitsea/Float/Float_opt/superfloat_cdom.py index 202e2a26..9ad2f4dd 100644 --- a/src/bitsea/Float/Float_opt/superfloat_cdom.py +++ b/src/bitsea/Float/Float_opt/superfloat_cdom.py @@ -1,4 +1,12 @@ import argparse +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of cdom. @@ -47,8 +55,7 @@ def dump_cdom_file(outfile, p, Pres, Value, Qc, metatata, mode='w'): ncOUT = NC.netcdf_file(outfile + ".tmp" ,mode) if mode=='w': # if not existing file, we'll put header, TEMP, PSAL - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, Pres, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/Float_opt/superfloat_chla.py b/src/bitsea/Float/Float_opt/superfloat_chla.py index e1c0ffd0..0fa27a47 100644 --- a/src/bitsea/Float/Float_opt/superfloat_chla.py +++ b/src/bitsea/Float/Float_opt/superfloat_chla.py @@ -1,4 +1,12 @@ import argparse +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of chla. @@ -47,8 +55,7 @@ def dump_chla_file(outfile, p, Pres, Value, Qc, metatata, mode='w'): ncOUT = NC.netcdf_file(outfile + ".tmp" ,mode) if mode=='w': # if not existing file, we'll put header, TEMP, PSAL - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, Pres, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/Float_opt/superfloat_irr380.py b/src/bitsea/Float/Float_opt/superfloat_irr380.py index 6b2db63e..0dd3a468 100644 --- a/src/bitsea/Float/Float_opt/superfloat_irr380.py +++ b/src/bitsea/Float/Float_opt/superfloat_irr380.py @@ -1,4 +1,12 @@ import argparse +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of irr_380. @@ -47,8 +55,7 @@ def dump_irr_380_file(outfile, p, Pres, Value, Qc, metatata, mode='w'): ncOUT = NC.netcdf_file(outfile + ".tmp" ,mode) if mode=='w': # if not existing file, we'll put header, TEMP, PSAL - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, Pres, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/Float_opt/superfloat_irr412.py b/src/bitsea/Float/Float_opt/superfloat_irr412.py index 712214fb..2b0c963f 100644 --- a/src/bitsea/Float/Float_opt/superfloat_irr412.py +++ b/src/bitsea/Float/Float_opt/superfloat_irr412.py @@ -1,4 +1,12 @@ import argparse +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of irr_412. @@ -47,8 +55,7 @@ def dump_irr_412_file(outfile, p, Pres, Value, Qc, metatata, mode='w'): ncOUT = NC.netcdf_file(outfile + ".tmp" ,mode) if mode=='w': # if not existing file, we'll put header, TEMP, PSAL - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, Pres, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/Float_opt/superfloat_irr490.py b/src/bitsea/Float/Float_opt/superfloat_irr490.py index 2804b36e..a49f4445 100644 --- a/src/bitsea/Float/Float_opt/superfloat_irr490.py +++ b/src/bitsea/Float/Float_opt/superfloat_irr490.py @@ -1,4 +1,12 @@ import argparse +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of irr_490. @@ -47,8 +55,7 @@ def dump_irr_490_file(outfile, p, Pres, Value, Qc, metatata, mode='w'): ncOUT = NC.netcdf_file(outfile + ".tmp" ,mode) if mode=='w': # if not existing file, we'll put header, TEMP, PSAL - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, Pres, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/Float_opt/superfloat_nitrate.py b/src/bitsea/Float/Float_opt/superfloat_nitrate.py index 4e15067f..f1042a7f 100644 --- a/src/bitsea/Float/Float_opt/superfloat_nitrate.py +++ b/src/bitsea/Float/Float_opt/superfloat_nitrate.py @@ -1,4 +1,12 @@ import argparse +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of chla. @@ -51,8 +59,7 @@ def dump_nitrate_file(outfile, p_pos, p, Pres, Value, Qc, metadata,mode='w'): os.system(command) ncOUT = NC.netcdf_file(outfile + ".tmp",mode) if mode=='w': # if not existing file, we'll put header, TEMP, PSAL - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, Pres, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/Float_opt/superfloat_oxygen.py b/src/bitsea/Float/Float_opt/superfloat_oxygen.py index e9776683..8091ee26 100644 --- a/src/bitsea/Float/Float_opt/superfloat_oxygen.py +++ b/src/bitsea/Float/Float_opt/superfloat_oxygen.py @@ -39,6 +39,14 @@ def argument(): import numpy as np import gsw +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def remove_bad_sensors(Profilelist,var): ''' @@ -70,8 +78,7 @@ def convert_oxygen(p,doxypres,doxyprofile): from micromol/Kg to mmol/m3 ''' if doxypres.size == 0: return doxyprofile - Pres, temp, Qc = p.read("TEMP",read_adjusted=False) - Pres, sali, Qc = p.read("PSAL",read_adjusted=False) + PresT, temp, QcT, Pres, sali, QcS = read_temp_psal(p) SA = gsw.SA_from_SP(sali, Pres, p.lon, p.lat) density = gsw.rho(SA, gsw.CT_from_t(SA, temp, Pres), Pres) density_on_zdoxy = np.interp(doxypres,Pres,density) @@ -85,8 +92,7 @@ def dump_oxygen_file(outfile, p, Pres, Value, Qc, metatata, mode='w'): ncOUT = NC.netcdf_file(outfile + ".tmp" ,mode) if mode=='w': # if not existing file, we'll put header, TEMP, PSAL - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, Pres, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/Float_opt/superfloat_par.py b/src/bitsea/Float/Float_opt/superfloat_par.py index 72f2cc68..33df5007 100644 --- a/src/bitsea/Float/Float_opt/superfloat_par.py +++ b/src/bitsea/Float/Float_opt/superfloat_par.py @@ -1,4 +1,12 @@ import argparse +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of par. @@ -47,8 +55,7 @@ def dump_par_file(outfile, p, Pres, Value, Qc, metatata, mode='w'): ncOUT = NC.netcdf_file(outfile + ".tmp" ,mode) if mode=='w': # if not existing file, we'll put header, TEMP, PSAL - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, Pres, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/Float_opt/superfloat_ph.py b/src/bitsea/Float/Float_opt/superfloat_ph.py index 9b765b02..a65c4acf 100644 --- a/src/bitsea/Float/Float_opt/superfloat_ph.py +++ b/src/bitsea/Float/Float_opt/superfloat_ph.py @@ -1,4 +1,12 @@ import argparse +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of ph. @@ -47,8 +55,7 @@ def dump_ph_file(outfile, p, Pres, Value, Qc, metatata, mode='w'): ncOUT = NC.netcdf_file(outfile + ".tmp" ,mode) if mode=='w': # if not existing file, we'll put header, TEMP, PSAL - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, Pres, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/MedBGCins.py b/src/bitsea/Float/MedBGCins.py index 636b8a4a..b29867ae 100644 --- a/src/bitsea/Float/MedBGCins.py +++ b/src/bitsea/Float/MedBGCins.py @@ -7,7 +7,7 @@ from bitsea.commons.grid import RegularGrid from bitsea.mhelpers.linear_shift import linear_shift -default_clim_n3n_nc = "/gss/gss_work/DRES_OGS_BiGe/Observations/CLIMATOLOGY/Nutrients/MedBGC-ins/Climatology_N3n_600_800m.nc" +default_clim_n3n_nc = "/g100_work/OGS_test2528/Observations/CLIMATOLOGY/Nutrients/MedBGC-ins/Climatology_N3n_600_800m.nc" clim_file = Path(os.getenv("CLIM_MED_BGC_INS_FILE", default_clim_n3n_nc)) if not clim_file.exists(): diff --git a/src/bitsea/Float/Superfloat_oxy/superfloat_generator.py b/src/bitsea/Float/Superfloat_oxy/superfloat_generator.py index b2437b76..e12c63ff 100644 --- a/src/bitsea/Float/Superfloat_oxy/superfloat_generator.py +++ b/src/bitsea/Float/Superfloat_oxy/superfloat_generator.py @@ -4,6 +4,13 @@ import os +def read_temp_adjusted(profile_obj): + PresT, Temp, QcT = profile_obj.read('TEMP', read_adjusted=True) + if PresT is None or Temp is None or len(PresT) == 0 or len(Temp) == 0: + PresT, Temp, QcT = profile_obj.read('TEMP', read_adjusted=False) + return PresT, Temp, QcT + + def general_quenching(profile_obj, Pres, Profile, Qc): ''' Manages the quenching for Coriolis profiles @@ -14,7 +21,7 @@ def general_quenching(profile_obj, Pres, Profile, Qc): In case of cutted profile (Z_top > 10.) the quench value is replied up to the sea surface on a set of extrapolation points, every 5 meters. ''' - PresT, Temp, _ = profile_obj.read('TEMP', read_adjusted=False) + PresT, Temp, _ = read_temp_adjusted(profile_obj) mld = calculated_depths.mld(Temp, PresT, zref=0, deltaTemp=0.1) if mld < 80: if Pres[0] <= mld+5: @@ -57,7 +64,7 @@ def quenching(profile_obj, PresChl, Chl, chl_lov_zero): ''' Quenched_Chl = Chl.copy() if profile_obj.time.month in range(1,13): - PresT, Temp, _ = profile_obj.read('TEMP', read_adjusted=False) + PresT, Temp, _ = read_temp_adjusted(profile_obj) mld = calculated_depths.mld(Temp, PresT, zref=0, deltaTemp=0.1) if mld > 200:mld=20 diff --git a/src/bitsea/Float/Superfloat_oxy/superfloat_oxy_function.py b/src/bitsea/Float/Superfloat_oxy/superfloat_oxy_function.py index d73d3aec..87377af8 100644 --- a/src/bitsea/Float/Superfloat_oxy/superfloat_oxy_function.py +++ b/src/bitsea/Float/Superfloat_oxy/superfloat_oxy_function.py @@ -8,6 +8,14 @@ import numpy as np import datetime +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def remove_bad_sensors(Profilelist,var): ''' @@ -44,8 +52,7 @@ def dump_oxygen_file(outfile, p, Pres, Value, Qc, metadata, mode='a'): if mode=='w': # if not existing file, we'll put header, TEMP, PSAL setattr(ncOUT, 'origin' , 'coriolis') setattr(ncOUT, 'file_origin', metadata.filename) - PresT, Temp, QcT = p.read('TEMP') #, read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL') # , read_adjusted=False) + PresT, Temp, QcT, PresS, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) diff --git a/src/bitsea/Float/TREND_ANALYSIS.py b/src/bitsea/Float/TREND_ANALYSIS.py index f4405d32..55e66511 100644 --- a/src/bitsea/Float/TREND_ANALYSIS.py +++ b/src/bitsea/Float/TREND_ANALYSIS.py @@ -48,7 +48,7 @@ def compute_trend(df): -def trend_conditions(WMO, days, Bool, DEPTH, min_d, max_d , TI_3, tmp): +def trend_conditions(WMO,cycle, days, Bool, DEPTH, min_d, max_d , TI_3, tmp): """ compute the trend if time-nans and vals for the selected day are satisfied. CASE_* is a Boolean @@ -60,11 +60,11 @@ def trend_conditions(WMO, days, Bool, DEPTH, min_d, max_d , TI_3, tmp): y = np.array(tmp.VAR) x = np.array(tmp.index) LIST_APP_TRENDS = Trend_RANSAC_ThSEN(x , y) - lst = [WMO ,DEPTH, days , min_d , max_d, LIST_APP_TRENDS[0][1], LIST_APP_TRENDS[1][1] ,np.nan,np.nan,np.nan] + lst = [WMO,cycle ,DEPTH, days , min_d , max_d, LIST_APP_TRENDS[0][1], LIST_APP_TRENDS[1][1] ,np.nan,np.nan,np.nan] elif CASE_2: - lst = [WMO ,DEPTH, days , np.nan , np.nan , np.nan , np.nan, np.nan,np.nan,np.nan] + lst = [WMO,cycle ,DEPTH, days, np.nan, np.nan, np.nan, np.nan, np.nan,np.nan,np.nan] else: - lst = [WMO ,DEPTH, days , min_d , max_d, np.nan , np.nan, np.nan,np.nan,np.nan] + lst = [WMO,cycle ,DEPTH, days, min_d, max_d, np.nan, np.nan, np.nan,np.nan,np.nan] return (lst) def drift_coding(WMO,Bool ,df0,df1): diff --git a/src/bitsea/Float/oxygen_saturation.py b/src/bitsea/Float/oxygen_saturation.py index a749e2b5..bc33f72b 100644 --- a/src/bitsea/Float/oxygen_saturation.py +++ b/src/bitsea/Float/oxygen_saturation.py @@ -14,13 +14,16 @@ def oxy_sat(p): * O2o * a concentration of oxygen in mmol/m3 nan if there are no TEMP and PSAL values with pressure less than 5m. ''' - if p.has_adjusted: - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False)#p._my_float.adjusted('TEMP')) - PresS, Sali, QcS = p.read('PSAL', read_adjusted=False)#p._my_float.adjusted('PSAL')) + try: + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + PresS, Sali, QcS = p.read('PSAL', read_adjusted=True) + except TypeError: # dovrebbe entrare se uso superfloat wrapper + PresT, Temp, QcT = p.read('TEMP') + PresS, Sali, QcS = p.read('PSAL') else: - PresT, Temp, QcT = p.read('TEMP') - PresS, Sali, QcS = p.read('PSAL') + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + PresS, Sali, QcS = p.read('PSAL', read_adjusted=False) ii = PresT<=depth_threshold if ii.sum() == 0: diff --git a/src/bitsea/Float/superfloat_bbp700.py b/src/bitsea/Float/superfloat_bbp700.py index 1f8a5d1e..b4c240ee 100644 --- a/src/bitsea/Float/superfloat_bbp700.py +++ b/src/bitsea/Float/superfloat_bbp700.py @@ -1,5 +1,13 @@ import argparse from bitsea.utilities.argparse_types import existing_dir_path +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + PresS, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (PresS is None or PresT is None or Temp is None or Sali is None or len(PresS) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + PresS, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, PresS, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of BBP700. @@ -65,8 +73,7 @@ def dump_bbp700_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): if mode=='w':# if not existing file, we'll put header, TEMP, PSAL setattr(ncOUT, 'origin' , 'coriolis') setattr(ncOUT, 'file_origin', metadata.filename) - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, PresS, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) @@ -90,13 +97,13 @@ def dump_bbp700_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): ncvar=ncOUT.createVariable('TEMP_QC','f',('nTEMP',)) ncvar[:]=QcT - ncvar=ncOUT.createVariable('PSAL','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PSAL','f',('nPSAL',)) ncvar[:]=Sali setattr(ncvar, 'variable' , 'SALI') setattr(ncvar, 'units' , "PSS78") - ncvar=ncOUT.createVariable('PRES_PSAL','f',('nTEMP',)) - ncvar[:]=PresT - ncvar=ncOUT.createVariable('PSAL_QC','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PRES_PSAL','f',('nPSAL',)) + ncvar[:]=PresS + ncvar=ncOUT.createVariable('PSAL_QC','f',('nPSAL',)) ncvar[:]=QcS print("dumping bbp700 on " + str(outfile), flush=True) @@ -130,7 +137,7 @@ def get_outfile(p,outdir): return filename def bbp_algorithm(pCor, outfile, metadata,writing_mode): - Pres, _, _ = pCor.read('TEMP', read_adjusted=False) + Pres, _, _ = read_temp_psal(pCor)[3:] if len(Pres)<5: print("few values in Coriolis TEMP in " + str(pCor._my_float.filename), flush=True) return diff --git a/src/bitsea/Float/superfloat_chla.py b/src/bitsea/Float/superfloat_chla.py index 54296b86..c0d507e7 100644 --- a/src/bitsea/Float/superfloat_chla.py +++ b/src/bitsea/Float/superfloat_chla.py @@ -1,5 +1,13 @@ from bitsea.utilities.argparse_types import existing_dir_path import argparse +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of chla. @@ -65,8 +73,7 @@ def get_outfile(p,outdir): return filename def dumpfile(outfile, p,Pres,chl_profile,Qc,metadata): - _ , Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, PresS, Sali, QcS = read_temp_psal(p) print("dumping chla on " , outfile, p.time.strftime(" %Y%m%d-%H:%M:%S"), flush=True) ncOUT = NC.Dataset(outfile,"w") @@ -77,7 +84,7 @@ def dumpfile(outfile, p,Pres,chl_profile,Qc,metadata): ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) - ncOUT.createDimension('nPSAL', len(PresT)) + ncOUT.createDimension('nPSAL', len(PresS)) ncOUT.createDimension('nCHLA', len(Pres )) ncvar=ncOUT.createVariable("REFERENCE_DATE_TIME", 'c', ("DATETIME",)) @@ -106,7 +113,7 @@ def dumpfile(outfile, p,Pres,chl_profile,Qc,metadata): setattr(ncvar, 'units' , "PSS78") ncvar=ncOUT.createVariable('PRES_PSAL','f',('nPSAL',)) - ncvar[:]=PresT + ncvar[:]=PresS ncvar=ncOUT.createVariable('PSAL_QC','f',('nPSAL',)) ncvar[:]=QcS @@ -137,7 +144,7 @@ def chla_algorithm(pCor,outfile): ''' background_value=0.005 - Pres, _, _ = pCor.read('TEMP', read_adjusted=False) + Pres, _, _ = read_temp_psal(pCor)[3:] if len(Pres)<5: print("few values in Coriolis TEMP in " + str(pCor._my_float.filename), flush=True) return diff --git a/src/bitsea/Float/superfloat_generator.py b/src/bitsea/Float/superfloat_generator.py index 343f30ae..d65c2042 100644 --- a/src/bitsea/Float/superfloat_generator.py +++ b/src/bitsea/Float/superfloat_generator.py @@ -4,6 +4,13 @@ import os +def read_temp_adjusted(profile_obj): + PresT, Temp, QcT = profile_obj.read('TEMP', read_adjusted=True) + if PresT is None or Temp is None or len(PresT) == 0 or len(Temp) == 0: + PresT, Temp, QcT = profile_obj.read('TEMP', read_adjusted=False) + return PresT, Temp, QcT + + def general_quenching(profile_obj, Pres, Profile, Qc): ''' Manages the quenching for Coriolis profiles @@ -14,7 +21,7 @@ def general_quenching(profile_obj, Pres, Profile, Qc): In case of cutted profile (Z_top > 10.) the quench value is replied up to the sea surface on a set of extrapolation points, every 5 meters. ''' - PresT, Temp, _ = profile_obj.read('TEMP', read_adjusted=False) + PresT, Temp, _ = read_temp_adjusted(profile_obj) mld = calculated_depths.mld(Temp, PresT, zref=0, deltaTemp=0.1) if mld < 80: if Pres[0] <= mld+5: @@ -57,7 +64,7 @@ def quenching(profile_obj, PresChl, Chl, chl_lov_zero): ''' Quenched_Chl = Chl.copy() if profile_obj.time.month in range(1,13): - PresT, Temp, _ = profile_obj.read('TEMP', read_adjusted=False) + PresT, Temp, _ = read_temp_adjusted(profile_obj) mld = calculated_depths.mld(Temp, PresT, zref=0, deltaTemp=0.1) if mld > 200:mld=20 diff --git a/src/bitsea/Float/superfloat_kd490.py b/src/bitsea/Float/superfloat_kd490.py index 61250bb5..5375e4ba 100644 --- a/src/bitsea/Float/superfloat_kd490.py +++ b/src/bitsea/Float/superfloat_kd490.py @@ -1,5 +1,13 @@ import argparse from bitsea.utilities.argparse_types import existing_dir_path +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + PresS, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (PresS is None or PresT is None or Temp is None or Sali is None or len(PresS) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + PresS, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, PresS, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of down irradiance at 490nm. @@ -68,8 +76,7 @@ def dump_kd_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): if mode=='w': # if not existing file, we'll put header, TEMP, PSAL setattr(ncOUT, 'origin' , 'coriolis') setattr(ncOUT, 'file_origin', metadata.filename) - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, PresS, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) @@ -93,13 +100,13 @@ def dump_kd_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): ncvar=ncOUT.createVariable('TEMP_QC','f',('nTEMP',)) ncvar[:]=QcT - ncvar=ncOUT.createVariable('PSAL','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PSAL','f',('nPSAL',)) ncvar[:]=Sali setattr(ncvar, 'variable' , 'SALI') setattr(ncvar, 'units' , "PSS78") - ncvar=ncOUT.createVariable('PRES_PSAL','f',('nTEMP',)) - ncvar[:]=PresT - ncvar=ncOUT.createVariable('PSAL_QC','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PRES_PSAL','f',('nPSAL',)) + ncvar[:]=PresS + ncvar=ncOUT.createVariable('PSAL_QC','f',('nPSAL',)) ncvar[:]=QcS print("dumping kd490 on " + str(outfile), flush=True) diff --git a/src/bitsea/Float/superfloat_nitrate.py b/src/bitsea/Float/superfloat_nitrate.py index 1d3d60ee..506bd637 100644 --- a/src/bitsea/Float/superfloat_nitrate.py +++ b/src/bitsea/Float/superfloat_nitrate.py @@ -1,5 +1,13 @@ import argparse from bitsea.utilities.argparse_types import existing_dir_path +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of nitrate. @@ -71,8 +79,7 @@ def convert_nitrate(p,pres,profile): from micromol/Kg to mmol/m3 ''' if pres.size == 0: return profile - _ , temp, _ = p.read("TEMP",read_adjusted=False) - Pres, sali, _ = p.read("PSAL",read_adjusted=False) + _, temp, _, Pres, sali, _ = read_temp_psal(p) SA = gsw.SA_from_SP(sali, Pres, p.lon, p.lat) density = gsw.rho(SA, gsw.CT_from_t(SA, temp, Pres), Pres) density_on_z = np.interp(pres,Pres,density) @@ -89,12 +96,11 @@ def dump_nitrate_file(outfile, p, Pres, Value, Qc, metadata,mode='w'): if mode=='w': # if not existing file, we'll put header, TEMP, PSAL setattr(ncOUT, 'origin' , 'coriolis') setattr(ncOUT, 'file_origin', metadata.filename) - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, PresS, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) - ncOUT.createDimension('nPSAL', len(PresT)) + ncOUT.createDimension('nPSAL', len(PresS)) ncvar=ncOUT.createVariable("REFERENCE_DATE_TIME", 'c', ("DATETIME",)) ncvar[:]=p.time.strftime("%Y%m%d%H%M%S") @@ -114,13 +120,13 @@ def dump_nitrate_file(outfile, p, Pres, Value, Qc, metadata,mode='w'): ncvar=ncOUT.createVariable('TEMP_QC','f',('nTEMP',)) ncvar[:]=QcT - ncvar=ncOUT.createVariable('PSAL','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PSAL','f',('nPSAL',)) ncvar[:]=Sali setattr(ncvar, 'variable' , 'SALI') setattr(ncvar, 'units' , "PSS78") - ncvar=ncOUT.createVariable('PRES_PSAL','f',('nTEMP',)) - ncvar[:]=PresT - ncvar=ncOUT.createVariable('PSAL_QC','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PRES_PSAL','f',('nPSAL',)) + ncvar[:]=PresS + ncvar=ncOUT.createVariable('PSAL_QC','f',('nPSAL',)) ncvar[:]=QcS print("dumping nitrate on " + str(outfile), flush=True) @@ -151,7 +157,7 @@ def dump_nitrate_file(outfile, p, Pres, Value, Qc, metadata,mode='w'): def nitrate_algorithm(p, outfile, metadata, writing_mode): F = p._my_float - Pres, _, _ = p.read('TEMP', read_adjusted=False) + Pres, _, _ = read_temp_psal(p)[3:] if len(Pres)<5: print("few values in Coriolis TEMP in " + str(F.filename), flush=True) return diff --git a/src/bitsea/Float/superfloat_oxygen.py b/src/bitsea/Float/superfloat_oxygen.py index d114a488..61142bde 100644 --- a/src/bitsea/Float/superfloat_oxygen.py +++ b/src/bitsea/Float/superfloat_oxygen.py @@ -63,7 +63,7 @@ def argument(): import TREND_ANALYSIS import bitsea.basins.OGS as OGS from bitsea.instruments.var_conversions import FLOATVARS -from commons_local import cross_Med_basins, save_report +from commons_local import cross_Med_basins df_clim = pd.read_csv('EMODNET_climatology.csv',index_col=0) df_cstd = pd.read_csv('EMODNET_stdev.csv',index_col=0) @@ -72,8 +72,18 @@ class Metadata(): def __init__(self, filename): self.filename = str(filename) self.status_var = 'n' - self.drift_code = -5 - self.offset = -999 + self.drift_code = np.nan + self.offset = np.nan + self.drift = np.nan + + +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS def convert_oxygen(p,doxypres,doxyprofile): @@ -81,13 +91,53 @@ def convert_oxygen(p,doxypres,doxyprofile): from micromol/Kg to mmol/m3 ''' if doxypres.size == 0: return doxyprofile - Pres, temp, Qc = p.read("TEMP",read_adjusted=False) - Pres, sali, Qc = p.read("PSAL",read_adjusted=False) + PresT, temp, Qc, Pres, sali, QcS = read_temp_psal(p) + if len(temp) != len(sali): + temp = np.interp(Pres, PresT, temp) SA = gsw.SA_from_SP(sali, Pres, p.lon, p.lat) density = gsw.rho(SA, gsw.CT_from_t(SA, temp, Pres), Pres) density_on_zdoxy = np.interp(doxypres,Pres,density) return doxyprofile * density_on_zdoxy/1000. +def save_df_report(out_dir, df_report, filename): + """ + Save or update a CSV report file by appending the first row of a DataFrame, + avoiding duplicate entries. + + Parameters + ---------- + out_dir : str or Path. Directory where the CSV file is stored. + df_report : pandas.DataFrame DataFrame from which the first row is extracted and saved. + filename : str. Name of the CSV file. + + Returns + ------- + None + + Notes + ----- + If the file already exists, the function updates it by adding the new entry + only if it is not already present. + """ + + # prepare directory + out_dir = Path(out_dir) + out_dir.mkdir(parents=True, exist_ok=True) + file_path = out_dir / filename + + # prepare report + df_new = df_report.iloc[[0]].copy() + if file_path.exists(): + df_existing = pd.read_csv(file_path) + df_combined = pd.concat([df_existing, df_new], ignore_index=True) + else: + df_combined = df_new + + # save report + df_combined = df_combined.drop_duplicates() + df_combined.to_csv(file_path, index=False) + #print(f"saved {file_path} with {len(df_combined)} rows.") + def dump_oxygen_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): nP=len(Pres) tmpfile=outfile.with_suffix('.nc.tmp') @@ -99,12 +149,11 @@ def dump_oxygen_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): if mode=='w': # if not existing file, we'll put header, TEMP, PSAL setattr(ncOUT, 'origin' , 'coriolis') setattr(ncOUT, 'file_origin', metadata.filename) - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, PresS, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) - ncOUT.createDimension('nPSAL', len(PresT)) + ncOUT.createDimension('nPSAL', len(PresS)) ncvar=ncOUT.createVariable("REFERENCE_DATE_TIME", 'c', ("DATETIME",)) ncvar[:]=p.time.strftime("%Y%m%d%H%M%S") @@ -125,13 +174,13 @@ def dump_oxygen_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): ncvar=ncOUT.createVariable('TEMP_QC','f',('nTEMP',)) ncvar[:]=QcT - ncvar=ncOUT.createVariable('PSAL','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PSAL','f',('nPSAL',)) ncvar[:]=Sali setattr(ncvar, 'variable' , 'SALI') setattr(ncvar, 'units' , "PSS78") - ncvar=ncOUT.createVariable('PRES_PSAL','f',('nTEMP',)) - ncvar[:]=PresT - ncvar=ncOUT.createVariable('PSAL_QC','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PRES_PSAL','f',('nPSAL',)) + ncvar[:]=PresS + ncvar=ncOUT.createVariable('PSAL_QC','f',('nPSAL',)) ncvar[:]=QcS print("dumping oxygen on " + str(outfile), flush=True) @@ -153,6 +202,7 @@ def dump_oxygen_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): ncvar[:]=Qc setattr(ncvar, 'status_var' , metadata.status_var) setattr(ncvar, 'drift_code' , metadata.drift_code) + setattr(ncvar, 'drift' , metadata.drift) setattr(ncvar, 'offset' , metadata.offset) setattr(ncvar, 'variable' , 'DOXY') setattr(ncvar, 'units' , "mmol/m3") @@ -176,15 +226,29 @@ def read_doxy(pCor): def Depth_interp(Profilelist, HistoricalDataset): """ - HistoricalDataset is a dictionary: keys are p.ID(), values are tuples (Pres,Value,Qc) - - data at 600m and 800m are created by interpolating data using layer 580-620m - and 780-820 m + Interpolate oxygen (DOXY) values at standard depths (600m and 800m) from historical profiles. + + Parameters + ---------- + Profilelist : list of ProfileList of float profiles to process. + HistoricalDataset : dictDictionary where keys are profile IDs and values are tuples (Pres, Value, Qc) containing pressure, measured oxygen values, and quality control flags. + + What it does + ------------ + 1. For each profile in Profilelist, checks if data exists within ±20m around 600m and 800m. + 2. Interpolates oxygen (DOXY) values at exactly 600m and 800m if sufficient data exists. + 3. Creates a DataFrame summarizing profile ID, time, location, depth, and interpolated values. + 4. Sets a condition flag True if at least one valid oxygen value exists at 600m. + + Returns + ------- + df : DataFrame with interpolated oxygen values at 600m and 800m for each profile. + CONDITION : bool= True if at least one valid oxygen value exists at 600m, False otherwise. """ - THRES = 20 + + THRES = 20 # data are interpolated starting from a layer +/- 20m es 580-620m LIST_DEPTH = [600,800] COUNT=0 - VARNAME ,varmod = 'DOXY' , 'O2o' columnlist = ['ID','time','lat','lon','name','Type','VAR', 'Depth'] df = pd.DataFrame(index=np.arange(0,2), columns=columnlist) for DEPTH in LIST_DEPTH: @@ -198,28 +262,92 @@ def Depth_interp(Profilelist, HistoricalDataset): df.Depth[COUNT] = DEPTH else: df.Depth[COUNT] = DEPTH + # interp. at 600/800m if data at layer +-20m df.VAR[COUNT] = np.interp(DEPTH , Pres[IDX], (Profile[IDX])) COUNT+=1 - CONDITION = df[df.Depth==600].VAR.notnull().any() # no data at 600m + # CONDITION is True if at least one valid value of VAR exists at 600 m + CONDITION = df[df.Depth==600].VAR.notnull().any() return(df, CONDITION) -def trend_analysis(p, Profilelist_hist, Dataset): - starttime = p.time - timedelta(days=365*3) - TI = TimeInterval.fromdatetimes(starttime, p.time) - Profilelist = [p for p in Profilelist_hist if TI.contains(p.time)] +def extract_oxy_timeseries_at_depth(p, Profilelist_hist, Dataset): + """ + Extract and prepare oxygen timeseries at 600m and 800m depths for a given profile. + + Parameters + ---------- + p : Profile Current float profile for which the timeseries is prepared. + Profilelist_hist : List of historical profiles from the same float. + Dataset : dict Dictionary with profile IDs as keys and tuples (Pres, Value, Qc) as values. + + What it does + ------------ + 1. Defines a time window of 3 years before the profile time up to 1 hour after. + 2. Filters historical profiles to include only those within the time window. + 3. Interpolates oxygen (DOXY) values at 600m and 800m from historical profiles. + 4. Identifies the corresponding Mediterranean basin of the profile. + 5. Checks whether there is at least one valid oxygen value at 600m for trend analysis. + + Returns + ------- + df : pandas.DataFrame DataFrame containing interpolated oxygen values at 600m and 800m depths. + NAME_BASIN : str Name of the Mediterranean basin where the profile is located. + condition1_to_detrend : bool True if there is at least one valid oxygen value at 600m, False otherwise. + """ + + starttime = p.time - timedelta(days=365*3) + endtime = p.time + timedelta(hours=1) + TI = TimeInterval.fromdatetimes(starttime, endtime) # fixed bug + Profilelist = [p for p in Profilelist_hist if TI.contains(p.time)] df, condition1_to_detrend = Depth_interp(Profilelist, Dataset) ARGO = Rectangle(float(p.lon) , float(p.lon) , float(p.lat) , float(p.lat)) NAME_BASIN , BORDER_BASIN = cross_Med_basins(ARGO) return df, NAME_BASIN, condition1_to_detrend -def get_trend_report(p, df): +def get_trend_report(p, df, Save_Repo_when_rejected): + """ + Compute trend diagnostics of dissolved oxygen time series at 600m and 800m + for a given float profile. + + Parameters + ---------- + p : Profile Current float profile containing metadata (time, cycle, float info). + df : pandas.DataFrame containing interpolated oxygen values (VAR) at different depths + Save_Repo_when_rejected : bool If True, the function returns early afterregistering basic data (eg cycle time..) + + What it does + ------------ + 1. Initializes a report DataFrame to store trend diagnostics for 600m and 800m. + 2. For each depth (600m and 800m): + - Filters the input DataFrame for the selected depth and float. + - Computes time series properties: + * duration (in days) + * minimum and maximum dates + * presence of valid (non-NaN) data + - Applies quality checks (length of time series and NaNs). + - Calls `trend_conditions` to compute preliminary trend metrics + (e.g. Theil-Sen and RANSAC estimators). + 3. If `Save_Repo_when_rejected` is True: + - Returns the partial report without performing full trend analysis. + 4. Otherwise: + - Combines Theil-Sen and RANSAC results. + - Evaluates the consistency of trend signs (`sign_analysis`). + - Assigns a drift classification using `drift_coding`. + - Computes and stores the trend per year. + - Updates metadata with the drift code. + 5. Repeats the process for both depths (600m and 800m). + + Returns + ------- + df_report : pandas.DataFrame DataFrame summarizing trend diagnostics for each depth + tmp : pandas.DataFrame Subset of the input DataFrame corresponding to 600m depth, used for + """ + LIST_DEPTH = [600,800] - COLUMNS = ['WMO','Depth','DURATION','min_date','max_date','Theil-Sen','RANSAC'] - df_report = pd.DataFrame(index=np.arange(0,2), columns=COLUMNS) - df_report['TREND_TIME_SERIES'] ,df_report['TREND_per_YEAR'] ,df_report['DRIFT_CODE'] = np.nan , np.nan , np.nan + COLUMNS = ['WMO','cycle','Depth','DURATION','min_date','max_date','Theil-Sen','RANSAC','TREND_TIME_SERIES', 'TREND_per_YEAR', 'DRIFT_CODE'] + df_report = pd.DataFrame(index=np.arange(0,2), columns=COLUMNS) COUNT=0 TI_3 = TimeInterval.fromdatetimes(p.time, p.time) @@ -227,33 +355,77 @@ def get_trend_report(p, df): for DEPTH in LIST_DEPTH: tmp = df[(df.Depth == DEPTH) & (df.name == wmo)] tmp.index = np.arange(0,len(tmp.index)) - days, min_d , max_d = CORIOLIS_checks.lenght_timeseries(tmp, 'time') - Bool = CORIOLIS_checks.nans_check(tmp, 'VAR') + days, min_d , max_d = CORIOLIS_checks.lenght_timeseries(tmp,'time') # for each depth calculate length timeserie + Bool_nans = CORIOLIS_checks.nans_check(tmp, 'VAR') #True -> the time series has enough valid data (few NaNs) tmp.dropna(inplace=True) - lst = trend_conditions(wmo,days, Bool , DEPTH,min_d, max_d , TI_3, tmp) + cycle = p._my_float.cycle + if Save_Repo_when_rejected: #save report even if no data at depth + lst = trend_conditions(wmo,cycle,days, Bool_nans , np.nan,min_d, max_d , TI_3, tmp) + df_report.iloc[COUNT,:] = pd.Series(lst, df_report.columns) # only cycle and wmo registered + return df_report + + lst = trend_conditions(wmo,cycle,days, Bool_nans , DEPTH,min_d, max_d , TI_3, tmp) df_report.iloc[COUNT,:] = pd.Series(lst, df_report.columns) serv = df_report.loc[df_report.WMO==wmo] A = np.append(np.array( serv['Theil-Sen']), np.array( serv['RANSAC'])) - if serv.DURATION.any()==0: - pass + if (serv.DURATION == 0).any(): # if duration is 0 days at 600 or 800m + pass # skip sign analysis and trend analysis else: Bool = sign_analysis(A) df_report = TREND_ANALYSIS.drift_coding(wmo, Bool, serv, df_report) + #arr = np.array(df_report.TREND_per_YEAR)[0] + #metadata.drift_code = np.round(arr, 3) COUNT+=1 tmp = df[(df.Depth == 600) & (df.name == wmo)] return df_report, tmp def clim_check(p, df_report, NAME_BASIN, tmp): + """ + Compute climatology-based offset and update trend report diagnostics. + + This function compares observed dissolved oxygen values with a basin + climatological reference and computes an offset, optionally correcting + for the estimated trend. It also updates the report DataFrame with + climatology-related metadata. + + Parameters + ---------- + p : Profile Current float profile containing metadata (e.g. WMO, cycle). + df_report : pandas.DataFrame Trend report containing results such as trend per year and time series trend. + NAME_BASIN : str Name of the Mediterranean basin associated with the profile. + tmp : pandas.DataFrame Time series subset at a given depth, containing oxygen values (VAR). + + Returns + ------- + OFFSET : float, Difference between observed (or trend-corrected) value and climatology. + threshold : float, Climatological threshold defined as 2 × standard deviation. + df_report : pandas.DataFrame Updated report including OFFSET, basin, climatology, and threshold. + + Notes + ----- + - If the trend is missing or < 1, the offset is computed directly from + the last observed value. + - Otherwise, the last value is corrected by removing the estimated trend + before computing the offset. + """ VALCLIM = float(df_clim.loc[df_clim.index==NAME_BASIN].iloc[:,0]) - TREND_null = df_report.TREND_TIME_SERIES.isnull().values.any() - if TREND_null: - OFFSET = float(tmp.VAR.iloc[-1]) - VALCLIM + TREND_null = df_report.TREND_per_YEAR.isnull().values.any() # bug fixed + + vals = tmp.loc[tmp.Depth == 600, "VAR"].dropna().tail(3) + if vals.empty: + raise ValueError("clim_check(): no valid 600m values available for OFFSET computation") + + observed = float(vals.iloc[-1]) + if TREND_null or abs(df_report.TREND_per_YEAR[0]) < 1: + OFFSET = observed - VALCLIM else: - Corrrected_val = float(tmp.VAR.iloc[-1]) - float(df_report.TREND_TIME_SERIES.iloc[0]) - OFFSET = Corrrected_val - VALCLIM + corrected_val = observed - float(df_report.TREND_TIME_SERIES.iloc[0]) + OFFSET = corrected_val - VALCLIM STDCLIM = float(df_cstd.loc[df_cstd.index==NAME_BASIN].iloc[:,0]) STDCLIM_2 = 2*STDCLIM + if STDCLIM_2 <10: + STDCLIM_2=10 wmo = p._my_float.wmo df_report.loc[ df_report.WMO == wmo , 'OFFSET'] = OFFSET df_report.loc[ df_report.WMO == wmo , 'basin'] = NAME_BASIN @@ -297,51 +469,88 @@ def doxy_algorithm(p, Profilelist_hist, Dataset, outfile, metadata,writing_mode) * Dataset * dictionary, provided by load_history() ''' Pres, Value, Qc = Dataset[p.ID()] - PresT, _, _ = p.read('TEMP', read_adjusted=False) - if len(PresT)<5: + PresT, _, _, _, _, _ = read_temp_psal(p) + if PresT is None or len(PresT) < 5: print("few values in Coriolis TEMP in " + str(p._my_float.filename), flush=True) return metadata.status_var = p._my_float.status_var('DOXY') - df, NAME_BASIN, condition1_to_detrend = trend_analysis(p, Profilelist_hist, Dataset) + # condition1_to_detrend data at 600 m must be available (condition1_to_detrend = True) + df, NAME_BASIN, condition1_to_detrend = extract_oxy_timeseries_at_depth(p, Profilelist_hist, Dataset) Oxy_Profile = Value if not condition1_to_detrend: DRIFT_CODE = -1 - print("no detrend possible", flush=True) + print("no detrend possible, no data at depths", flush=True) metadata.drift_code = DRIFT_CODE + df_report = get_trend_report(p, df, True) # true it write basic info and exit + df_report['basin']=NAME_BASIN else: - df_report, tmp = get_trend_report(p, df) - if tmp[tmp.Depth==600].VAR.isnull().values.any() : + df_report, tmp = get_trend_report(p, df, False) # false it write basic info and go on in writing drift info + df_report['basin'] = NAME_BASIN + # No drift correction if timeseries is short or last 3 values at 600 m are NaN + if len(tmp)<5 or tmp.loc[tmp.Depth == 600, 'VAR'].tail(3).isnull().all(): #bug fixed DRIFT_CODE = -1 metadata.drift_code = DRIFT_CODE else: OFFSET, threshold, df_report = clim_check(p,df_report, NAME_BASIN, tmp) - if abs(OFFSET) >= threshold: + if abs(OFFSET) >= threshold: # discard profiles wmo = p._my_float.wmo df_report.loc[ (df_report.WMO == wmo) & (df_report.Depth== 600), 'Black_list'] = 'True' - timenum = int(p.time.strftime("%Y%m%d")) - save_report( OUT_META / "Blacklist_wmo.csv", 1,['WMO', 'DATE_DAY' , 'OFFSET' , 'STDCLIM_2'],[int(wmo), timenum, OFFSET , threshold]) + save_df_report(OUT_META , df_report , "Floats_rejected.csv") return - else: - Oxy_Profile = apply_detrend(Pres, Value, df_report) - metadata.drift_code = df_report['DRIFT_CODE'].iloc[0] - metadata.offset = OFFSET - # high freq.csv - + else: + if df_report.DRIFT_CODE[0] ==1: + Oxy_Profile = apply_detrend(Pres, Value, df_report) + metadata.drift_code = df_report['DRIFT_CODE'].iloc[0] + metadata.offset = OFFSET + metadata.drift = df_report.TREND_TIME_SERIES.iloc[0] + else: + metadata.drift_code = df_report['DRIFT_CODE'].iloc[0] + metadata.offset = OFFSET + + + # posso metterlo opzionale + save_df_report(OUT_META , df_report , "Floats_accepted.csv") Path.mkdir(outfile.parent, exist_ok=True) dump_oxygen_file(outfile, p, Pres, Oxy_Profile, Qc, metadata,mode=writing_mode) -def load_history(wmo): +def load_history(wmo, write_report=False): ''' Replicates superfloat dataset without detrend - the previous doxy_algorithm + Adds optional logging to a CSV for rejected profiles + + Input: wmo : str , WMO code of the float to load + write_report : bool, If True, rejected profiles (status 'R', missing data, saturation test fail) are logged to a CSV in OUTDIR. + + Output: Profilelist : List of profiles that passed Data Mode and Saturation checks + Dataset : Dictionary of data keyed by profile ID, values = (Pres, Value, Qc) ''' + + global OUTDIR + df_reject = None + csv_file = OUTDIR / "DataMode_and_Saturation_rejection_doxy.csv" + + if write_report: + if csv_file.exists(): + df_reject = pd.read_csv(csv_file) + else: + df_reject = pd.DataFrame(columns=["wmo","cycle","date","reasons"]) + + def log_reject(wmo_name, cycle,date, reason, NAME_BASIN): + nonlocal df_reject + df_reject = pd.concat( + [df_reject, pd.DataFrame({"wmo": [wmo_name], "cycle": [cycle] , "date": [date],"reasons": [reason], "basin": [NAME_BASIN] })], + ignore_index=True) + df_reject.drop_duplicates(inplace=True) + print("Loading dataset for float", wmo, "...", flush=True) + TI = TimeInterval("1950","2050",'%Y') R = Rectangle(-6,36,30,46) PROFILES_COR =bio_float.FloatSelector('DOXY', TI, R) @@ -349,24 +558,40 @@ def load_history(wmo): Profilelist=[] Dataset={} for p in Profilelist_wmo: - if p._my_float.status_var('DOXY')=='R': continue - Pres, Value, Qc = read_doxy(p) - if Pres is None: continue - if p._my_float.status_var('DOXY')=='D': + + ARGO = Rectangle(float(p.lon) , float(p.lon) , float(p.lat) , float(p.lat)) + NAME_BASIN , BORDER_BASIN = cross_Med_basins(ARGO) + status = p._my_float.status_var('DOXY') + + if status == 'R': + if write_report: + log_reject(p.name(), p._my_float.cycle, p.time, "RT" ,NAME_BASIN ) + continue + + Pres, Value, Qc = read_doxy(p) # convert doxy and filter len pres<5 + if Pres is None: + if write_report: + log_reject(p.name(), p._my_float.cycle, p.time, "PresNone", NAME_BASIN) + continue + + if status == 'D': condition_to_write = True - else: - condition_to_write = oxygen_saturation.oxy_check(Pres,Value,p) + else: #ADJUSTED MODE + condition_to_write = oxygen_saturation.oxy_check(Pres, Value, p) + if not condition_to_write: - print(p._my_float.filename, "Saturation Test not passed", flush=True) + if write_report: + log_reject(p.name(), p._my_float.cycle, p.time, "SaturationTest", NAME_BASIN) continue - Profilelist.append(p) Dataset[p.ID()] = (Pres,Value,Qc) + + if write_report and df_reject is not None: + df_reject.to_csv(csv_file, index=False) print("...done", flush=True) return Profilelist, Dataset - OUTDIR = Path(args.outdir) OUT_META = Path(args.outdiag) input_file=args.update_file @@ -382,11 +607,14 @@ def load_history(wmo): for wmo in wmo_list: print (wmo, flush=True) - - Hist_filtered_Profilelist, Dataset = load_history(wmo) + #if wmo != '6903818': continue + Hist_filtered_Profilelist, Dataset = load_history(wmo,True) #true --> save a repo of discarded prechecks Selected_Profilelist=bio_float.filter_by_wmo(PROFILES_COR, wmo) Profilelist= [p for p in Selected_Profilelist if p in Hist_filtered_Profilelist] for ip, p in enumerate(Profilelist): + #import sys + #if p._my_float.cycle == 250: + # sys.exit('carol') outfile = get_outfile(p,OUTDIR) writing_mode=superfloat_generator.writing_mode(outfile) diff --git a/src/bitsea/Float/superfloat_par.py b/src/bitsea/Float/superfloat_par.py index 5880cb3b..6dfacf0f 100644 --- a/src/bitsea/Float/superfloat_par.py +++ b/src/bitsea/Float/superfloat_par.py @@ -1,5 +1,13 @@ import argparse from bitsea.utilities.argparse_types import existing_dir_path +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of downwelling PAR. @@ -68,12 +76,11 @@ def dump_par_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): if mode=='w': # if not existing file, we'll put header, TEMP, PSAL setattr(ncOUT, 'origin' , 'coriolis') setattr(ncOUT, 'file_origin', metadata.filename) - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, PresS, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) - ncOUT.createDimension('nPSAL', len(PresT)) + ncOUT.createDimension('nPSAL', len(PresS)) ncvar=ncOUT.createVariable("REFERENCE_DATE_TIME", 'c', ("DATETIME",)) ncvar[:]=p.time.strftime("%Y%m%d%H%M%S") @@ -93,13 +100,13 @@ def dump_par_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): ncvar=ncOUT.createVariable('TEMP_QC','f',('nTEMP',)) ncvar[:]=QcT - ncvar=ncOUT.createVariable('PSAL','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PSAL','f',('nPSAL',)) ncvar[:]=Sali setattr(ncvar, 'variable' , 'SALI') setattr(ncvar, 'units' , "PSS78") - ncvar=ncOUT.createVariable('PRES_PSAL','f',('nTEMP',)) - ncvar[:]=PresT - ncvar=ncOUT.createVariable('PSAL_QC','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PRES_PSAL','f',('nPSAL',)) + ncvar[:]=PresS + ncvar=ncOUT.createVariable('PSAL_QC','f',('nPSAL',)) ncvar[:]=QcS print("dumping par on " + str(outfile), flush=True) @@ -134,7 +141,7 @@ def get_outfile(p,outdir): def par_algorithm(pCor, outfile, metadata,writing_mode): - Pres, _, _ = pCor.read('TEMP', read_adjusted=False) + Pres, _, _ = read_temp_psal(pCor)[3:] if len(Pres)<5: print("few values in Coriolis TEMP in " + str(pCor._my_float.filename), flush=True) return diff --git a/src/bitsea/Float/superfloat_ph.py b/src/bitsea/Float/superfloat_ph.py index 298641a9..9ffedeab 100644 --- a/src/bitsea/Float/superfloat_ph.py +++ b/src/bitsea/Float/superfloat_ph.py @@ -1,5 +1,13 @@ import argparse from bitsea.utilities.argparse_types import existing_dir_path +def read_temp_psal(p): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=True) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=True) + if (Pres is None or PresT is None or Temp is None or Sali is None or len(Pres) < 5 or len(PresT) < 5): + PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) + Pres, Sali, QcS = p.read('PSAL', read_adjusted=False) + return PresT, Temp, QcT, Pres, Sali, QcS + def argument(): parser = argparse.ArgumentParser(description = ''' Creates superfloat files of PH_IN_SITU_TOTAL. @@ -66,12 +74,11 @@ def dump_ph_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): if mode=='w': # if not existing file, we'll put header, TEMP, PSAL setattr(ncOUT, 'origin' , 'coriolis') setattr(ncOUT, 'file_origin', metadata.filename) - PresT, Temp, QcT = p.read('TEMP', read_adjusted=False) - PresT, Sali, QcS = p.read('PSAL', read_adjusted=False) + PresT, Temp, QcT, PresS, Sali, QcS = read_temp_psal(p) ncOUT.createDimension("DATETIME",14) ncOUT.createDimension("NPROF", 1) ncOUT.createDimension('nTEMP', len(PresT)) - ncOUT.createDimension('nPSAL', len(PresT)) + ncOUT.createDimension('nPSAL', len(PresS)) ncvar=ncOUT.createVariable("REFERENCE_DATE_TIME", 'c', ("DATETIME",)) ncvar[:]=p.time.strftime("%Y%m%d%H%M%S") @@ -91,13 +98,13 @@ def dump_ph_file(outfile, p, Pres, Value, Qc, metadata, mode='w'): ncvar=ncOUT.createVariable('TEMP_QC','f',('nTEMP',)) ncvar[:]=QcT - ncvar=ncOUT.createVariable('PSAL','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PSAL','f',('nPSAL',)) ncvar[:]=Sali setattr(ncvar, 'variable' , 'SALI') setattr(ncvar, 'units' , "PSS78") - ncvar=ncOUT.createVariable('PRES_PSAL','f',('nTEMP',)) - ncvar[:]=PresT - ncvar=ncOUT.createVariable('PSAL_QC','f',('nTEMP',)) + ncvar=ncOUT.createVariable('PRES_PSAL','f',('nPSAL',)) + ncvar[:]=PresS + ncvar=ncOUT.createVariable('PSAL_QC','f',('nPSAL',)) ncvar[:]=QcS print("dumping ph on " + str(outfile), flush=True) @@ -130,7 +137,7 @@ def get_outfile(p,outdir): return filename def ph_algorithm(pCor, outfile, metadata,writing_mode): - Pres, _, _ = pCor.read('TEMP', read_adjusted=False) + Pres, _, _ = read_temp_psal(pCor)[3:] if len(Pres)<5: print("few values in Coriolis TEMP in " + str(pCor._my_float.filename), flush=True) return diff --git a/src/bitsea/basins/basin.py b/src/bitsea/basins/basin.py index c985f856..5c9325ab 100644 --- a/src/bitsea/basins/basin.py +++ b/src/bitsea/basins/basin.py @@ -1,9 +1,8 @@ # Copyright (c) 2015 eXact Lab srl # Author: Stefano Piani import importlib -from collections.abc import Sequence from inspect import currentframe -from typing import Tuple +from typing import Sequence, Tuple import matplotlib.pyplot as plt import numpy as np diff --git a/src/bitsea/instruments/float_canyonmed.py b/src/bitsea/instruments/float_canyonmed.py new file mode 100644 index 00000000..af60fa60 --- /dev/null +++ b/src/bitsea/instruments/float_canyonmed.py @@ -0,0 +1,505 @@ +import netCDF4 +import numpy as np +import datetime +import os +import matplotlib.pyplot as pl +from bitsea.commons.utils import addsep +from bitsea.instruments.instrument import Instrument, Profile +from scipy.optimize import curve_fit +from bitsea.instruments.var_conversions import CANYONVAR as conversion + +mydtype= np.dtype([ + ('file_name','U200'), + ('lat',np.float32), + ('lon',np.float32), + ('time','U17'), + ('parameters','U200'), + ('prof_flags' ,'U200')] + ) +GSS_DEFAULT_LOC ='/g100_scratch/userexternal/camadio0/ARGOPY_TESTS/Climatologies_Argopy/' +ONLINE_REPO = addsep(os.getenv("ONLINE_REPO",GSS_DEFAULT_LOC)) +FloatIndexer=addsep(ONLINE_REPO) + "CANYON_MED_QC/Float_Index.txt" +INDEX_FILE=np.loadtxt(FloatIndexer,dtype=mydtype, delimiter=",",ndmin=1) + +nFiles=INDEX_FILE.size +for iFile in range(nFiles): + INDEX_FILE['parameters'][iFile] = INDEX_FILE['parameters'][iFile].replace("NITRATE_PPCON","NITRATE").replace("CHLA_PPCON","CHLA").replace("BBP700_PPCON","BBP700") + + +class BioFloatProfile(Profile): + def __init__(self, time, lon, lat, my_float, available_params , profile_flags ,mean=None): + self.time = time + self.lon = lon + self.lat = lat + self._my_float = my_float + if available_params[0]==" ": + self.available_params = available_params[1:] + else: + self.available_params = available_params + self.profile_flags= profile_flags + self.mean = mean + self.has_adjusted=False + + def __eq__(self, other): + if isinstance(other, BioFloatProfile): + if (self.lon == other.lon) & (self.lat == other.lat) & (self.time == other.time): + return self._my_float == other._my_float + else: + return False + else: + return False + + def read(self,var,var_mod=None,sourcedata='best', force_reading=None): + ''' + Reads profile data from file. Wrapper for BioFloat.read() + + Arguments: + * var * string + * read_adjusted * IS NOT USED, but we leave it here because there is a lot of calls using it + Once all bit.sea code will use superfloat instead of lovbio_float, we'll remove it + + Returns 3 numpy arrays: Pres, Profile, Qc + ProfileType codes= "best" take all profiles || "insitu" take insitu only || "ppcon" take ppcon only + ProfileType default=> "best" + ''' + + return self._my_float.read(var,var_mod,sourcedata ,force_reading=None) + + + def read_fitted(self,var, func): + ''' + This is used to calculate the radiometric profile values at the surface. + Use an exponential fit and extrapolate values to the surface. + Returns 3 numpy arrays: Pres, Profile, Qc with values extrapolated to 0 m + ''' + + return self._my_float.read_fitted(var, func, mean=self.mean) + + + def name(self): + '''returns a string, the wmo of the associated BioFloat. + ''' + return self._my_float.wmo + + def ID(self): + return self.name() + "_" + self.time.strftime("%Y%m%d-%H:%M:%S_") + str(self.lon) + "_"+ str(self.lat) + def reference_var(self,var): + ''' + Returns the reference varname, for a given profile object and + a ogstm model varname + For BioFloats p.reference_var('O2o') returns 'DOXY' + ''' + return conversion[var] + +class BioFloat(Instrument): + + default_mean = None + + def __init__(self,lon,lat,time,filename,available_params, profile_flags): + self.lon = lon + self.lat = lat + self.time = time + self.filename = filename + if available_params[0]==" ": + self.available_params = available_params[1:] + else: + self.available_params = available_params + self.profile_flags = profile_flags.replace(" ","") + wmo, cycle = os.path.basename(filename).rsplit("_") + self.wmo = wmo[2:] + self.cycle = int(cycle[:3]) + + def __eq__(self,other): + if isinstance(other, BioFloat): + if (self.filename == other.filename): + return (self.lon == other.lon ) & (self.lat == other.lat) & (self.time == other.time) + else: + return False + else: + return False + + + + def status_profile(self,var): + """return code or extended nomenclature for insitu, ppcon, both + """ + li = list(self.available_params.split(" ")) + lf = self.profile_flags + iParam = li.index(var) + flags_code = lf[iParam] + return (flags_code) + + + def has_insitu(self, var): + if self.status_profile(var) in ['I']: return True + return False + def has_tris(self, var): + if self.status_profile(var) in['T']: return True + return False + def has_double(self, var): + if self.status_profile(var) == 'D': return True + return False + def has_issue(self, var): + if self.status_profile(var) == 'X': return True + return False + + + def read_raw(self , var , flag_code, force_reading=None): + if force_reading: + flag_code = False + + ncIN= netCDF4.Dataset(self.filename, 'r') + vars_available = ncIN.variables.keys() + + if flag_code in ['I', 'B'] or force_reading == 'II': + pres_name = 'PRES_' + var + var_name = var + qc_name = var + '_QC' + + elif flag_code in ['P'] or force_reading == 'PP': + pres_name = 'PRES_' + var + '_PPCON' + var_name = var + '_PPCON' + qc_name = var + '_PPCON_QC' + + elif flag_code in ['C'] or force_reading == 'CC': + pres_name = 'PRES_' + var + '_CANYONMED' + var_name = var + '_CANYONMED' + qc_name = var + '_CANYONMED_QC' + + else: + raise ValueError(f"Invalid flag_code ({flag_code}) or force_reading ({force_reading})") + + for name in [pres_name, var_name, qc_name]: + if name not in vars_available: + raise KeyError(f"Variable '{name}' not found in {self.filename}") + + Pres = np.array(ncIN.variables[pres_name]) + Profile = np.array(ncIN.variables[var_name]) + Qc = np.array(ncIN.variables[qc_name]) + + ncIN.close() + return Pres, Profile, Qc + + + def read(self,var,var_mod=None,sourcedata='best',force_reading=None): + ''' + With sourcedata='best' gives priority to insitu + ''' + assert sourcedata in ["best","insitu","ppcon",'canyon_med'] + + if sourcedata=='best': + flag_code = self.status_profile(var) # get flag insitu ppcon both none + Pres, Profile, Qc = self.read_raw(var, flag_code,force_reading=None) + + elif sourcedata=="insitu": + assert self.has_insitu(var) + Pres, Profile, Qc = self.read_raw(var, 'I',force_reading=None) + if (Qc <0).any() : + raise ValueError("Incorrect QC flag for canyon_med") + + elif sourcedata=="ppcon": + assert self.has_ppcon(var) + Pres, Profile, Qc = self.read_raw(var, 'P',force_reading='PP') + if (Qc >0).any() : + raise ValueError("Incorrect QC flag for canyon_med") + + elif sourcedata=="canyon_med": + print(self.filename) + flag_code = self.status_profile(var) + Pres, Profile, Qc = self.read_raw(var, 'X',force_reading='CC') + if (Qc >0).any() : + raise ValueError("Incorrect QC flag for canyon_med") + + else: + raise ValueError(f"Invalid sourcedata option: {sourcedata}") + + if var_mod is None: + return Pres, Profile, Qc + + #assert var_mod in ['P_c', 'POC'] + ii=(Pres >= 180) & (Pres <= 200) + if (var_mod=='P_c'): + bbp470 = Profile * ( 470.0/ 700)**(-0.78)# [m-1] + Profile = 12128 * bbp470 + 0.59 # Conversion by Bellacicco 201? + if ii.sum() > 0 : + shift=Profile[ii].mean() + print( var_mod +" : adding a shift of " + str(shift)) + Profile = Profile - shift + ii=Profile<=0 + Profile[ii] = 0.0 + + if (var_mod == "POC"): + POC = Profile * 52779.37 - 3.57 # Bellacicco 2019 + if ii.sum() > 0 : + shift=POC[ii].mean() + print( "POC: adding a shift of " + str(shift)) + Profile = POC - shift + ii=Profile<=0 + Profile[ii] = 0.0 + + + return Pres, Profile, Qc + + + def read_fitted(self, var, func, mean=None): + ''' + This is used to calculate the radiometric profile values at the surface. + Use an exponential fit and extrapolate values to the surface. + Returns the irradiance values at the surface. + + Reads profile data from file, applies a rarefaction and optionally a filter to the data + + Returns the measurement's surface value + + ''' + + raw_pres, raw_prof = self.read_raw(var) + pres,index = np.unique(raw_pres,return_index=True) + prof = raw_prof[index] + + if pres[0] < 1.5 and pres[3] < 10.: + + popt,_ = curve_fit(func, pres, prof) + pres_new = np.insert(pres, 0, 0.) if not pres[0] == 0. else pres # Add z=0 m + prof_new = func(pres_new, *popt) + + else: + pres_new = pres + prof_new = prof + + qc = np.ones_like(pres_new)*2 + + if pres_new.size ==0: + return pres_new, prof_new, qc + + if mean == None: + if BioFloat.default_mean != None: + return pres_new, BioFloat.default_mean.compute(prof_new, pres_new), qc + else: + return pres_new, prof_new, qc + else: + return pres_new, mean.compute(prof_new, pres_new), mean.compute(qc,pres_new) + + + def origin(self,var): + ''' + Arguments: + * var * string + Returns: + a tuple of two strings: + * origin * string, "lov" or "Coriolis" + * file_origin * string , path of the source file + ''' + class Info(): + def __init__(self): + self.file_orig = "None" + self.correction = "None" + self.shift = "None" + self.status_var = 'n' + def __repr__(self): + return "status_var = %s" %(self.status_var) + + info = Info() + + ncIN=netCDF4.Dataset(self.filename,'r') + info.status_var=ncIN.variables[var].status_var + info.file_orig=ncIN.file_origin + + if var=='NITRATE': + info.shift=ncIN.variables[var].bottomshift + info.correction=ncIN.variables[var].correction + ncIN.close() + + return info + + + def basicplot(self,Pres,profile): + pl.figure() + pl.plot(profile,Pres) + pl.gca().invert_yaxis() + pl.show(block=False) + + def plot(self,Pres,profile,fig=None,ax=None,**kwargs): + ''' + Args: + - * Pres * : array of pressure + - * profile * : array of profile values + - * fig * : a reference to a Figure object, if None a new Figure will be created. + - * ax * : a reference to an Axes object, if None a new axis will be created. + + Returns : + A figure and an axis object + + Examples: + fig, ax = f.plot(pres,profile) + fig, ax = f.plot(pres,profile,fig,ax) + fig, ax = f.plot(pres,profile,fig,ax, linestyle = 'dashed', linewidth = 2, color='green') + + ''' + if (fig is None) or (ax is None): + fig , ax = pl.subplots() + ax.plot(profile,Pres, **kwargs) + if not ax.yaxis_inverted(): ax.invert_yaxis() + return fig,ax + + def profiles(self, var, mean=None): + return [BioFloatProfile(var, self.time, self.lon, self.lat, self, mean)] + + @staticmethod + def from_file(filename): + ''' + Returns the single Bio_Float instance corresponding to filename + ''' + + nFiles=INDEX_FILE.size + for iFile in range(nFiles): + timestr = INDEX_FILE['time'][iFile] + lon = INDEX_FILE['lon' ][iFile] + lat = INDEX_FILE['lat' ][iFile] + thefilename = INDEX_FILE['file_name'][iFile] + available_params = INDEX_FILE['parameters'][iFile] + profile_flags = INDEX_FILE['prof_flags'][iFile] + float_time = datetime.datetime.strptime(timestr,'%Y%m%d-%H:%M:%S') + if filename.endswith(thefilename): + return BioFloat(lon,lat,float_time,filename,available_params, profile_flags) + return None + + +def FloatSelector(var, T, region): + ''' + Arguments: + * var * is a string indicating variable, + if var is None, no selection is done about variable + * T * is a TimeInterval instance or a timerequestors.Clim_season instance + or whatever object having a contains() method working as the TimeInteval one does. + * region * is an instance of Region or its derived (Polygon, Basin, ...) + + Returns: + a list of BioFloatProfile objects. + Caveats: + In order to work on dataset different from the cineca DRES archive + /gss/gss_work/DRES_OGS_BiGe/Observations/TIME_RAW_DATA/ONLINE/ + remember to define the environment variable ONLINE_REPO + export ONLINE_REPO=/some/path/with/ COPERNICUS/ FLOAT_BIO/ FLOAT_LOVBIO/ SAT/ SUPERFLOAT/ + ''' + + nFiles=INDEX_FILE.size + selected = [] + for iFile in range(nFiles): + timestr = INDEX_FILE['time'][iFile] + lon = INDEX_FILE['lon' ][iFile] + lat = INDEX_FILE['lat' ][iFile] + filename = INDEX_FILE['file_name'][iFile] + available_params = INDEX_FILE['parameters'][iFile] + profile_flags = INDEX_FILE['prof_flags'][iFile] + float_time = datetime.datetime.strptime(timestr,'%Y%m%d-%H:%M:%S') + filename = ONLINE_REPO + 'CANYON_MED_QC/' + filename + + if var is None : + VarCondition = True + else: + VarCondition = var in available_params + + if VarCondition: + if T.contains(float_time) and region.is_inside(lon, lat): + thefloat = BioFloat(lon,lat,float_time,filename,available_params, profile_flags) + selected.append(BioFloatProfile(float_time,lon,lat, thefloat,available_params,profile_flags,None)) + + return selected + +def get_wmo_list(Profilelist): + ''' + Argument: + * Profilelist * list of Profile objects + + Returns: + a list of wmo strings + ''' + wmo_set=set() + for p in Profilelist: + wmo_set.add(p._my_float.wmo) + return list(wmo_set) + +def filter_by_wmo(Profilelist,wmo): + ''' + + Subsetter, filtering by wmo + + Arguments: + * Profilelist * list of Profile objects + * wmo * string + + Returns: + a list of Profile Objects + ''' + + return [p for p in Profilelist if p._my_float.wmo == wmo] +def remove_bad_sensors(Profilelist,var): + ''' + + Subsetter, filtering out bad sensors for that var. + At the moment that method does not do anything, because bad sensors have been already removed + by superfloat generators. + + In case of new expert evaluation, if we decide to remove another bad sensor, + we can add here the part of code to have quickly the desired result. + Then, the procedure to avoid bad sensors in superfloat dataset is + - adding bad float in superfloat_${var}.py remove_bad_sensor() + - remove files from SUPERFLOAT DATASET + - launch dump_index.py without Float_Indexer.txt input + + + Arguments: + * Profilelist * list of Profile objects + * var * string + + Returns: + a list of Profile Objects + ''' + + return Profilelist + + +if __name__ == '__main__': + filename="/gss/gss_work/DRES_OGS_BiGe/Observations/TIME_RAW_DATA/ONLINE_V10C/CANYON_MED_QC/5907088/SR5907088_016.nc" + F=BioFloat.from_file(filename) + print("profile_flags=",F.profile_flags) + print("status_profile=", F.status_profile('NITRATE')) + print("has_insitu", F.has_insitu('NITRATE')) + F.read('DOXY') + F.read('DOXY',sourcedata='insitu') + Pres_P,Value_P, _ =F.read('CHLA',sourcedata='ppcon') + Pres_I,Value_I, Qc=F.read('CHLA',sourcedata='insitu') + Pres,Value, Qc=F.read('CHLA',sourcedata='best') + assert (Value == Value_I).all() + + + from bitsea.basins.region import Rectangle + from bitsea.commons.time_interval import TimeInterval + + var = 'NITRATE' + TI = TimeInterval('20240201','20240301','%Y%m%d') + R = Rectangle(-6,36,30,46) + + PROFILE_LIST=FloatSelector(var, TI, R) + for p in PROFILE_LIST: + print(p._my_float.has_insitu('NITRATE')) + Pres, Values, Qc = p.read('NITRATE') + + import sys + sys.exit() + + Pres,V, Qc = F.read(var) + import sys + sys.exit() + +# print len(PROFILE_LIST) + + for p in PROFILE_LIST[100:200]: + Pres,V, Qc = p.read(var) +# if Pres.min()>0: +# print Pres.min() + + wmo_list= get_wmo_list(PROFILE_LIST) + for wmo in wmo_list: + sublist = filter_by_wmo(PROFILE_LIST, wmo) diff --git a/src/bitsea/instruments/var_conversions.py b/src/bitsea/instruments/var_conversions.py index d205e83a..43ec4267 100644 --- a/src/bitsea/instruments/var_conversions.py +++ b/src/bitsea/instruments/var_conversions.py @@ -10,7 +10,16 @@ 'P_c':'BBP700', 'BBP532':'BBP532', 'pH' : 'PH_IN_SITU_TOTAL' } - + +CANYONVAR={'O2o':'DOXY', \ + 'P_l':'CHLA', \ + 'N3n':'NITRATE', \ + 'N1p' :'PO4', \ + 'DIC' :'DIC', \ + 'N5s' :'SiOH4', \ + 'ALK' :'AT', \ + 'pH' :'PH_IN_SITU_TOTAL' } + #DOXY:units = "micromole/kg" ; #CHLA:units = "mg/m3" ; #NITRATE:units = "micromole/kg" ; diff --git a/src/bitsea/static/DatasetExtractor.py b/src/bitsea/static/DatasetExtractor.py index e0d1103d..1675fe68 100644 --- a/src/bitsea/static/DatasetExtractor.py +++ b/src/bitsea/static/DatasetExtractor.py @@ -45,6 +45,7 @@ def profileGenerator(self,TIME,lon,lat,values,depth,dataset): lon = uniques[i,1].astype(np.float32) lat = uniques[i,2].astype(np.float32) inthisprofile = ib==i + inthisprofile= inthisprofile.ravel() Values=values[inthisprofile] Depth = depth[inthisprofile]