Add MTG-IRS

[ADCollard]

MTG-IRS data will arrive as netCDF files of principal components.

In JEDI/UFO we will reconstruct and (optionally) re-apodize the data. So the data ingest phase will involve (again optionally)) thinning the data (as the data volume is huge!) and converting the observations to IODA format.

From Emily:
Here is what we think about how to proceed:
Create YAML for input and output variables (example: iasi.yaml)

Create a Python script to do the following:

Read in YAML

Read data based on the input section of YAML

Generate numpy arrays for each output variable

Create an empty container and add the following (container.add)

variable name

variable data array

data path

category

Call the netCDF encoder

The process involves three parts:

1. Read IRS data into numpy arrays
2. Create a container and add data arrays and other information about the data into the container
3. Encode the data in the container into IODA format and generate an output file.

The first and the third steps are easy. You probably will need some help for the second part. The bufr-query documentation contains the functions you will need to write the converter.
You can also take a look at the existing BUFR converters in the SPOC repository (YAMLs and Python Scripts) to get familiar with the data conversion process.

[jack-woollen]

The process involves three parts:
1 - Read IRS data into numpy arrays
2 - Create a container and add data arrays and other information about the data into the container
3 - Encode the data in the container into IODA format and generate an output file.

Step 1 involves copying IRS netcdf arrays directly into numpy's. Use python, no need for yaml here.
Step 2, use the yaml specs to create a container from numpy's. Need something like bufr.Parser below?
Step 3 call the encoder as in the bufr example below.

Steps 2 and 3 example:
container = bufr.Parser(input_path, YAML_PATH).parse()
datasets = netcdf.Encoder(YAML_PATH).encode(container, OUTPUT_PATH)[category]
obs_temp = dataset["ObsValue/brightnessTemperature"][:]

[ADCollard]

Adding comments from recent e-mail exchange:

From @jack-woollen:
The following files propose a layout for storing irs data into ioda format netcdf files. Hopefully we can discuss at the tag up today.
Jack

irs-ioda-template.txt outline of the ioda template file
iodapy.py.txt python code to create the template file
iodaout.txt ncdump of the template file

From @CoryMartin-NOAA:
Linking the IODA conventions documentation here for future reference:
IODA Conventions
IODA Conventions Tables

From @rmclaren:
Hi Jack,

Thanks for your efforts!

I would start by checking your code into a github repository. Since its really early days that repo can just be a personal repo (under your account) which you can then share with us. This will make it easier for us to look at the code, supply comments, issues etc... Maybe Cory could chime in if he already has an idea for where this code should live (in this case you could branch an existing project).

Question for Others: Do we anticipate processing a lot of different NetCDF files (possibly more in the future) or is this just a one-off?

Comments Python Code:
1) If we do plan to do more NetCDF files it might make sense to make more of a framework where you supply a YAML file with the necessary attributes and automate this process a little more. I attached an example as an illustration (please see the YAML file). The YAML file is for BUFR, but a netcdf converter and YAML could follow the same basic pattern ( Input [bufr:] -> transformations (optional) -> output [encoder:] ). In this way we would not be manually coding each variable and its attributes leading to simpler "user" level code (or no code at all, just YAML).
Others: Please chime in if you think this would be useful, or if it's overkill.
2) It may be that these files are very small (and you can get away with it) but in general I would strongly discourage code like the following:

n=0
d=0
m=2


for i in range(m,159,5):
    for j in range(m,159,5):
        if loca.variables['latitude'][i][j] < 100:
            print(n,i,j,loca.variables['latitude'][i][j])
            pystore(i,j,n,d)
            n=n+1

Looping through all the values like this in Python is extremely slow. You should use NUMPY arrays instead, and use the provided vectorized operations (very fast). Also the ranges for i and j are identical (may be intentional, again I don't know the data).
3) Python imports should all be at the top of the file.

Netcdf comments (see attached example):
1) I'm not familiar with this dataset but usually we don't nest groups in the way that they are nested in the example (ex: ObsValue/mwir/... ).
2) All variables should share the same primary dimension (Location) except for "global" attributes.
3) Global attributes should not be variables on the sub groups, but variables under the "root" group (please see example).
4) Missing dimension variables (see example).
Note: Bufr-Query netcdf encoder handles this stuff, but don't know if you want to go there right now...

Ron

mhs_netcdf_header.txt
bufrtest_mhs_basic_mapping.yaml

From @CoryMartin-NOAA

Hi Ron,

We can chat more this afternoon but to address your first 2 questions:

1. we should put it in SPOC
2. Most of your encoder has focused on BUFR files because that's what we use in GSI, but a lot of data we convert from HDF5/netCDF4 to BUFR, and then convert back to HDF5 for IODA... So, yes, we should probably think of a future where we will be converting HDF5/netCDF4 files into IODA format (in fact we are already doing this for marine and aerosol data using a C++ encoder) https://github.com/NOAA-EMC/obsForge/tree/main/utils/preproc

-Cory

From @jack-woollen

Hi Ron

Thanks for your very thorough comments. My focus for now is just to come up with a template for the IRS data, which is a little different from most of the satellites we work with. The strawman template I made stores data thinned out of 1 IRS file. The instrument makes 25600 soundings every 10 sec, delivered in a package (one netcdf file) called a dwell. There are 280 dwells on the earth, divided into 4 LACs. As far as I can tell it scans the earth about every half hour. We need to store 6 hours of data in a dump file. Some of the data is indexed by location, some is indexed by, and pertinent for, each particular dwell. I'm hoping we can work out the required contents of the IODA/IRS files, and then find the most effective way to put them together. I think the template is just a start in that direction. There's always details to deal with, but first there's this elephant in the room.

Jack

[ADCollard]

@jack-woollen Thanks for this. I think we need to not get too hung up on how to deal with the dwell structure. Any meta data that is dwell-specific should just be included for each individual ob. Or am I missing the issue?

The local PCs complicate matters (as they are per dwell). Furthermore, if we are subsetting so we have only a few obs per dwell, the data compression advantage goes away if we use these. Can I suggest that we ignore the local PCs for now? It is very likely that we never use them.

[ADCollard]

Use this comment to refine output (IODA) format, I have indicated the fields we should keep in bold.

/

netcdf ioda {

dimensions:
location = UNLIMITED ; // (185 currently)
dwell_index = UNLIMITED ; // (1 currently)
global_pc = 150 ;
local_pc = 5 ;
mw_number = 1079 ;
lw_number = 881 ;

group: ObsValue {

group: mwir {
variables:
int global_pc_scores(location, global_pc) ;
ubyte global_pc_comp_err(dwell_index) ; (This is the overflow error flag - we probably just want to filter these out)
float global_pcr_scores(location) ;
ubyte global_pcrs_quality(location) ;
int local_pc_scores(location, local_pc) ;
ubyte local_pc_comp_err(dwell_index) ;
float local_eigenvalues(dwell_index, local_pc) ;
float local_pcr_scores(location) ;
float residual_energy(dwell_index) ;
float local_pcr_operator(dwell_index, mw_number, local_pc) ;
ushort dc_image_radiance(location) ; (I am not sure what this is!)
short spatial_sample_quality(location) ;
byte detector_sample_quality(location) ;
} // group mwir

group: lwir {
variables:
int global_pc_scores(location, global_pc) ;
ubyte global_pc_comp_err(dwell_index) ; (This is the overflow error flag - we probably just want to filter these out)
float global_pcr_scores(location) ;
ubyte global_pcrs_quality(location) ;
int local_pc_scores(location, local_pc) ;
ubyte local_pc_comp_err(dwell_index) ;
float local_eigenvalues(dwell_index, local_pc) ;
float local_pcr_scores(location) ;
float residual_energy(dwell_index) ;
float local_pcr_operator(dwell_index, lw_number, local_pc) ;
ushort dc_image_radiance(location) ; (I am not sure what this is!)
short spatial_sample_quality(location) ;
byte detector_sample_quality(location) ;
} // group lwir
} // group ObsValue

group: MetaData {
variables:
double time(dwell_index) ;
ubyte dwell_type(dwell_index) ;
ushort dwell_number(dwell_index) ;
ubyte stroke_direction(dwell_index) ;
ushort latitude(location) ;
ushort longitude(location) ;
short satellite_azimuth_angle(location) ;
short satellite_zenith_angle(location) ;
short solar_azimuth_angle(location) ;
short solar_zenith_angle(location) ;
ubyte dust_warning(location) ;
ushort cloud_signal(location) ;
ubyte cloud_fraction(location) ;

group: mwir { The first part of this section is probably not needed at all as these are fixed properties
variables:
float spatial_sampling_distance ;
float start_wavenumber ;
float end_wavenumber ;
float wavenumber_step ;
float wavenumber(mw_number) ;
int start_number ;
int end_number ;
float maximum_opd ;
Per the data guide, the following are not implemented in IDPF-S v2.1
ubyte instrument_quality_warning(location) ;
ubyte radiometric_quality_warning(location) ;
ubyte spectral_quality_warning(location) ;
ubyte rad_rts_noise_warning(location) ;
ubyte phase_average_warning(location) ;
ubyte phase_std_warning(location) ;
int number_of_expected_earthview_spatial_samples(dwell_index) ;
int number_of_valid_earthview_spatial_samples(dwell_index) ;
int number_of_spaceview_spatial_samples(dwell_index) ;
int number_of_saturation_warning_spatial_samples(dwell_index) ;
int number_of_excluded_detector_sub_pixels(dwell_index) ;
int number_of_invalid_detector_sub_pixels(dwell_index) ;
int number_of_straylight_warning_spatial_samples(dwell_index) ;
int number_of_straylight_correction_warning_spatial_samples(dwell_index) ;
} // group mwir

group: lwir { The first part of this section is probably not needed at all as these are fixed properties
variables:
float spatial_sampling_distance ;
float start_wavenumber ;
float end_wavenumber ;
float wavenumber_step ;
float wavenumber(lw_number) ;
int start_number ;
int end_number ;
float maximum_opd ;
Per the data guide, the following are not implemented in IDPF-S v2.1
ubyte instrument_quality_warning(location) ;
ubyte radiometric_quality_warning(location) ;
ubyte spectral_quality_warning(location) ;
ubyte rad_rts_noise_warning(location) ;
ubyte phase_average_warning(location) ;
ubyte phase_std_warning(location) ;
int number_of_expected_earthview_spatial_samples(dwell_index) ;
int number_of_valid_earthview_spatial_samples(dwell_index) ;
int number_of_spaceview_spatial_samples(dwell_index) ;
int number_of_saturation_warning_spatial_samples(dwell_index) ;
int number_of_excluded_detector_sub_pixels(dwell_index) ;
int number_of_invalid_detector_sub_pixels(dwell_index) ;
int number_of_straylight_warning_spatial_samples(dwell_index) ;
int number_of_straylight_correction_warning_spatial_samples(dwell_index) ;
} // group lwir
} // group MetaData

group: ObsError {
} // group ObsError

group: PreQc {
} // group PreQc
}

[CoryMartin-NOAA]

What does this file look like? https://github.com/JCSDA-internal/ufo-data/blob/develop/testinput_tier_1/PCscores_transform_iasi_2024120112.nc4

netcdf PCscores_transform_iasi_2024120112 {
dimensions:
	Location = 5 ;
	Channel = 15 ;
variables:
	int64 Location(Location) ;
		Location:_FillValue = -9223372036854775801LL ;
	int Channel(Channel) ;
		Channel:_FillValue = -2147483643LL ;

// global attributes:
		string :_ioda_layout = "ObsGroup" ;
		:_ioda_layout_version = 0 ;
		:history = "Wed Jun 11 09:50:05 2025: ncks --fix_rec_dmn Location -d Location,0,4 jopa_iasi_20241201T1200Z_out.rr.nc4 jopa_iasi_20241201T1200Z_out.rr.nc4.snip" ;
		:NCO = "netCDF Operators version 4.7.5 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ;

group: MetaData {
  variables:
  	int64 dateTime(Location) ;
  		dateTime:_FillValue = -9223372036854775801LL ;
  		string dateTime:units = "seconds since 1970-01-01T00:00:00Z" ;
  	float latitude(Location) ;
  		latitude:_FillValue = -3.368795e+38f ;
  	float longitude(Location) ;
  		longitude:_FillValue = -3.368795e+38f ;
  	float principalComponentScore1(Location) ;
  		principalComponentScore1:_FillValue = -3.368795e+38f ;
  	float principalComponentScore10(Location) ;
  		principalComponentScore10:_FillValue = -3.368795e+38f ;
  	float principalComponentScore100(Location) ;
  		principalComponentScore100:_FillValue = -3.368795e+38f ;
  	float principalComponentScore101(Location) ;
  		principalComponentScore101:_FillValue = -3.368795e+38f ;
                   |
                   |
                   |
                   |
  		principalComponentScore98:_FillValue = -3.368795e+38f ;
  	float principalComponentScore99(Location) ;
  		principalComponentScore99:_FillValue = -3.368795e+38f ;
  	int satelliteIdentifier(Location) ;
  		satelliteIdentifier:_FillValue = -2147483643 ;
  	int scanLineNumber(Location) ;
  		scanLineNumber:_FillValue = -2147483643 ;
  	float sensorAzimuthAngle(Location) ;
  		sensorAzimuthAngle:_FillValue = -3.368795e+38f ;
  	int sensorScanPosition(Location) ;
  		sensorScanPosition:_FillValue = -2147483643 ;
  	float sensorZenithAngle(Location) ;
  		sensorZenithAngle:_FillValue = -3.368795e+38f ;
  } // group MetaData

[rmclaren]

@CoryMartin-NOAA Did you mean to post the IASI PCA dataset information to this thread?

[CoryMartin-NOAA]

@rmclaren yes, just as an example of the PCA packing that the UKMO is doing (I'm surprised they aren't using multidimensional arrays...)

[rmclaren]

Thanks for the context in that slide @jackwoolen. I'm just thinking out loud here. So whenever you process one type of file into another format you need an intermediate format to bridge between the two formats. I think that format should be xarray (https://xarray.dev/). Not only does it efficiently support multiple dimensional arrays but there are a lot of extra libraries and tools built on it (it opens a lot of possibilities). XArray has native support for NetCDF and ZARR as well so it would be very easy to use it for this purpose.

[jack-woollen]

@jack-woollen Thanks for this. I think we need to not get too hung up on how to deal with the dwell structure. Any meta data that is dwell-specific should just be included for each individual ob. Or am I missing the issue?

The local PCs complicate matters (as they are per dwell). Furthermore, if we are subsetting so we have only a few obs per dwell, the data compression advantage goes away if we use these. Can I suggest that we ignore the local PCs for now? It is very likely that we never use them.

@ADCollard The lists below represent data items which are not stored by location in the irs file, so need special treatment.
I put marks in the margin to indicate dispositions which seem reasonable from our discussion today.

G - global scalar variables to be included only once in each file
L - dwell specific information which should be indexed by location
X - unnecessary information which can be ignored

G float spatial_sampling_distance ;
G float start_wavenumber ;
G float end_wavenumber ;
G float wavenumber_step ;
G float wavenumber(wavenumber) ;
G int start_number ;
G int end_number ;
G float maximum_opd ;

L double time ;
L ubyte dwell_type ;
L ushort dwell_number ;
L ubyte stroke_direction ;

L ubyte global_pc_comp_err ;
L float residual_energy ;
X ubyte local_pc_comp_err ;
X float local_eigenvalues(local_pc) ;
X float local_pcr_operator(wavenumber, local_pc) ;

L int number_of_expected_earthview_spatial_samples ;
L int number_of_valid_earthview_spatial_samples ;
L int number_of_spaceview_spatial_samples ;
L int number_of_saturation_warning_spatial_samples ;
L int number_of_excluded_detector_sub_pixels ;
L int number_of_invalid_detector_sub_pixels ;
L int number_of_straylight_warning_spatial_samples ;
L int number_of_straylight_correction_warning_spatial_samples ;

[CoryMartin-NOAA]

See https://github.com/NOAA-EMC/spoc/blob/feature/simple-example/dump/scripts/atmosphere/example_for_jack.py
which produces:

netcdf out {
dimensions:
        Location = 100 ;
variables:
        int Location(Location) ;
                Location:_FillValue = 2147483647 ;

// global attributes:
                :dummy = "dummy str" ;
data:

 Location = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0 ;

group: ObsValue {
  variables:
        double sine(Location) ;
                sine:_FillValue = 1.79769313486232e+308 ;
                sine:long_name = "Sine Wave" ;
                sine:units = "unitless" ;
  data:

   sine = 0, 0.0634239196565645, 0.126592453573749, 0.18925124436041, 
      0.251147987181079, 0.312033445698487, 0.371662455660328, 
      0.429794912089172, 0.486196736100469, 0.540640817455598, 
      0.59290792905464, 0.642787609686539, 0.690079011482112, 
      0.734591708657533, 0.776146464291757, 0.814575952050336, 
      0.849725429949514, 0.881453363447582, 0.909631995354518, 
      0.934147860265107, 0.954902241444074, 0.971811568323542, 
      0.984807753012208, 0.993838464461254, 0.998867339183008, 
      0.999874127673875, 0.996854775951942, 0.989821441880933, 
      0.978802446214779, 0.963842158559942, 0.945000818714668, 
      0.922354294104581, 0.895993774291336, 0.866025403784438, 
      0.832569854634771, 0.795761840530832, 0.755749574354258, 
      0.712694171378863, 0.666769000516292, 0.618158986220605, 
      0.567059863862771, 0.513677391573406, 0.458226521727411, 
      0.400930535406614, 0.342020143325669, 0.28173255684143, 
      0.220310532786541, 0.15800139597335, 0.0950560433041824, 
      0.0317279334980677, -0.0317279334980679, -0.0950560433041826, 
      -0.15800139597335, -0.220310532786541, -0.28173255684143, 
      -0.342020143325669, -0.400930535406614, -0.45822652172741, 
      -0.513677391573406, -0.567059863862771, -0.618158986220605, 
      -0.666769000516292, -0.712694171378863, -0.755749574354258, 
      -0.795761840530832, -0.832569854634771, -0.866025403784439, 
      -0.895993774291336, -0.922354294104581, -0.945000818714668, 
      -0.963842158559942, -0.978802446214779, -0.989821441880933, 
      -0.996854775951942, -0.999874127673875, -0.998867339183008, 
      -0.993838464461254, -0.984807753012208, -0.971811568323542, 
      -0.954902241444074, -0.934147860265107, -0.909631995354518, 
      -0.881453363447582, -0.849725429949514, -0.814575952050336, 
      -0.776146464291757, -0.734591708657533, -0.690079011482112, 
      -0.64278760968654, -0.59290792905464, -0.540640817455597, 
      -0.486196736100469, -0.429794912089172, -0.371662455660327, 
      -0.312033445698487, -0.251147987181079, -0.18925124436041, 
      -0.126592453573749, -0.0634239196565645, -2.44929359829471e-16 ;
  } // group ObsValue
}

[jack-woollen]

@CoryMartin-NOAA I like that! I have this problem:

Traceback (most recent call last):
File "/scratch3/NCEPDEV/global/Jack.Woollen/IRSPP/IRSPP/IRSPPv1.3_test_cases/output/irs2ioda/container/container.py", line 2, in
import bufr
ModuleNotFoundError: No module named 'bufr'

[CoryMartin-NOAA]

@jack-woollen you have to have an ObsForge build available, here's a branch with the proper SPOC hash + changes to load modules:
https://github.com/NOAA-EMC/obsForge/tree/feature/simple-spoc-example

git clone -b feature/simple-spoc-example --recursive https://github.com/NOAA-EMC/obsForge
cd obsForge
./build.sh
source ush/of_setup.sh
cd sorc/spoc/dump/scripts/atmosphere/
python example_for_jack.py

Let me know if you run into any issues.