arts calibration with new dropsonde data, attributes finalized

Author: JakobDeutloff

scripts/arts_calibration/analyse_bt_diffs.py

@@ -2,7 +2,8 @@
 import pandas as pd
 import matplotlib.pyplot as plt
 import xarray as xr
-from src import readwrite_functions as rwfuncs
+import yaml
+import numpy as np
 
 # %% define frequencies
 freq_k = [22.24, 23.04, 23.84, 25.44, 26.24, 27.84, 31.40]
@@ -43,35 +44,14 @@
 
 
 # %% Read csv BTs
-dates = [
-    "20240811",
-    "20240813",
-    "20240816",
-    "20240818",
-    "20240821",
-    "20240822",
-    "20240825",
-    "20240827",
-    "20240829",
-    "20240831",
-    "20240903",
-    "20240906",
-    "20240907",
-    "20240909",
-    "20240912",
-    "20240914",
-    "20240916",
-    "20240919",
-    "20240921",
-    "20240923",
-    "20240924",
-    "20240926",
-    "20240928",
+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 dates:
+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)
     )
@@ -80,18 +60,25 @@
     )
 
 # load dropsonde data
-configfile = "config_ipns.yaml"
-cfg = rwfuncs.extract_config_params(configfile)
-ds_dropsonde = xr.open_dataset(cfg["path_dropsondes"], engine="zarr")
+with open("process_config.yaml") as f:
+    cfg = yaml.safe_load(f)
+
+ds_dropsonde = xr.open_dataset(
+    "ipns://latest.orcestra-campaign.org/products/HALO/dropsondes/Level_3/PERCUSION_Level_3.zarr",
+    engine="zarr",
+)
+ds_dropsonde = ds_dropsonde.assign_coords(sonde=np.arange(ds_dropsonde.sonde.size))
 
-# %% restructure data
+# 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
+launch_time = ds_dropsonde.sel(
+    sonde=[int(float(x)) for x in TB_arts.columns.values]
+).sonde_time.values
 TB_arts.columns = launch_time
 TB_hamp.columns = launch_time
-TB_arts = TB_arts.T
-TB_hamp = TB_hamp.T
+TB_arts = TB_arts.T.dropna()
+TB_hamp = TB_hamp.T.dropna()
 
 # %% calculate statistics for each flight
 diffs = TB_arts - TB_hamp

scripts/arts_calibration/arts_bt_calculation.py

@@ -1,14 +1,13 @@
 # %%
 import os
 import xarray as xr
-import src.readwrite_functions as rwfuncs
-from src.post_processed_hamp_data import PostProcessedHAMPData
 import yaml
 import pandas as pd
 import typhon
 import pyarts
 from tqdm import tqdm
 import FluxSimulator as fsm
+import numpy as np
 from src.arts_functions import (
     Hamp_channels,
     basic_setup,
@@ -20,9 +19,10 @@
     is_complete,
 )
 
+# from src.plots_functions import plot_dropsonde, plot_TB_comparison
+
 
 # %% get ARTS data
-print("Download ARTS data")
 pyarts.cat.download.retrieve(verbose=True)
 
 # %% setup sensor
@@ -34,9 +34,20 @@
 # %% setup workspace
 ws = basic_setup([], sensor_description=sensor_description)
 
+# %% load data
+configfile = "process_config.yaml"
+with open(configfile, "r") as file:
+    cfg = yaml.safe_load(file)
+
+ds_dropsonde = xr.open_dataset(
+    "ipns://latest.orcestra-campaign.org/products/HALO/dropsondes/Level_3/PERCUSION_Level_3.zarr",
+    engine="zarr",
+)
+ds_dropsonde = ds_dropsonde.assign_coords(sonde=np.arange(ds_dropsonde.sonde.size))
+
 
 # %% define function
-def calc_arts_bts(date, flightletter):
+def calc_arts_bts(date, flightletter, ds_dropsonde, cfg):
     """
     Calculates brightness temperatures for the radiometer frequencies with
     ARTS based on the dropsonde profiles for the flight on date.
@@ -50,84 +61,64 @@ def calc_arts_bts(date, flightletter):
     None. Data is saved in arts_comparison folder.
     """
 
-    print("Read Config")
-    configfile = "config_ipns.yaml"
-    with open(configfile, "r") as file:
-        config_yml = yaml.safe_load(file)
-    config_yml["date"] = date
-    config_yml["flightletter"] = flightletter
-    with open(configfile, "w") as file:
-        yaml.dump(config_yml, file)
-    cfg = rwfuncs.extract_config_params(configfile)
-
-    print("Load Bahamas Data")
-    ds_bahamas = xr.open_dataset(cfg["path_position_attitude"], engine="zarr")
-
-    print("Load Dropsonde Data")
-    ds_dropsonde = xr.open_dataset(cfg["path_dropsondes"], engine="zarr")
-    ds_dropsonde = ds_dropsonde.where(
-        (ds_dropsonde["interp_time"] > pd.to_datetime(cfg["date"]))
-        & (
-            ds_dropsonde["interp_time"]
-            < pd.to_datetime(cfg["date"]) + pd.DateOffset(hour=23)
-        )
-    ).dropna(dim="sonde_id", how="all")
-
-    print("Load HAMP Data")
-    hampdata = PostProcessedHAMPData(
-        xr.open_dataset(cfg["path_radar"], engine="zarr"),
-        xr.open_dataset(cfg["path_radiometers"], engine="zarr"),
-        xr.open_dataset(cfg["path_iwv"], engine="zarr"),
+    flightname = f"HALO-{date}{flightletter}"
+    ds_radar = xr.open_dataset(
+        f"{cfg['save_dir']}/radar/{flightname}_radar.zarr", engine="zarr"
     )
+    ds_radiometers = xr.open_dataset(
+        f"{cfg['save_dir']}/radiometer/{flightname}_radio.zarr", engine="zarr"
+    )
+    ds_dropsonde = ds_dropsonde.where(
+        (ds_dropsonde["interp_time"] > pd.to_datetime(date))
+        & (ds_dropsonde["interp_time"] < pd.to_datetime(date) + pd.DateOffset(hour=23))
+    ).dropna(dim="sonde", how="all")
 
     print("Calculate Cloud Mask")
     ds_dropsonde = ds_dropsonde.assign(
         radar_cloud_flag=(
-            hampdata.radar.sel(time=ds_dropsonde.launch_time, method="nearest").sel(
+            ds_radar.sel(time=ds_dropsonde.sonde_time, method="nearest").sel(
                 height=slice(200, None)
             )["dBZe"]
             > -30
         ).max("height")
         * 1
     )
     cloud_free_idxs = (
-        ds_dropsonde["sonde_id"]
-        .where(ds_dropsonde["radar_cloud_flag"] == 0, drop=True)
+        ds_dropsonde["sonde"]
+        .where((ds_dropsonde["radar_cloud_flag"] == 0).compute(), drop=True)
         .values
     )
 
-    freqs_hamp = hampdata.radiometers.frequency.values
+    freqs_hamp = ds_radiometers.frequency.values
     TBs_arts = pd.DataFrame(index=freqs / 1e9, columns=cloud_free_idxs)
     TBs_hamp = pd.DataFrame(index=freqs_hamp, columns=cloud_free_idxs)
 
     print("Setup Folders")
-    if not os.path.exists(f"Data/arts_calibration/{cfg['flightname']}"):
-        os.makedirs(f"Data/arts_calibration/{cfg['flightname']}")
-    if not os.path.exists(f"Data/arts_calibration/{cfg['flightname']}/plots"):
-        os.makedirs(f"Data/arts_calibration/{cfg['flightname']}/plots")
-
-    print(f"Running {cloud_free_idxs.size} dropsondes for {cfg['flightname']}")
-    for sonde_id in tqdm(cloud_free_idxs):
-        ds_dropsonde_loc, hampdata_loc, height, drop_time = get_profiles(
-            sonde_id, ds_dropsonde, hampdata
+    if not os.path.exists(f"Data/arts_calibration/{flightname}"):
+        os.makedirs(f"Data/arts_calibration/{flightname}")
+    if not os.path.exists(f"Data/arts_calibration/{flightname}/plots"):
+        os.makedirs(f"Data/arts_calibration/{flightname}/plots")
+
+    print(f"Running {cloud_free_idxs.size} dropsondes for {date}")
+    for sonde in tqdm(cloud_free_idxs):
+        ds_dropsonde_loc, radiometers_loc, height, drop_time = get_profiles(
+            sonde, ds_dropsonde, ds_radiometers
         )
 
-        if not is_complete(ds_dropsonde_loc, hampdata_loc, drop_time, height, sonde_id):
+        if not is_complete(ds_dropsonde_loc, radiometers_loc, drop_time, height, sonde):
             continue
 
         surface_temp = get_surface_temperature(ds_dropsonde_loc)
         surface_ws = get_surface_windspeed(ds_dropsonde_loc)
 
-        ds_dropsonde_extrap = extrapolate_dropsonde(
-            ds_dropsonde_loc, height, ds_bahamas
-        )
+        ds_dropsonde_extrap = extrapolate_dropsonde(ds_dropsonde_loc, height)
 
         # plot_dropsonde(ds_dropsonde_extrap, ds_dropsonde_loc)
 
         profile_grd = fsm.generate_gridded_field_from_profiles(
             pyarts.arts.Vector(ds_dropsonde_extrap["p"].values),
             ds_dropsonde_extrap["ta"].values,
-            ds_dropsonde_extrap["alt"].values,
+            ds_dropsonde_extrap["altitude"].values,
             gases={
                 "H2O": typhon.physics.specific_humidity2vmr(
                     ds_dropsonde_extrap["q"].values
@@ -145,24 +136,23 @@ def calc_arts_bts(date, flightletter):
             height,
         )
 
-        TBs_arts[sonde_id] = pd.DataFrame(data=BTs, index=freqs / 1e9)
-        TBs_hamp[sonde_id] = hampdata_loc["TBs"]
+        TBs_arts[sonde] = pd.DataFrame(data=BTs, index=freqs / 1e9)
+        TBs_hamp[sonde] = radiometers_loc["TBs"]
 
-        # plot_TB_comparison(TBs_arts[sonde_id], TBs_hamp[sonde_id], sonde_id)
+        # plot_TB_comparison(TBs_arts[sonde], TBs_hamp[sonde], sonde)
 
-    TBs_arts.to_csv(f"Data/arts_calibration/{cfg['flightname']}/TBs_arts.csv")
-    TBs_hamp.to_csv(f"Data/arts_calibration/{cfg['flightname']}/TBs_hamp.csv")
+    TBs_arts.to_csv(f"Data/arts_calibration/{flightname}/TBs_arts.csv")
+    TBs_hamp.to_csv(f"Data/arts_calibration/{flightname}/TBs_hamp.csv")
 
 
 # %% loop over flights
-
 flights = pd.read_csv("flights.csv", index_col=0)
 flights_processed = flights[
     (flights["location"] == "sal") | (flights["location"] == "barbados")
 ]
 
 for date in flights_processed.index:
-    calc_arts_bts(str(date), flights_processed.loc[date]["flightletter"])
-
-
+    calc_arts_bts(
+        str(date), flights_processed.loc[date]["flightletter"], ds_dropsonde, cfg
+    )
 # %%