v0.3.0 Testing Feedback

[tomvothecoder]

You can post any feedback, questions, or comments in this post. GitHub discussions is threads based, so it will be easy to follow the individual conversations.

UPDATE (8/31/22)
Just an FYI that v0.3.1 was released which includes additional bug fixes that you might have come across in v0.3.0.
Check out the GitHub discussions post here: #322

[msahn]

Feedback

I have tested xCDAT v0.3.0 following the test instruction. A metric measuring precipitation variability across timescales (code, demo), one of the metrics implemented in PMP, has been converted to use xCDAT instead of CDAT.

Below is the library script I modified to convert from CDAT to xCDAT. The APIs are changed for open&write nc files, horizontal_regridding, data concatenating, land-sea masking, spatial average, and others.

• library using CDAT
• library using xCDAT

The differences in JSON results are negligibly small when I test it with demo data.

• using CDAT: demo notebook, JSON result
• using xCDAT: demo notebook, JSON result

Question

While converting my script to use xCDAT, I noticed that most xCDAT APIs are based on "xarray.Dataset". Could it be possible to make it based on "xarray.DataArray" as well?

[tomvothecoder]

Thank you @msahn for posting your experimentation results here!

While converting my script to use xCDAT, I noticed that most xCDAT APIs are based on "xarray.Dataset". Could it be possible to make it based on "xarray.DataArray" as well?

We spent a few months in the early development stage prototyping APIs with xarray.DataArray because it is similar to the cdms2.TransientVariable. However, we shifted to basing xCDAT APIs on xarray.Dataset because of limitations with the DataArray API (inability to hold bounds as a data variable). This is mentioned in the API Reference Overview here: https://xcdat.readthedocs.io/en/latest/api.html#overview

Here's helpful documentation on xarray.Dataset: https://docs.xarray.dev/en/stable/user-guide/data-structures.html#dataset

[mzelinka]

I just started playing around with xCDAT and had this same question. If I want to compute a spatial average on a quantity that is not a dataset (because of previous operations that have changed it from a dataset to a dataarray), how do I do it? Is there a way to assign the spatial coordinate bounds to a dataarray, based on its parent dataset?

[pochedls]

It seems like this isn't possible (at least not now). Considerable time went into figuring out how to do this, including adding accessors to modify the xarray dataarray structure to hold bounds. It didn't work (I think the issue was that user added accessors are not persistent). Others have raised this issue (e.g., here).

I think this is why xarray libraries sometimes pass the dataset around in function calls (even for a function that you think might be able to operate on just the dataarray, e.g., here)

[durack1]

@mzelinka did this example (of a subset region) work? https://xcdat.readthedocs.io/en/latest/examples/spatial-average.html#3.-Tropical-Region

[mzelinka]

yes, as long as I am working with a Dataset and not a DataArray it works fine.

[mzelinka]

I am having an issue reading in a netcdf file with xcdat, but xarray does not throw an error:

[mzelinka]

I think the issue is with the empty "height" coordinate....

[pochedls]

Sorry about this, Mark. I believe this issue and a bunch of other issues were fixed in version 0.3.1.

I think you can upgrade with conda install -c conda-forge xcdat=0.3.1

[mzelinka]

Thanks - that worked!

[mzelinka]

Ok, I have a new one (this may be related to my lack of understanding of xarray). I have a Dataset from one experiment called 'pert' for which I can effectively compute the spatial average. I have a Dataset from another control experiment called 'mocon' for which I can also effectively compute the spatial average. All good so far. When I difference these two Datasets to compute an anomaly Dataset called 'anom', it loses all the attributes. So I create a new Dataset 'anom2' that has the attributes from 'pert' assigned to 'anom'. But now when I compute a spatial average, only NaNs are returned. What am I doing wrong? See screenshots below:

[pochedls]

Thanks for catching this @mzelinka. I see why this is confusing.

anom = pert['tas'] - mocon is subtracting two datasets from one another. When you do this xarray performs a subtraction operation on all dataarrays, including lat_bnds and lon_bnds. These arrays thus become [[0., 0.][0., 0.], ..., [0., 0.]]. These are what are used to calculate the spatial weights (the weights are thus zero, which ends up resulting in NaN values in the spatial average operation).

A way to do this:

# get anomaly dataarray
daa = ds2.tas - ds1.tas
# convert to a dataset
dsa = daa.to_dataset()
# the bounds were dropped - add them back
# note that ds1.lat.bounds returns the name of the bounds for lat in ds1 (i.e., a string)
# note that ds1.bounds.get_bounds('Y') gets the latitude bounds object (i.e., a dataarray)
dsa[ds1.lat.bounds] = ds1.bounds.get_bounds('Y')
dsa[ds1.lon.bounds] = ds1.bounds.get_bounds('X')
# do spatial average
dsaa = dsa.spatial.average('tas')
# check values
dsaa.tas.values

array([0.50558825, 0.71024598, 0.64418223, ..., 6.66370024, 6.93816916,
7.04828337])

I wonder if it is worth adding an xcdat.dataset.to_dataset() function that ensures the bounds are included in the new dataset? Or another way to handle this kind of operation more easily.

[mzelinka]

Fascinating. Seems like xarray should know not to difference the data variables with '_bnds' in their name. Seems like strange default behavior.

[tomvothecoder]

Here's a simpler way to do this, which stays on the Dataset level to maintain bounds and avoid the additional steps for creating a new anomalies Dataset.

# 1. Calculate anomalies using the "tas" var from both datasets as new data variable, "tas_anom"
ds2["tas_anom"] = ds2.tas - ds1.tas

# 2. Perform spatial average on "tas_anom" 
ds2 = ds2.spatial.average("tas_anom")

Fascinating. Seems like xarray should know not to difference the data variables with '_bnds' in their name. Seems like strange default behavior.

In a Dataset, xarray typically stores bounds as data variables and not coordinate variables. Arithmetic at the Dataset level applies to all data variables, including bounds stored as data variables.

To avoid this behavior, you have to perform arithmetic on specific data variables within the Dataset.

# This performs arithmetic on all data variables, including bounds.
ds_new = ds2 - ds1

# This performs arithmetic only on `tas`.
ds_new["tas"] = ds2["tas"] - ds1["tas"]

[mzelinka]

I'm curious if the following behavior is expected / desired or if it is a bug. If load in a dataset, the data_var I want has the name 'sst'. But if I swap the longitudes using xcdat, this variable has a new name 'sea_surface_temperature'. I'm using xcdat v0.3.3. See screen shot below:

[mzelinka]

Looks fixed on my end too. Thanks!

[mzelinka]

Using xcdat v0.4.0, I just got an error when trying to to swap the longitude axes. I wonder if the issue is that coord_keys = get_dim_coords(ds, "X").coords.keys() on L236 of swap_lon_axis() returns not just the coordinate key for the "X" dimension but also "height" since this is a 2-meter air temperature.

Manually choosing only the longitude key on L260 resolves the error. That is: ds= ds.sortby(list(coord_keys)[0], ascending=True) does not throw this error.

See screen shots below.

[pochedls]

Just logging some notes...this doesn't seem to happen with all datasets (e.g., this dataset doesn't seem to have an issue/p/css03/esgf_publish/CMIP6/CMIP/NCAR/CESM2/amip/r1i1p1f1/Amon/tas/gn/v20190218/tas_Amon_CESM2_amip_r1i1p1f1_gn_195001-201412.nc).

Is this a one-off or happening to a lot of daily datasets?

Maybe this should be logged as a bug.

[mzelinka]

Thanks @pochedls . My example may be an edge case. I can confirm that for the file you listed, the get_dim_coords(ds, "X").coords.keys() command does not return height, just lon, as expected / desired (see screen shot).

[tomvothecoder]

Using xcdat v0.4.0, I just got an error when trying to to swap the longitude axes. I wonder if the issue is that coord_keys = get_dim_coords(ds, "X").coords.keys() on L236 of swap_lon_axis() returns not just the coordinate key for the "X" dimension but also "height" since this is a 2-meter air temperature.

Hey @mzelinka and @pochedls, this issue is related to height being a singleton coordinate. Singleton coordinate variables are included alongside other coordinate variables (refer to code example below). There is an open xarray issue regarding singletons here: pydata/xarray#6196

import xcdat as xc

fnames = "/p/css03/esgf_publish/CMIP6/CMIP/BCC/BCC-CSM2-MR/historical/r1i1p1f1/day/tasmax/gn/v20181126/*.nc"

ds = xc.open_mfdataset(fnames)

print(xc.get_dim_coords(ds, axis="X"))

# Notice how `height` is included in the `lon` coordinate var.

# <xarray.DataArray 'lon' (lon: 320)>
# array([  0.   ,   1.125,   2.25 , ..., 356.625, 357.75 , 358.875])
# Coordinates:
#   * lon      (lon) float64 0.0 1.125 2.25 3.375 4.5 ... 355.5 356.6 357.8 358.9
#     height   float64 2.0
# Attributes:
#     bounds:         lon_bnds
#     units:          degrees_east
#     axis:           X
#     long_name:      Longitude
#     standard_name:  longitude

How it affects swap_lon_axis()

In swap_lon_axis(), we are getting all of the coordinates keys using lon, which includes height.
This causes swap_lon_axis(ds, sort_ascending=True) to break on this file.

xcdat/xcdat/axis.py

Lines 235 to 236 in d5d4fd7

	ds = dataset.copy()
	coord_keys = get_dim_coords(ds, "X").coords.keys()

xcdat/xcdat/axis.py

Lines 259 to 260 in d5d4fd7

	if sort_ascending:
	ds = ds.sortby(list(coord_keys), ascending=True)

Temporary workaround:

• Set sort_ascending=False in swap_lon_axis() and perform sorting manually

The fix for swap_lon_axis() is to ignore singleton coordinate keys.

I will open a new issue to fix this.

[tomvothecoder]

I am going to lock this thread because v0.3.0 is outdated. Please refer to #391 for v0.4.0 testing feedback.