regression analysis of arts calibration

Author: JakobDeutloff

config_ipns.yaml

@@ -1,4 +1,4 @@
-date: '20240929'
+date: '20240928'
 flightletter: a
 flightname: HALO-{date}{flightletter}
 iwv: ipns://latest.orcestra-campaign.org/products/HALO/iwv/HALO-{date}{flightletter}.zarr

flights.csv

@@ -0,0 +1,33 @@
+,flightletter,location
+20240809,a,transfer
+20240811,a,sal
+20240813,a,sal
+20240816,a,sal
+20240818,a,sal
+20240821,a,sal
+20240822,a,sal
+20240825,a,sal
+20240827,a,sal
+20240829,a,sal
+20240831,a,sal
+20240903,a,sal
+20240906,a,transfer
+20240907,a,barbados
+20240909,a,barbados
+20240912,a,barbados
+20240914,a,barbados
+20240916,a,barbados
+20240919,a,barbados
+20240921,a,barbados
+20240923,a,barbados
+20240924,a,barbados
+20240926,a,barbados
+20240928,a,barbados
+20240929,a,transfer
+20241105,a,op
+20241107,a,op
+20241110,a,op
+20241112,b,op
+20241114,b,op
+20241116,a,op
+20241119,a,op

scripts/arts_calibration/arts_bt_calculation.py

@@ -154,41 +154,15 @@ def calc_arts_bts(date, flightletter):
     TBs_hamp.to_csv(f"Data/arts_calibration/{cfg['flightname']}/TBs_hamp.csv")
 
 
-# %% call function
-dates = [
-    "20240811",
-    "20240813",
-    "20240816",
-    "20240818",
-    "20240821",
-    "20240822",
-    "20240825",
-    "20240827",
-    "20240829",
-    "20240831",
-    "20240903",
-    "20240906",
-    "20240907",
-    "20240909",
-    "20240912",
-    "20240914",
-    "20240916",
-    "20240919",
-    "20240921",
-    "20240923",
-    "20240924",
-    "20240926",
-    "20240928",
-    "20240929",
-]
+# %% loop over flights
 
-flightletters = {}
-for date in dates:
-    flightletters[date] = "a"
+flights = pd.read_csv("flights.csv", index_col=0)
+flights_processed = flights[
+    (flights["location"] == "sal") | (flights["location"] == "barbados")
+]
 
-# %%
-for date in dates:
-    calc_arts_bts(date, flightletters[date])
+for date in flights_processed.index:
+    calc_arts_bts(str(date), flights_processed.loc[date]["flightletter"])
 
 
 # %%

scripts/arts_calibration/regressions.py

@@ -0,0 +1,86 @@
+# %%
+import numpy as np
+from scipy.stats import linregress
+import pandas as pd
+import xarray as xr
+from src import readwrite_functions as rwfuncs
+from src.plots_functions import plot_regression
+
+# %% Read csv BTs
+
+flights = pd.read_csv("flights.csv", index_col=0)
+flights_processed = flights[
+    (flights["location"] == "sal") | (flights["location"] == "barbados")
+]
+
+TB_arts_list = []
+TB_hamp_list = []
+for date in flights_processed.index:
+    TB_arts_list.append(
+        pd.read_csv(f"Data/arts_calibration/HALO-{date}a/TBs_arts.csv", index_col=0)
+    )
+    TB_hamp_list.append(
+        pd.read_csv(f"Data/arts_calibration/HALO-{date}a/TBs_hamp.csv", index_col=0)
+    )
+
+# load dropsonde data
+configfile = "config_ipns.yaml"
+cfg = rwfuncs.extract_config_params(configfile)
+ds_dropsonde = xr.open_dataset(cfg["path_dropsondes"], engine="zarr")
+
+# restructure data
+TB_arts = pd.concat(TB_arts_list, axis=1)
+TB_hamp = pd.concat(TB_hamp_list, axis=1)
+launch_time = ds_dropsonde.sel(sonde_id=TB_arts.columns).launch_time.values
+TB_arts.columns = launch_time
+TB_hamp.columns = launch_time
+TB_arts = TB_arts.T.dropna()
+TB_hamp = TB_hamp.T.dropna()
+
+# %% drop extreme outliers
+TB_arts_std = TB_arts.std("index")
+TB_hamp_std = TB_hamp.std("index")
+TB_arts_med = TB_arts.median("index")
+TB_hamp_med = TB_hamp.median("index")
+TB_arts = TB_arts.where((TB_arts - TB_arts_med).abs() < 3 * TB_arts_std)
+TB_hamp = TB_hamp.where((TB_hamp - TB_hamp_med).abs() < 3 * TB_hamp_std)
+
+# %% count NaN values
+mask_arts = ((TB_arts - TB_arts_med).abs() < 2 * TB_arts_std).mean("index")
+mask_hamp = ((TB_hamp - TB_hamp_med).abs() < 2 * TB_hamp_std).mean("index")
+mask_arts
+
+# %%
+mask_hamp
+
+
+# %%
+multi_index = pd.MultiIndex.from_product(
+    [TB_arts.columns, ["slope", "intercept"]], names=["frequency", "parameter"]
+)
+regression_coeffs = pd.DataFrame(
+    index=np.unique(TB_arts.index.date), columns=multi_index
+)
+
+for date in flights_processed.index:
+    date = pd.to_datetime(date, format="%Y%m%d").date()
+    TB_arts_day = TB_arts.loc[TB_arts.index.date == date]
+    TB_hamp_day = TB_hamp.loc[TB_hamp.index.date == date]
+
+    for f in TB_arts_day.columns:
+        # Flatten the arrays and filter out NaN values
+        x = TB_hamp_day[f].values.flatten()
+        y = TB_arts_day[f].values.flatten()
+        mask = ~np.isnan(x) & ~np.isnan(y)
+        x = x[mask]
+        y = y[mask]
+
+        if len(x) > 0 and len(y) > 0:  # Ensure there are values left after filtering
+            regression = linregress(x, y)
+            regression_coeffs.loc[date, (f, "slope")] = regression.slope
+            regression_coeffs.loc[date, (f, "intercept")] = regression.intercept
+
+# %% plot regressions
+plot_regression(regression_coeffs, TB_arts, TB_hamp, date)
+
+# %%