Work with JEDI IODA files

Work with JEDI IODA files#

IODA stands for Interface for Observational Data Access.

JEDI can only ingest observatios stored in the IODA format. Therefore, observations from other formats (such as text, PrepBUFR, BUFR, etc.) must first be converted into IODA files before they can be used by JEDI. These observation files are usually referred to as input IODA files.

After data assimilation (DA), JEDI writes both the original observation information and DA-related quantities (e.g., innovations, H(x), quality control flags, and bias corrections) into output IODA files.
In the MPAS and JEDI based wx-workflow (including rrfs-workflow, rtma-workflow, ar-worklfow, etc), these output IODA files are referred to as jdiag files.

This notebook demonstrates how to process both input and output IODA files using the netCDF4 Python package as well as the DAmonitor obsSpace class, which provides a more convenient and concise interface than working directly with the raw netCDF4 API.

0. Environment Setup#

%%time 

# autoload external python modules if they changed
%load_ext autoreload
%autoreload 2

import os, sys
pyDAmonitor_ROOT=os.getenv("pyDAmonitor_ROOT")
if pyDAmonitor_ROOT is None:
    print("!!! pyDAmonitor_ROOT is NOT set. Run `source ush/load_pyDAmonitor.sh`")
else:
    print(f"pyDAmonitor_ROOT={pyDAmonitor_ROOT}\n")
sys.path.insert(0, pyDAmonitor_ROOT)
    
from DAmonitor.base import query_dataset, query_data, query_obj

pyDAmonitor_ROOT=/ncrc/home2/Guoqing.Ge/pyDAmonitor

CPU times: user 7.26 s, sys: 117 ms, total: 7.38 s
Wall time: 320 ms

1. Use the raw netCDF4 interface to process ioda files#

1.1 Inspecting the Structure of IODA Files#

%%time
from netCDF4 import Dataset

dcObs={
    "obs": "aircar",
    "jdiag": "t133",
      }
print("-------- obs_file --------")
obs_file = f"{pyDAmonitor_ROOT}/data/ioda/ioda_{dcObs['obs']}.nc"
ds = Dataset(obs_file, mode='r')
print(ds)
query_dataset(ds)     # query_dataset is a function provided by DAmonitor to loop through and print all groups/variables information of a NetCDF4 dataset

print("\n\n-------- jdiag_file --------")
jdiag_file = f"{pyDAmonitor_ROOT}/data/mpasjedi/jdiag_{dcObs['obs']}_{dcObs['jdiag']}.nc"
ds = Dataset(jdiag_file, mode='r')
print(ds)
query_dataset(ds)

-------- obs_file --------

<class 'netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
    _ioda_layout: ObsGroup
    _ioda_layout_version: 0
    Orig_obs_file: ioda_aircar.nc
    Grid_file: /scratch3/BMC/rtwrfruc/RRFSv2X/code/rrfs-workflow.20260611/fix/conus3km/conus3km.static.nc
    Shrink_factor: 0.005
    dimensions(sizes): Location(65847)
    variables(dimensions): int64 Location(Location)
    groups: MetaData, ObsError, ObsType, ObsValue, QualityMarker
*MetaData
    aircraftFlightNumber, aircraftFlightPhase, dateTime, dumpReportType, height, latitude, longitude, prepbufrDataLvlCat, prepbufrReportType, pressure, stationElevation, stationIdentification, timeOffset, longitude_latitude_pressure
*ObsError
    airTemperature, pressure, relativeHumidity, specificHumidity, windEastward, windNorthward
*ObsType
    airTemperature, specificHumidity, windEastward, windNorthward
*ObsValue
    airTemperature, dewpointTemperature, specificHumidity, windEastward, windNorthward
*QualityMarker
    airTemperature, height, pressure, specificHumidity, windEastward, windNorthward
Location


-------- jdiag_file --------
<class 'netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
    Grid_file: /scratch3/BMC/rtwrfruc/RRFSv2X/code/rrfs-workflow.20260611/fix/conus3km/conus3km.static.nc
    Orig_obs_file: ioda_aircar.nc
    Shrink_factor: 0.005
    _ioda_layout: ObsGroup
    _ioda_layout_version: 0
    dimensions(sizes): Location(27732)
    variables(dimensions): int64 Location(Location)
    groups: EffectiveError0, EffectiveError1, EffectiveError2, EffectiveQC0, EffectiveQC1, EffectiveQC2, MetaData, ObsBias0, ObsBias1, ObsBias2, ObsError, ObsType, ObsValue, QualityMarker, TempObsErrorData, hofx0, hofx1, hofx2, innov1, oman, ombg
*EffectiveError0
    airTemperature
*EffectiveError1
    airTemperature
*EffectiveError2
    airTemperature
*EffectiveQC0
    airTemperature
*EffectiveQC1
    airTemperature
*EffectiveQC2
    airTemperature
*MetaData
    stationIdentification, aircraftFlightPhase, timeOffset, prepbufrDataLvlCat, longitude, stationElevation, latitude, prepbufrReportType, pressure, longitude_latitude_pressure, dateTime, height, dumpReportType, aircraftFlightNumber
*ObsBias0
    airTemperature
*ObsBias1
    airTemperature
*ObsBias2
    airTemperature
*ObsError
    airTemperature, relativeHumidity, specificHumidity, windEastward, pressure, windNorthward
*ObsType
    specificHumidity, windEastward, airTemperature, windNorthward
*ObsValue
    dewpointTemperature, specificHumidity, airTemperature, windEastward, windNorthward
*QualityMarker
    airTemperature, pressure, windEastward, specificHumidity, height, windNorthward
*TempObsErrorData
    airTemperature
*hofx0
    airTemperature
*hofx1
    airTemperature
*hofx2
    airTemperature
*innov1
    airTemperature
*oman
    airTemperature
*ombg
    airTemperature
Location
CPU times: user 11.6 ms, sys: 7.99 ms, total: 19.6 ms
Wall time: 221 ms

Understanding the IODA File Structure#

An IODA file is a standard NetCDF4 file and can be opened directly with ds = Dataset(obs_file, mode="r").

Calling print(ds) provides a high-level summary of the file, including its dimensions (e.g., 65,847 locations in the input observation file, and 28,039 locations in the jdiag file), groups, and variables.

From the output above, we can see that an input IODA file (observation file) contains groups such as MetaData, ObsError, ObsType, ObsValue, QualityMarker along with dimensions such as Location (and Channel for satellite observations).

After data assimilation, JEDI writes an output IODA file (commonly called a jdiag file in wx-workflow), which contains the original observation information plus several additional groups that store DA diagnostics, including EffectiveErrorX, EffectiveQCX, ObsBiasX, hofxX, innov1, oman, ombg etc.

These groups have the following meanings:

ombg:  observation - background
oman:  observation - analysis
innova1:  observation - first_guess after the first outer loop

# X = 0, 1, 2, ... indicates the DA outer loops: 
# 0: before any outer loop 
# 1: after the first outer loop 
# 2: after the second outer loop 
# ...

EffectiveErrorX:  effective obs error used in the inner minimization process, see the above for the meaning of X (0, 1, 2)
EffectiveQCX:  effective QC flag
ObsBiasX:  observation bias
hofxX:  H(x), the non-linear forward operatior result, i.e. the model counterpart to the observation

An input IODA file often contains observations for multiple variables (e.g., t, q, uv, etc.),
whereas a jdiag file typically contains DA diagnostics for a single variable. For example, jdiag_aircar_t133.nc contains diagnostics only for temperature observations.

A jdiag file may also contains diagnostics for multiple variables if they are assimilated within the same ObsSpace. However, in wx-workflow (including rrfs-workflow, rtma-workflow, ar-workflow, etc), the current convention is to configure one ObsSpace for each PrepBUFR observation type (i.e., one entry in the convinfo file) or for each satellite sensor. As a result, most jdiag files produced by these workflows contain diagnostics for only one single observation variable.

1.2 Access the metadata and acutal observation data#

# print the "Location" dimension and its actual content
print(ds.dimensions['Location'])
print(len(ds.variables['Location']))
print(ds.variables['Location'][:])

"<class 'netCDF4.Dimension'>" (unlimited): name = 'Location', size = 27732
27732
[    0  8555  8557 ... 30767 37440 41184]

# print the ObsValue of t (airTemperature) and other variables
print(ds.groups["ObsValue"].variables["airTemperature"])
print("----------------")
print(ds.groups["ObsValue"].variables["airTemperature"][:])  # actual content
print(ds.groups["ObsValue"].variables["specificHumidity"][:])

<class 'netCDF4.Variable'>
float32 airTemperature(Location)
    _FillValue: -3.3687953e+38
    coordinates: longitude latitude
    long_name: Temperature
    units: K
path = /ObsValue
unlimited dimensions: Location
current shape = (27732,)
filling on
----------------
[279.25 273.55 279.25 ... 282.05 292.75 304.25]
[-- -- -- ... -- -- --]

Note#

The file jdiag_aircar_t133.nc contains data assimilation (DA) results only for the aircar_t133 ObsSpace. Consequently, only ObsValue/airTemperature contains valid observation values. All other variables under the ObsValue group contain only missing values (i.e., [-- -- -- ... -- -- --]).

2. Use the DAmonitor class to process ioda files#

As shown above, accessing metadata and variables in an IODA file using the netCDF4 interface requires navigating through the group/variable hierarchy.
For example, to read the air temperature observations, you would write: ds.groups["ObsValue"].variables["airTemperature"][:]
This expression contains 52 characters and can become cumbersome when we access data repeatedly.

The DAmonitor class simplifies this process by flattening the IODA group/variable structure and providing direct access to commonly used metadata and variables.
For example, instead of writing print(ds.groups["ObsValue"].variables["specificHumidity"][:]), we may simply write print(t133.q.ObsValue)
This reduces the access expression from 52 characters to just a few characters, making the code more concise, readable, and easier to maintain.

2.1 DAmonitor Overview#

DAmonitor consists of several modules, including obs, mpas, base, and shapes, each serving a different purpose:

obs: Provides classes and functions for observation-space processing, such as obsSpace and fit_rate.
mpas: Provides classes and functions for model-space processing, including visualization utilities such as hcross_contour and vcross_contour.
base: Provides common utilities, including source, get_run_directory, get_inv_bkg_ana_files, load_inv_bkg_ana, query_dataset, query_data, query_obj, to_dataframe, load_figs, and aircraft_reject_list_to_df.
shapes: Provides geographic features to facilitate plotting, including coast_lines, state_lines, counties_gdf, and county_lines.

In this notebook, we focus on the obsSpace class from the obs module.

We can easily load an IODA file into an obsSpace object using obs = obsSpace(filename) and then inspect the contents (attributes, data, etc) of the object using query_obj(obs)

As shown below, an obsSpace object may contain the following attributes (the exact set depends on the type of observations in the IODA file):

Attribute	Description
`filepath`	The input file path.
`ds`	The original `netCDF4.Dataset` object (`self.ds = Dataset(filepath, mode='r')`). You can access the dataset directly just as in Step 1 if needed.
`nlocs`	The number of locations (records) in the IODA file.
`locations`	The indices of all locations.
`nchannels`	The number of channels for a satellite sensor.
`channels`	The channel numbers.
`groups`	The group names contained in the IODA file.
`metadata`	Stores all variables from the original `Metadata` group for convenient access, in addition to direct dataset access.
`t`, `td`, `tv`, `ts`, `u`, `v`, `q`, `rh`, `ps`, `bt`, `rad`, etc.	Observation data for the corresponding variable. These attributes are available only if the variable exists in the IODA file (e.g., `rad` (radiance) is present only in satellite IODA files).
Observation data attributes	Each observation data in the above may contain the following attributes: all variables from the `Metadata` group, along with `ObsError`, `ObsType`, `ObsValue`, `QualityMarker`, `EffectiveErrorX`, `EffectiveQCX`, `ObsBiasX`, `hofxX`, `innov1`, `oman`, `ombg`, and others, depending on the contents of the IODA file.

# use obsSpace to access an ioda file

from DAmonitor.obs import obsSpace
t133=obsSpace(jdiag_file)   # t133 is an obsSpace object
print(t133.filepath)
query_obj(t133)

/ncrc/home2/Guoqing.Ge/pyDAmonitor/data/mpasjedi/jdiag_aircar_t133.nc

['ds',
 'filepath',
 'groups',
 'locations',
 'metadata',
 'nlocs',
 'q',
 'rh',
 't',
 'td',
 'u',
 'v']

# query the data  (note: t is an internal object containing a `data` attribute and query_data is intelligent enough to query the `data` attribute directly)

query_data(t133.t)

EffectiveError0, EffectiveError1, EffectiveError2, EffectiveQC0, EffectiveQC1, EffectiveQC2, stationIdentification, aircraftFlightPhase, timeOffset, prepbufrDataLvlCat, longitude, stationElevation, latitude, prepbufrReportType, pressure, dateTime, height, dumpReportType, aircraftFlightNumber, ObsBias0, ObsBias1, ObsBias2, ObsError, ObsType, ObsValue, QualityMarker, TempObsErrorData, hofx0, hofx1, hofx2, innov1, oman, ombg, Location

# access the actual data

print(t133.t.stationIdentification[10])
print(t133.t.aircraftFlightNumber[10])
print(t133.t.aircraftFlightNumber)
print(t133.t.ObsValue)

t133.t.ObsValue

VTTDSIZA
C4GBTSZA
['ISQBOZRA' '0BIDSQRA' '0BIDSQRA' ... 'EPGCFGRA' 'EPGCFGRA' 'EPGCFGRA']
[279.25 273.55 279.25 ... 282.05 292.75 304.25]

masked_array(data=[279.25, 273.55, 279.25, ..., 282.05, 292.75, 304.25],
             mask=False,
       fill_value=np.float64(1e+20),
            dtype=float32)

# check what's inside EffectiveQc and EffectiveError

print(t133.t.EffectiveQC0)
print(t133.t.EffectiveError0)

[0 0 0 ... 0 0 0]
[0.8210167288780212 0.7708320617675781 0.7883864045143127 ...
 0.827965259552002 1.0205713510513306 7.202971935272217]

# check whether EffectiveError0, EffectiveError1 and EffectiveError2 are the same
# Note: They are same for this ioda file and this DA configuration. They may be different in other situations.

import numpy as np
a = t133.t.EffectiveError0
b = t133.t.EffectiveError1
c = t133.t.EffectiveError2
print(np.array_equal(a, b))
print(np.array_equal(a, c))

True
True

# check whether EffectiveQC0, EffectiveQC1 and EffectiveQC2 are the same
# Note: They are same for this ioda file and this DA configuration. They may be different in other situations.

import numpy as np
a = t133.t.EffectiveQC0
b = t133.t.EffectiveQC1
c = t133.t.EffectiveQC2
print(np.array_equal(a, b))
print(np.array_equal(a, c))

True
True

3. What method to use? DAmonitor or the raw method?#

DAmonitor is recommended for working with IODA files, as it provides a higher-level interface that simplifies common tasks and improves code readability. However, the raw method introduced in Section 1 is also accepted.