ARTS calibration with new setup running

Author: JakobDeutloff

config_ipns.yaml

@@ -1,8 +1,9 @@
-date: '20240811'
+date: '20240929'
 flightletter: a
 flightname: HALO-{date}{flightletter}
 iwv: ipns://latest.orcestra-campaign.org/products/HALO/iwv/HALO-{date}{flightletter}.zarr
 path_dropsondes: ipns://latest.orcestra-campaign.org/products/HALO/dropsondes/Level_3/PERCUSION_Level_3.zarr
 path_saveplots: ''
+position_attitude: ipns://latest.orcestra-campaign.org/products/HALO/position_attitude/HALO-{date}{flightletter}.zarr
 radar: ipns://latest.orcestra-campaign.org/products/HALO/radar/moments/HALO-{date}{flightletter}.zarr
 radiometer: ipns://latest.orcestra-campaign.org/products/HALO/radiometer/HALO-{date}{flightletter}.zarr

scripts/arts_calibration/analyse_bt_diffs.py

@@ -41,7 +41,6 @@
     "20240928",
 ]
 
-next_flight_183 = ["202408"]
 
 # %% Read csv BTs
 dates = [

scripts/arts_calibration/arts_bt_calculation.py

@@ -1,6 +1,14 @@
 # %%
 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
 from src.arts_functions import (
     Hamp_channels,
     basic_setup,
@@ -9,19 +17,8 @@
     get_surface_temperature,
     get_surface_windspeed,
     forward_model,
-    check_data,
+    is_complete,
 )
-from src.ipfs_helpers import read_nc
-from orcestra.postprocess.level0 import bahamas
-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
-from src.plots_functions import plot_dropsonde
 
 
 # %% get ARTS data
@@ -39,7 +36,7 @@
 
 
 # %% define function
-def calc_arts_bts(date):
+def calc_arts_bts(date, flightletter):
     """
     Calculates brightness temperatures for the radiometer frequencies with
     ARTS based on the dropsonde profiles for the flight on date.
@@ -54,27 +51,18 @@ def calc_arts_bts(date):
     """
 
     print("Read Config")
-    # change the date in config.yaml to date
     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)
-    # read config
     cfg = rwfuncs.extract_config_params(configfile)
 
-    # load bahamas data from ipfs
     print("Load Bahamas Data")
-    ds_bahamas = (
-        read_nc(
-            f"ipns://latest.orcestra-campaign.org/raw/HALO/bahamas/{cfg['flightname']}/QL_*.nc"
-        )
-        .pipe(bahamas)
-        .interpolate_na("time")
-    )
+    ds_bahamas = xr.open_dataset(cfg["path_position_attitude"], engine="zarr")
 
-    # read dropsonde data
     print("Load Dropsonde Data")
     ds_dropsonde = xr.open_dataset(cfg["path_dropsondes"], engine="zarr")
     ds_dropsonde = ds_dropsonde.where(
@@ -85,15 +73,13 @@ def calc_arts_bts(date):
         )
     ).dropna(dim="sonde_id", how="all")
 
-    # read HAMP post-processed data
     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"),
     )
 
-    # cloud mask from radar
     print("Calculate Cloud Mask")
     ds_dropsonde = ds_dropsonde.assign(
         radar_cloud_flag=(
@@ -110,48 +96,38 @@ def calc_arts_bts(date):
         .values
     )
 
-    # initialize result arrays
     freqs_hamp = hampdata.radiometers.frequency.values
-    TBs_arts = pd.DataFrame(index=freqs_hamp, columns=cloud_free_idxs)
-    TBs_hamp = TBs_arts.copy()
+    TBs_arts = pd.DataFrame(index=freqs / 1e9, columns=cloud_free_idxs)
+    TBs_hamp = pd.DataFrame(index=freqs_hamp, columns=cloud_free_idxs)
 
-    # setup folders
     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")
 
-    # loop over cloud free sondes
     print(f"Running {cloud_free_idxs.size} dropsondes for {cfg['flightname']}")
     for sonde_id in tqdm(cloud_free_idxs):
-        # get profiles
         ds_dropsonde_loc, hampdata_loc, height, drop_time = get_profiles(
             sonde_id, ds_dropsonde, hampdata
         )
 
-        # check if dropsonde is broken (contains only nan values)
-        if not check_data(
-            ds_dropsonde_loc, ds_bahamas, hampdata_loc, drop_time, height, sonde_id
-        ):
+        if not is_complete(ds_dropsonde_loc, hampdata_loc, drop_time, height, sonde_id):
             continue
 
-        # get surface values
         surface_temp = get_surface_temperature(ds_dropsonde_loc)
         surface_ws = get_surface_windspeed(ds_dropsonde_loc)
 
-        # extrapolate dropsonde profiles
         ds_dropsonde_extrap = extrapolate_dropsonde(
             ds_dropsonde_loc, height, ds_bahamas
         )
 
-        plot_dropsonde(ds_dropsonde_extrap)
+        # plot_dropsonde(ds_dropsonde_extrap, ds_dropsonde_loc)
 
-        # convert xarray to ARTS gridded field
         profile_grd = fsm.generate_gridded_field_from_profiles(
             pyarts.arts.Vector(ds_dropsonde_extrap["p"].values),
             ds_dropsonde_extrap["ta"].values,
-            ds_dropsonde_extrap["gpsalt"].values,
+            ds_dropsonde_extrap["alt"].values,
             gases={
                 "H2O": typhon.physics.specific_humidity2vmr(
                     ds_dropsonde_extrap["q"].values
@@ -161,7 +137,6 @@ def calc_arts_bts(date):
             },
         )
 
-        # run arts
         BTs, _ = forward_model(
             ws,
             profile_grd,
@@ -170,55 +145,50 @@ def calc_arts_bts(date):
             height,
         )
 
-        # except (ValueError, KeyError, RuntimeError) as e:
-        #    print(
-        #        f"ARTS or extrapolation failed for dropsonde {sonde_id} with error: {e}, skipping"
-        #    )
-        #    pass
-
-        # Store arts BTs
         TBs_arts[sonde_id] = pd.DataFrame(data=BTs, index=freqs / 1e9)
-        # get according hamp data
-        TBs_hamp[sonde_id] = hampdata_loc.radiometers.TBs.values
+        TBs_hamp[sonde_id] = hampdata_loc["TBs"]
+
+        # plot_TB_comparison(TBs_arts[sonde_id], TBs_hamp[sonde_id], sonde_id)
 
-    # save results
     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")
 
 
 # %% call function
-# date = str(sys.argv[1])
-calc_arts_bts("20240811")
-
-# %% test turning extrapolated dropsonde profiles to ARTS input
-# read config
-cfg = rwfuncs.extract_config_params("config_ipns.yaml")
-ds_dropsonde = xr.open_dataset(cfg["path_dropsondes"], engine="zarr")
-profile = ds_dropsonde.isel(sonde_id=100)
-profile = profile.dropna("gpsalt")
-
-# %% convert xarray to ARTS gridded field
-profile_grd = fsm.generate_gridded_field_from_profiles(
-    pyarts.arts.Vector(profile["p"].values),
-    profile["ta"].values,
-    profile["alt"].values,
-    gases={
-        "H2O": typhon.physics.specific_humidity2vmr(profile["q"].values),
-        "N2": 0.78,
-        "O2": 0.21,
-    },
-)
+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",
+]
+
+flightletters = {}
+for date in dates:
+    flightletters[date] = "a"
 
 # %%
+for date in dates:
+    calc_arts_bts(date, flightletters[date])
 
-BTs, _ = forward_model(
-    ws,
-    profile_grd,
-    10,
-    300,
-    1e4,
-)
-# %%
-freqs = sensor_description[:, 0] + sensor_description[:, 1] + sensor_description[:, 2]
 
 # %%

scripts/process_level_1.py

@@ -0,0 +1,19 @@
+# %%
+import src.readwrite_functions as rwfuncs
+from src.post_processed_hamp_data import PostProcessedHAMPData
+from src.plots_functions import testplot_hamp
+import xarray as xr
+
+# %% load data
+cfg = rwfuncs.extract_config_params("config_ipns.yaml")
+
+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"),
+)
+
+# %% plot
+testplot_hamp(hampdata.radar, hampdata.radiometers)
+
+# %%