ADD: add possibility to retrieve soundings over MCH DWH

Author: wolfidan

src/pyrad_proc/pyrad/io/__init__.py

@@ -111,8 +111,8 @@
     read_ml_ts
     read_windmills_data
     read_radiosounding_wyoming
-    read_radiosounding_igra
-    read_fzl_igra
+    retrieve_radiosounding_igra
+    retrieve_fzl
     read_meteoswiss_xml_vad
 
 Writing data
@@ -255,10 +255,11 @@
 from .read_data_sensor import read_coord_sensors, read_disdro_parsivel  # noqa
 from .read_data_sensor import read_knmi  # noqa
 from .read_data_sensor import (  # noqa
-    read_radiosounding_wyoming,
-    read_radiosounding_igra,
-)
-from .read_data_sensor import read_fzl_igra  # noqa
+    retrieve_radiosounding_wyoming,  # noqa
+    retrieve_radiosounding_igra,  # noqa
+    retrieve_radiosounding_mch,  # noqa
+)  # noqa
+from .read_data_sensor import retrieve_fzl  # noqa
 
 from .read_data_sun import read_sun_hits_multiple_days, read_sun_hits  # noqa
 from .read_data_sun import read_sun_retrieval, read_solar_flux  # noqa

src/pyrad_proc/pyrad/io/read_data_sensor.py

@@ -31,8 +31,8 @@
     read_disdro
     read_disdro_parsivel
     read_radiosounding_wyoming
-    read_radiosounding_igra
-    read_fzl_igra
+    retrieve_radiosounding_igra
+    retrieve_fzl
 """
 
 import os
@@ -58,6 +58,22 @@
 
 urllib3.disable_warnings(InsecureRequestWarning)
 
+COLUMNS_SOUNDING = [
+    "LVLTYPE1",
+    "LVLTYPE2",
+    "ETIME",
+    "PRESS",
+    "PFLAG",
+    "GPH",
+    "ZFLAG",
+    "TEMP",
+    "TFLAG",
+    "RH",
+    "DPDP",
+    "WDIR",
+    "WSPD",
+]
+
 
 def read_windmills_data(fname):
     """
@@ -2517,22 +2533,24 @@ def read_disdro_parsivel(fname):
     return df_disdro
 
 
-def read_radiosounding_wyoming(station, dtime):
+def retrieve_radiosounding_wyoming(station, dtime):
     """
-    Download radiosounding data from the University of Wyoming website.
+    Download radiosounding data from the University of Wyoming website and
+    return as DataFrame with IGRA2-standard columns.
 
-    Parameters:
-        station : str
-            Radiosounding station 5 digits WMO code (e.g., "72776").
-        dtime : datetime.datetime
-            Datetime object specifying the desired date and time.
+    Parameters
+    ----------
+    station : str
+        Radiosounding station 5-digit WMO code (e.g., "72776").
+    dtime : datetime.datetime
+        Datetime object specifying the desired date and time.
 
-    Returns:
-        pandas.DataFrame: A DataFrame containing the radiosounding data with
-        columns: ['PRES', 'HGHT', 'TEMP', 'DWPT', 'RELH', 'MIXR', 'DRCT',
-                  'SKNT', 'THTA', 'THTE', 'THTV]
+    Returns
+    -------
+    pandas.DataFrame or None
+        A DataFrame containing the radiosounding data in IGRA2 format,
+        or None if retrieval was not successful.
     """
-
     base_url = "https://weather.uwyo.edu/cgi-bin/sounding"
     query_params = {
         "region": "naconf",
@@ -2545,36 +2563,68 @@ def read_radiosounding_wyoming(station, dtime):
     }
 
     try:
-        # Try a quick HEAD request to test SSL verification
+        # Quick HEAD request to handle SSL verification issues
         requests.head(base_url, timeout=5)
         verify_setting = True
     except requests.exceptions.SSLError:
         verify_setting = False
 
     response = requests.get(base_url, params=query_params, verify=verify_setting)
-
     if response.status_code != 200:
+        warn(f"Failed to retrieve data: HTTP {response.status_code}")
         return None
 
-    # Extract and parse the data using pandas
+    # Extract text block
     start_idx = response.text.find("PRES")
     end_idx = response.text.find("Station information and sounding indices")
-    data_str = response.text[start_idx:end_idx][0:-10]
+    if start_idx == -1 or end_idx == -1:
+        warn("Could not locate sounding data block in response")
+        return None
 
+    data_str = response.text[start_idx:end_idx].strip()
     data_io = StringIO(data_str)
-    data_df = pd.read_csv(
-        data_io,
-        sep=r"\s+",
-        header=0,
-        skiprows=[1, 2],
-        on_bad_lines="warn",
-    )
+
+    # Parse table (skip line with units)
+    try:
+        data_df = pd.read_csv(
+            data_io,
+            sep=r"\s+",
+            header=0,
+            skiprows=[1],
+            on_bad_lines="skip",
+        )
+    except Exception as e:
+        warn(f"Parsing error: {e}")
+        return None
+
+    # Convert to numeric where possible
     for col in data_df.columns:
-        data_df[col] = pd.to_numeric(data_df[col])
-    return data_df
+        data_df[col] = pd.to_numeric(data_df[col], errors="coerce")
+
+    # Map Wyoming columns to IGRA2 columns
+    df = pd.DataFrame(
+        {
+            "LVLTYPE1": "",  # Not available from Wyoming
+            "LVLTYPE2": "",
+            "ETIME": dtime.strftime("%Y%m%d%H"),  # assign launch time to all rows
+            "PRESS": data_df.get("PRES"),
+            "PFLAG": "",
+            "GPH": data_df.get("HGHT"),
+            "ZFLAG": "",
+            "TEMP": data_df.get("TEMP"),
+            "TFLAG": "",
+            "RH": data_df.get("RELH"),
+            "DPDP": data_df.get("DWPT"),
+            "WDIR": data_df.get("DRCT"),
+            "WSPD": data_df.get("SKNT") * 0.514444,  # knots → m/s
+        },
+        columns=COLUMNS_SOUNDING,
+    )
 
+    return df
 
-def read_radiosounding_igra(station, dtime):
+
+def retrieve_radiosounding_igra(station, dtime):
     """
     Download radiosounding data from  from the Integrated Global Radiosonde Archive
     (IGRA).
@@ -2594,22 +2644,6 @@ def read_radiosounding_igra(station, dtime):
     https://www.ncei.noaa.gov/data/integrated-global-radiosonde-archive/doc/igra2-data-format.txt
     """
     print("Downloading sounding from IGRA database...")
-
-    COLUMNS_SOUNDING = [
-        "LVLTYPE1",
-        "LVLTYPE2",
-        "ETIME",
-        "PRESS",
-        "PFLAG",
-        "GPH",
-        "ZFLAG",
-        "TEMP",
-        "TFLAG",
-        "RH",
-        "DPDP",
-        "WDIR",
-        "WSPD",
-    ]
     COLUMNS_STATION_LIST = [
         "CODE",
         "LATITUDE",
@@ -2752,10 +2786,82 @@ def read_radiosounding_igra(station, dtime):
     return sounding
 
 
-def read_fzl_igra(station, dtime, dscfg=None):
+def retrieve_radiosounding_mch(station, dtime):
     """
-    Get an estimation of the freezing level height from the
-    Integrated Global Radiosonde Archive (IGRA)
+    Retrieve sounding data from MCH DWH service and return as pandas DataFrame.
+    It requires to set the environment variable $JRETRIEVE_PATH that points to a
+    valid jretrieve executable.
+
+    Parameters:
+        station : str
+            Radiosounding station 5 digits WMO code (e.g., "72776").
+        dtime : datetime.datetime
+            Datetime object specifying the desired date and time.
+
+    Returns:
+        pandas.DataFrame: A DataFrame containing the radiosounding data with
+        ['LVLTYPE1', 'LVLTYPE2','ETIME', 'PRESS', 'PFLAG', 'GPH','ZFLAG', 'TEMP',
+        'TFLAG', 'RH','DPDP','WDIR','WSPD']
+
+    """
+    print("Downloading sounding from MCH DWH database...")
+
+    dtime_closest = _closest_sounding_time(dtime)
+    dtime_closest = dtime_closest.strftime("%Y%m%d%H%M%S")
+
+    # Variables to retrieve (DWH IDs mapped to IGRA2 columns)
+    id_variables = "744,742,745,747,748,743"  # p, geopotential, T, Td, ws, wd
+
+    # Build DWH command
+    ret_dwh = (
+        f'{os.path.expandvars("$JRETRIEVE_PATH")} -s profile -w 22 --verbose position '
+        f"-i int_ind,{station} "
+        f"-t {dtime_closest} -p {id_variables}"
+    )
+    try:
+        with os.popen(ret_dwh) as pipefile:
+            lines = pipefile.readlines()
+
+        rows = []
+        for line in lines:
+            if dtime_closest in line:
+                vals = line.split("|")
+                # Extract relevant fields by index (adjust if needed)
+                rows.append(
+                    {
+                        "LVLTYPE1": vals[4].strip(),
+                        "LVLTYPE2": vals[3].strip(),
+                        "ETIME": vals[1].strip(),
+                        "PRESS": float(vals[8].strip()),
+                        "PFLAG": "",  # no flag info available from DWH
+                        "GPH": float(vals[9].strip()),
+                        "ZFLAG": "",
+                        "TEMP": float(vals[10].strip()),
+                        "TFLAG": "",
+                        "RH": "",  # not retrieved
+                        "DPDP": float(vals[11].strip()),
+                        "WDIR": float(vals[13].strip()),
+                        "WSPD": float(vals[12].strip()),
+                    }
+                )
+
+        if not rows:
+            warn(f"No valid data returned by DWH call. Return: {lines}")
+            return None
+
+        df = pd.DataFrame(rows, columns=COLUMNS_SOUNDING)
+        return df
+
+    except EnvironmentError:
+        warn(f"Unable to read DWH data with call {ret_dwh}")
+        return None
+
+
+def retrieve_fzl(station, dtime, dscfg=None, source="IGRA"):
+    """
+    Get an estimation of the freezing level height either the
+    Integrated Global Radiosonde Archive (IGRA), University of Wyoming archive
+    or the MCH digital warehouse (DWH)
 
     Parameters:
         station : str
@@ -2765,10 +2871,11 @@ def read_fzl_igra(station, dtime, dscfg=None):
         dscfg : dictionary of dictionaries, Optional
             If provided will try to read the fzl from the pyrad config
             dictionary instead of fetching it from IGRA
+        source: str
+            either 'IGRA', "WYOMING" or 'MCH'
     Returns:
         fzl : float
-            Interpolated freezing level height (a.s.l) obtained
-            from the nearest in time IGRA sounding
+            Interpolated freezing level height (a.s.l)
     """
     if dscfg is None:
         dscfg = {}
@@ -2780,7 +2887,13 @@ def read_fzl_igra(station, dtime, dscfg=None):
     except (KeyError, TypeError):
         pass
 
-    sounding = read_radiosounding_igra(station, dtime)
+    if source == "IGRA":
+        sounding = retrieve_radiosounding_igra(station, dtime)
+    elif source == "MCH":
+        sounding = retrieve_radiosounding_mch(station, dtime)
+    elif source == "WYOMING":
+        sounding = retrieve_radiosounding_wyoming(station, dtime)
+
     height = sounding["GPH"]
     temp = sounding["TEMP"]
 

src/pyrad_proc/pyrad/proc/process_Doppler.py

@@ -46,7 +46,7 @@
     get_fieldname_pyart,
     convert_pydda_to_pyart_grid,
 )
-from ..io import read_radiosounding_igra
+from ..io import retrieve_radiosounding_igra
 from .process_grid import process_grid
 from ..util import compute_average_vad
 
@@ -1067,7 +1067,7 @@ def process_dda(procstatus, dscfg, radar_list=None):
     sounding = dscfg.get("sounding", None)
     if sounding:
         dtime = pyart.graph.common.generate_radar_time_begin(radar_list[0])
-        sounding_data = read_radiosounding_igra(sounding, dtime)
+        sounding_data = retrieve_radiosounding_igra(sounding, dtime)
     else:
         sounding_data = None
 

src/pyrad_proc/pyrad/proc/process_echoclass.py

@@ -43,7 +43,7 @@
 from ..io.io_aux import get_datatype_fields, get_fieldname_pyart
 from ..io.io_aux import get_file_list, get_datetime
 from ..io.read_data_other import read_centroids
-from ..io.read_data_sensor import read_fzl_igra
+from ..io.read_data_sensor import retrieve_fzl
 
 if importlib.util.find_spec("sklearn_extra") and importlib.util.find_spec("sklearn"):
     _SKLEARN_AVAILABLE = True
@@ -1500,6 +1500,9 @@ def process_hydroclass(procstatus, dscfg, radar_list=None):
             compute the freezing level, if no temperature field name is specified,
             if the temperature field is not in the radar object or if no
             fzl is explicitely defined.
+        sounding_source: str. Dataset keyword
+            Web source from which to get the sounding data, either "IGRA", "WYOMING", or
+            "MCH". If not provided, IGRA is used.
         use_dualpol: Bool. Dataset keyword
             Used with HYDRO_METHOD UKMO. If false no radar data is used and
             the classification is performed using temperature information
@@ -1585,8 +1588,11 @@ def process_hydroclass(procstatus, dscfg, radar_list=None):
                     use_debug=False,
                 )
                 sounding_code = dscfg["sounding"]
+                sounding_source = dscfg.get("sounding_source", "IGRA")
                 t0 = pyart.util.datetime_utils.datetime_from_radar(radar)
-                freezing_level = read_fzl_igra(sounding_code, t0, dscfg=dscfg)
+                freezing_level = retrieve_fzl(
+                    sounding_code, t0, dscfg=dscfg, source=sounding_source
+                )
             else:
                 warn(
                     "iso0 or temperature fields or sounding needed to create hydrometeor "
@@ -1747,8 +1753,11 @@ def process_hydroclass(procstatus, dscfg, radar_list=None):
             freezing_level = dscfg["freezing_level"]
         elif "sounding" in dscfg:
             sounding_code = dscfg["sounding"]
+            sounding_source = dscfg.get("sounding_source", "IGRA")
             t0 = pyart.util.datetime_utils.datetime_from_radar(radar)
-            freezing_level = read_fzl_igra(sounding_code, t0, dscfg=dscfg)
+            freezing_level = retrieve_fzl(
+                sounding_code, t0, dscfg=dscfg, source=sounding_source
+            )
 
         use_dualpol = dscfg.get("use_dualpol", True)
         use_temperature = dscfg.get("use_temperature", True)