Latest version drops vital coordinates of CMIP6 output #305

jbusecke · 2020-12-09T22:44:28Z

Description

I have run into trouble with the latest intake-esm version and the google CMIP6 data. When loading datasets via intake esm

import intake
import fsspec
import xarray as xr
from cmip6_preprocessing.preprocessing import combined_preprocessing

col = intake.open_esm_datastore(
        "https://raw.githubusercontent.com/NCAR/intake-esm-datastore/master/catalogs/pangeo-cmip6.json"
    )
cat = col().search(
        source_id='CanESM5',
        experiment_id='historical',
        variable_id='thetao',
        member_id="r1i1p1f1",
        table_id="Omon",
        grid_label='gn',
    )
cat

zarr_kwargs = {'consolidated':True}
ddict = cat.to_dataset_dict(
            zarr_kwargs=zarr_kwargs,
            preprocess=combined_preprocessing,
            storage_options={"token": "anon"},
        )
_, ds = ddict.popitem()
ds

This is missing a bunch of the coordinates, which breaks many of my workflows (and testing) using cmip6 data.

I assume this is an intake-esm issue since I can load and preprocess the file manually

# now load the 'raw' zarr file and preprocess in the same way
mm = fsspec.get_mapper(
            cat.df["zstore"][0]
        )
ds_raw = xr.open_zarr(mm, **zarr_kwargs)
ds_manual = combined_preprocessing(ds_raw)
ds_manual

and the coordinates are present (e.g. lon_bounds/lat_bounds).

This only happens with the current version. I tested with v2020.08.15 and it does not happen.

Version information: output of `intake_esm.show_versions()`

Paste the output of intake_esm.show_versions() here:

INSTALLED VERSIONS
------------------

cftime: 1.2.1
dask: 2.30.0
fastprogress: 0.2.7
fsspec: 0.8.4
gcsfs: 0.7.1
intake: 0.6.0
intake_esm: 2020.11.4
netCDF4: 1.5.4
pandas: 1.1.4
requests: 2.24.0
s3fs: 0.4.2
xarray: 0.16.1
zarr: 2.5.0

The text was updated successfully, but these errors were encountered:

andersy005 · 2020-12-10T19:59:16Z

@jbusecke, thank you for opening this issue! The only major change that may have introduced this issue is #287. I'm going to look into this and will get back to you ASAP.

andersy005 · 2020-12-10T20:35:20Z

@jbusecke, I'm having trouble reproducing the issue when I use the code snippet you provided. I don't see lon_bounds and lat_bounds in the raw zarr store:

# I extracted the store path via cat.df["zstore"][0]
In [17]: store = 'gs://cmip6/CMIP/CCCma/CanESM5/historical/r1i1p1f1/Omon/thetao/gn/'

In [18]: mm = fsspec.get_mapper(store)

In [19]: ds_raw = xr.open_zarr(mm, **zarr_kwargs)

In [20]: ds_raw
Out[20]: 
<xarray.Dataset>
Dimensions:             (bnds: 2, i: 360, j: 291, lev: 45, time: 1980, vertices: 4)
Coordinates:
  * i                   (i) int32 0 1 2 3 4 5 6 ... 353 354 355 356 357 358 359
  * j                   (j) int32 0 1 2 3 4 5 6 ... 284 285 286 287 288 289 290
    latitude            (j, i) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
  * lev                 (lev) float64 3.047 9.454 16.36 ... 5.375e+03 5.625e+03
    lev_bnds            (lev, bnds) float64 dask.array<chunksize=(45, 2), meta=np.ndarray>
    longitude           (j, i) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
  * time                (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:0...
    time_bnds           (time, bnds) object dask.array<chunksize=(1980, 2), meta=np.ndarray>
Dimensions without coordinates: bnds, vertices
Data variables:
    thetao              (time, lev, j, i) float32 dask.array<chunksize=(6, 45, 291, 360), meta=np.ndarray>
    vertices_latitude   (j, i, vertices) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>
    vertices_longitude  (j, i, vertices) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>

In [21]: ds_manual = combined_preprocessing(ds_raw)
/Users/abanihi/opt/miniconda3/envs/intake-esm-dev/lib/python3.8/site-packages/xarray/core/indexing.py:1361: PerformanceWarning: Slicing is producing a large chunk. To accept the large
chunk and silence this warning, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    ...     array[indexer]

To avoid creating the large chunks, set the option
    >>> with dask.config.set(**{'array.slicing.split_large_chunks': True}):
    ...     array[indexer]
  return self.array[key]

In [22]: ds_manual
Out[22]: 
<xarray.Dataset>
Dimensions:             (bnds: 2, lev: 45, time: 1980, vertex: 4, x: 360, y: 291)
Coordinates:
  * x                   (x) float32 0.5 1.5 2.5 3.5 ... 356.5 357.5 358.5 359.5
  * y                   (y) float32 -78.3935 -78.19057 ... 71.618454 71.65057
    lat                 (y, x) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
  * lev                 (lev) float64 3.047 9.454 16.36 ... 5.375e+03 5.625e+03
    lev_bounds          (lev, bnds) float64 dask.array<chunksize=(45, 2), meta=np.ndarray>
    lon                 (y, x) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
  * time                (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:0...
    time_bounds         (time, bnds) object dask.array<chunksize=(1980, 2), meta=np.ndarray>
    vertices_latitude   (y, x, vertex) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>
    vertices_longitude  (y, x, vertex) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>
  * bnds                (bnds) int64 0 1
  * vertex              (vertex) int64 0 1 2 3
Data variables:
    thetao              (time, lev, y, x) float32 dask.array<chunksize=(6, 45, 291, 360), meta=np.ndarray>

Am I missing anything? Is there any other information you can provide?

jbusecke · 2020-12-10T20:37:57Z

This looks like a problem with the combined_preprocessing function. What's the version of cmip6_preprocessing that you are using?

andersy005 · 2020-12-10T20:41:43Z

I'm running:

In [26]: cmip6_preprocessing.__version__
Out[26]: 'v0.1.4'

Should I upgrade?

andersy005 · 2020-12-10T20:57:09Z

Should I upgrade?

Looks like 'v0.1.4' is the latest release

jbusecke · 2020-12-10T23:04:57Z

I think I was using the upstream master, could you try with that? If that fixes this issue Ill release soon.

jbusecke · 2020-12-11T03:32:46Z

I will have a more detailed overview about this issue once cmip6_preprocessing/62 is up and running. This is where I discovered the issue, for now I have pinned the older intake-esm versions.

jbusecke · 2020-12-18T01:35:02Z

cmip6_preprocessing/67 shows the issue. Just by updating intake-esm I get a bunch of failing tests due to missing coordinates. Any idea what could cause this?

andersy005 · 2020-12-18T01:55:53Z

@jbusecke,

Your assessment is dead on. When I turn the aggregations off, I can see that there are coordinate information that get thrown away when the aggregations are enabled:

Aggreations OFF:

 Dimensions:             (bnds: 2, i: 360, j: 291, lev: 45, time: 1980, vertices: 4)
 Coordinates:
   * i                   (i) int32 0 1 2 3 4 5 6 ... 353 354 355 356 357 358 359
   * j                   (j) int32 0 1 2 3 4 5 6 ... 284 285 286 287 288 289 290
     latitude            (j, i) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
   * lev                 (lev) float64 3.047 9.454 16.36 ... 5.375e+03 5.625e+03
     lev_bnds            (lev, bnds) float64 dask.array<chunksize=(45, 2), meta=np.ndarray>
     longitude           (j, i) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
   * time                (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:0...
     time_bnds           (time, bnds) object dask.array<chunksize=(1980, 2), meta=np.ndarray>
 Dimensions without coordinates: bnds, vertices
 Data variables:
     thetao              (time, lev, j, i) float32 dask.array<chunksize=(6, 45, 291, 360), meta=np.ndarray>
     vertices_latitude   (j, i, vertices) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>
     vertices_longitude  (j, i, vertices) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>

Aggregations ON:

 Dimensions:    (i: 360, j: 291, lev: 45, member_id: 1, time: 1980)
 Coordinates:
   * i          (i) int32 0 1 2 3 4 5 6 7 8 ... 352 353 354 355 356 357 358 359
   * j          (j) int32 0 1 2 3 4 5 6 7 8 ... 283 284 285 286 287 288 289 290
     latitude   (j, i) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
   * lev        (lev) float64 3.047 9.454 16.36 ... 5.126e+03 5.375e+03 5.625e+03
     longitude  (j, i) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
   * time       (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
   * member_id  (member_id) <U8 'r1i1p1f1'
 Data variables:
     thetao     (member_id, time, lev, j, i) float32 dask.array<chunksize=(1, 6, 45, 291, 360), meta=np.ndarray>

I'm going to look into this and will ping you when I have a patch for this bug...

andersy005 · 2020-12-18T04:50:58Z

@jbusecke, thank you for the minimal example too! I can confirm that #306 seems to have fixed this issue:

In [1]: import intake

In [2]: import xarray as xr

In [3]: import fsspec

In [4]: from cmip6_preprocessing.preprocessing import combined_preprocessing

In [5]: col = intake.open_esm_datastore(
   ...:         "https://raw.githubusercontent.com/NCAR/intake-esm-datastore/master/catalogs/pangeo-cmip6.jso
   ...: n"
   ...:     )

In [6]: cat = col.search(
   ...:         source_id='CanESM5',
   ...:         experiment_id='historical',
   ...:         variable_id='thetao',
   ...:         member_id="r1i1p1f1",
   ...:         table_id="Omon",
   ...:         grid_label='gn',
   ...:     )

In [7]: cat
Out[7]: <pangeo-cmip6 catalog with 1 dataset(s) from 1 asset(s)>

In [8]: zarr_kwargs = {'consolidated':True}

intake-esm

In [9]: ddict = cat.to_dataset_dict(
   ...:             zarr_kwargs=zarr_kwargs,
   ...:             preprocess=combined_preprocessing,
   ...:             storage_options={"token": "anon"},
   ...:         )
   ...: _, ds = ddict.popitem()
   ...: ds

--> The keys in the returned dictionary of datasets are constructed as follows:
        'activity_id.institution_id.source_id.experiment_id.table_id.grid_label'

Out[9]: ███████████████████████████████████████████| 100.00% [1/1 00:00<00:00]
<xarray.Dataset>
Dimensions:        (bnds: 2, lev: 45, member_id: 1, time: 1980, vertex: 4, x: 360, y: 291)
Coordinates:
  * x              (x) float32 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
  * y              (y) float32 -78.3935 -78.19058 ... 89.36653 89.74177
    lat            (y, x) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
  * lev            (lev) float64 3.047 9.454 16.36 ... 5.375e+03 5.625e+03
    lev_bounds     (lev, bnds) float64 dask.array<chunksize=(45, 2), meta=np.ndarray>
    lon            (y, x) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
  * time           (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
    time_bounds    (time, bnds) object dask.array<chunksize=(1980, 2), meta=np.ndarray>
    lat_verticies  (y, x, vertex) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>
    lon_verticies  (y, x, vertex) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>
  * bnds           (bnds) int64 0 1
  * vertex         (vertex) int64 0 1 2 3
    lon_bounds     (bnds, y, x) float32 dask.array<chunksize=(1, 291, 360), meta=np.ndarray>
    lat_bounds     (bnds, y, x) float32 dask.array<chunksize=(1, 291, 360), meta=np.ndarray>
  * member_id      (member_id) <U8 'r1i1p1f1'
Data variables:
    thetao         (member_id, time, lev, y, x) float32 dask.array<chunksize=(1, 6, 45, 291, 360), meta=np.ndarray>

xarray + fsspec

In [10]: mm = fsspec.get_mapper(
    ...:             cat.df["zstore"][0]
    ...:         )
    ...: ds_raw = xr.open_zarr(mm, **zarr_kwargs)
    ...: ds_manual = combined_preprocessing(ds_raw)
    ...: ds_manual
    ...: 

Out[10]: 
<xarray.Dataset>
Dimensions:        (bnds: 2, lev: 45, time: 1980, vertex: 4, x: 360, y: 291)
Coordinates:
  * x              (x) float32 0.5 1.5 2.5 3.5 4.5 ... 356.5 357.5 358.5 359.5
  * y              (y) float32 -78.3935 -78.19058 ... 89.36653 89.74177
    lat            (y, x) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
  * lev            (lev) float64 3.047 9.454 16.36 ... 5.375e+03 5.625e+03
    lev_bounds     (lev, bnds) float64 dask.array<chunksize=(45, 2), meta=np.ndarray>
    lon            (y, x) float32 dask.array<chunksize=(291, 360), meta=np.ndarray>
  * time           (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
    time_bounds    (time, bnds) object dask.array<chunksize=(1980, 2), meta=np.ndarray>
    lat_verticies  (y, x, vertex) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>
    lon_verticies  (y, x, vertex) float32 dask.array<chunksize=(291, 360, 4), meta=np.ndarray>
  * bnds           (bnds) int64 0 1
  * vertex         (vertex) int64 0 1 2 3
    lon_bounds     (bnds, y, x) float32 dask.array<chunksize=(1, 291, 360), meta=np.ndarray>
    lat_bounds     (bnds, y, x) float32 dask.array<chunksize=(1, 291, 360), meta=np.ndarray>
Data variables:
    thetao         (time, lev, y, x) float32 dask.array<chunksize=(6, 45, 291, 360), meta=np.ndarray>

andersy005 · 2020-12-18T07:29:03Z

Also, it appears that with #306, everything is back to working (including jbusecke/xMIP#67) with an exception of one test

-- Docs: https://docs.pytest.org/en/stable/warnings.html
========================== short test summary info ===========================
FAILED cmip6_preprocessing/tests/test_preprocessing_cloud.py::test_check_dim_coord_values[CESM2-thetao-ssp585-gn]
= 1 failed, 644 passed, 2101 skipped, 38 xfailed, 9030 warnings in 8754.16s (2:25:54) =

P.S. That's plenty of long running tests... I wasn't sure how long they were going to take to finish running 😅

jbusecke · 2020-12-18T15:36:16Z

Yeah it's not great to run for 1 hour...I will try to improve this in the next days, but I wanted to get this running first.

Thanks for the fix! I am rerunning the tests, and will investigate what is happening with that one model...

jbusecke · 2020-12-18T15:36:44Z

Also thanks for the quick fix!

andersy005 · 2020-12-18T18:01:45Z

Of course! Just merged #306 into master. When you get a chance, can you test jbusecke/xMIP#67 against intake-esm's master branch and let me know how it goes? I plan to release a new version once I have your confirmation that everything is working.

jbusecke · 2020-12-18T18:29:03Z

Just switched the CI to the upstream master. Now we just need to wait an hour...🤪

jbusecke · 2020-12-18T22:11:27Z

It seems that everything is back to normal with the newest master (jbusecke/xMIP#67). Your single model failure could have been caused by some hickup with the cloud data?

andersy005 · 2020-12-18T22:22:00Z

Your single model failure could have been caused by some hickup with the cloud data?

Perhaps... Thank you for the update! I am going to cut a new release by the end of the day today.

andersy005 · 2020-12-19T01:27:40Z

Closing this as it has been addressed in #306

Please reopen if you still encounter this issue with the latest stable (v2020.12.18) version

andersy005 mentioned this issue Dec 18, 2020

🐛 Bug Fix: Disable _requested_variables for single variable assets #306

Merged

3 tasks

andersy005 added the bug Issues that present a reasonable conviction there is a reproducible bug. label Dec 18, 2020

andersy005 closed this as completed Dec 19, 2020

jbusecke mentioned this issue Dec 21, 2020

Activate single rerun with 10s delay for flakey tests. jbusecke/xMIP#68

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Latest version drops vital coordinates of CMIP6 output #305

Latest version drops vital coordinates of CMIP6 output #305

jbusecke commented Dec 9, 2020

andersy005 commented Dec 10, 2020

andersy005 commented Dec 10, 2020

jbusecke commented Dec 10, 2020

andersy005 commented Dec 10, 2020

andersy005 commented Dec 10, 2020

jbusecke commented Dec 10, 2020

jbusecke commented Dec 11, 2020 •

edited

Loading

jbusecke commented Dec 18, 2020

andersy005 commented Dec 18, 2020 •

edited

Loading

andersy005 commented Dec 18, 2020

andersy005 commented Dec 18, 2020

jbusecke commented Dec 18, 2020

jbusecke commented Dec 18, 2020

andersy005 commented Dec 18, 2020

jbusecke commented Dec 18, 2020

jbusecke commented Dec 18, 2020

andersy005 commented Dec 18, 2020

andersy005 commented Dec 19, 2020

Latest version drops vital coordinates of CMIP6 output #305

Latest version drops vital coordinates of CMIP6 output #305

Comments

jbusecke commented Dec 9, 2020

Description

Version information: output of intake_esm.show_versions()

andersy005 commented Dec 10, 2020

andersy005 commented Dec 10, 2020

jbusecke commented Dec 10, 2020

andersy005 commented Dec 10, 2020

andersy005 commented Dec 10, 2020

jbusecke commented Dec 10, 2020

jbusecke commented Dec 11, 2020 • edited Loading

jbusecke commented Dec 18, 2020

andersy005 commented Dec 18, 2020 • edited Loading

andersy005 commented Dec 18, 2020

andersy005 commented Dec 18, 2020

jbusecke commented Dec 18, 2020

jbusecke commented Dec 18, 2020

andersy005 commented Dec 18, 2020

jbusecke commented Dec 18, 2020

jbusecke commented Dec 18, 2020

andersy005 commented Dec 18, 2020

andersy005 commented Dec 19, 2020

Version information: output of `intake_esm.show_versions()`

jbusecke commented Dec 11, 2020 •

edited

Loading

andersy005 commented Dec 18, 2020 •

edited

Loading