Modidy ioda_obs_space to add cycle time

[EdwardSafford-NOAA]

Description

Several types of monitoring plots need to use the select time transform to set data groupings of 1, 4, and 120 cycles. To support this we need to add a dimension and coordinate of Time to the collection for each loaded data file. The value of Time coordinate will be derived from the file name.

[EdwardSafford-NOAA]

Found some interesting things in data_collections.py and wanted to make sure I'm reading this correctly.

I ran into a problem when I was testing my changes to ioda_obs_space.py with testIodaObsSpaceAircraft.yaml. When I added the Time dimension I ran into this problem:

File "/scratch1/NCEPDEV/da/Edward.Safford/noscrub/eva/opt/eva/data/data_collections.py", line 394, in  nan_float_values_outside_threshold  
data_var_value[np.abs(data_var_value) > threshold] = np.nan
ValueError: assignment destination is read-only

That's because when the first dataset is added to the data_collections object it's done so thusly:

if collection_name not in self._collections:
            self._collections[collection_name] = collection.copy(deep=False)

That shallow copy means we're getting a view in collection, not a writeable version.

So when we get to data_collections.nan_float_values_outside_threshold(threshold) it does this:

# Loop over the variables and set to nan outside of threshold
# -----------------------------------------------------------
for group_variable in groups_variables:

    # Split name into group and variable
    [group, variable] = group_variable.split('::')

    # Get the data
    data_var_value = self.get_variable_data(collection, group, variable)
   
    # For float data sceen outside threshold
    if 'float' in str(data_var_value.dtype):
        data_var_value[np.abs(data_var_value) > threshold] = np.nan

data_var_value is a slice of a view and isn't writeable. Apparently if the dataset is a single dimension we get away with it owing to some sort of subtle dark force deep in the bowels of Xarray.

Edited to remove 2nd question that I found an answer for.

I'm looking in to how I can make data_var_value writeable. I think that's a better solution than a deep copy for everything.

@CoryMartin-NOAA, @Dooruk , @kevindougherty-noaa @ADCollard

[EdwardSafford-NOAA]

I've found two ways to fix the problem in data_collections.py.

First way is in create_or_add_to_collection to use a deep copy when creating a collections object. This makes all the variable data writeable.

 88         if collection_name not in self._collections:
 89             self._collections[collection_name] = collection.copy(deep=True)

The second solution is to modify the get_variable_data method to make the returned array writeable. This avoids making the non-float data writeable:

303         # Extract the data array
304         variable_array = variable_array.data.copy()
305         variable_data = variable_array.data

I don't see a lot of difference in timing. They are both all over the place ranging from ~16s to > 40s.

[EdwardSafford-NOAA]

Additional timing tests demonstrate a consistent advantage to using the deep copy, so I'm going to go that way.

[CoryMartin-NOAA]

How many deep copies are needed? This doesn't significantly increase memory usage does it?

[kevindougherty-noaa]

@EdwardSafford-NOAA just curious, how does that affect memory with a deep copy? It may be fine with lower res datasets but if we add a hyperspectral dataset I worry memory will become an issue.

[CoryMartin-NOAA]

Ha, @kevindougherty-noaa and I with the same question in parallel

[EdwardSafford-NOAA]

Not sure. I can run some more tests with the available test yaml files. Any one(s) in particular you'd like to have me try?

[CoryMartin-NOAA]

The test datasets might not be able to show much, it may require one of the files I staged previously that were ops GSI files converted to IODA format. Maybe a CrIS or IASI?

[EdwardSafford-NOAA]

After some additional testing you can scratch all my timing conclusions. The test runs are all over the place and unless I ran a few hundred I don't think it's possible to conclude which solution is faster.

In my simple tests using CrIS data from your sample_ioda fileset I don't see any memory difference, but my yaml script is pretty simple. I'm going to reverse myself, err on the side of caution, and use the 'make writable only as needed' solution to move forward.

I'm still intrigued how it was that we were getting away with write access using a shallow copy in a single dimension dataset. ChatGPT says that's a known 'feature' of Xarray:DataSet, so take that for what it's worth.

[EdwardSafford-NOAA]

Ratz. The 'make writeable as needed' solution breaks mon_data_space. The copy() used to accomplish getting writeable data runs afoul of a python limitation: datetime64 values aren't natively representable in a raw memory buffer (mon_data_space encodes cycle_tm as datetime64). mon_data_space does work with the deep copy though. Before I look into how mon_data_space might be able to work around this I'm going to run all the test scripts that use other data objects so see the exact dimensions of what I'm working with.

[EdwardSafford-NOAA]

I determined that theses objects -- csv_space, jedi_variational_bias_correction, and mon_data_space -- all contain datetime64 datatypes and therefore failed my 'writeable' fix as per above. The solution I've found is to replace

# Extract the actual data array
variable_array = variable_array.data.copy()
variable_data = variable_array.data

with a check on the datatype. data_collections only tries to overwrite float* types (replacing them with NAN if they are out of range), so a slice works for everything else.

 # Extract the actual data array
 if np.issubdtype(variable_array.dtype, np.floating):
    variable_array = variable_array.data.copy()
    variable_data = variable_array.data
 else:
    variable_data = variable_array.data

So that's my working solution and I'm back to working this issue proper.

[karpob]

Just out of curiousity is there any reason you must have time from the filename? If we're talking about the IODA obs files, each ob is tagged with a time that can be decoded to a cycle YYYY-MM-DD-THH:MM:SS. If you could manage to encode time, I would imagine you could potentially do a timeseries with a variable window to see how obs are being used within a 6 hour cycle, which could be useful. I would think xarray should be able to handle the widowing/averaging for 6 hours or any window you choose.

[EdwardSafford-NOAA]

What you describe would certainly be useful and we might well eventually add something like that for an in-depth analysis of specific instruments. I'm working on reproducing existing graphics from the legacy DA monitoring suite which use various groupings of cycle data in the summary and time series plots. These plots don't break things down any further than cycle time. Cycle time for each file could be derived from the obs dateTime but for this purpose it seemed more straightforward to just grab it from the filename.

[karpob]

Yeah, I assued. My thought/2 cents is it might be easier to not try to bend xarray to your will and add it from the filename, but grab it using xarray and manipulate xarray after the fact.