# Xarray and collections of forecasts

**URL:** <https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054>\
**Category:** Science\
**Created:** [January 6, 2023, 2:55am UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054 "2023-01-06T02:55:40Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![abkfenris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/abkfenris/32/507_2.png) [@abkfenris](https://discourse.pangeo.io/u/abkfenris)\
**Post date:** [January 6, 2023, 2:55am UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/1 "2023-01-06T02:55:41Z")

</div>

For a while I’ve been pondering how to represent collections of forecasts in Xarray, and I know others have been too.

While THREDDS has a nice way to assemble multiple model runs as Forecast Model Run Collections, the staggered nature of the coordinates makes them messy to fit into a single `xr.Dataset`.

 ![Forecast Model Run Collections](https://canada1.discourse-cdn.com/flex030/uploads/pangeo/original/2X/a/a083c34b76effa485b0429cec9ff9dae85a05083.png)

After some conversation [during scheeming around Xpublish and ZarrDAP](https://github.com/xarray-contrib/xpublish/issues/138) (I’ll get back to working on evolving Xpublish soon!), it’s been even more on my mind to the point I couldn’t sleep one night until I wrote a [mock up](https://gist.github.com/abkfenris/ebab2db91fdc91c14084380334767273) of an API (and sent it to @rsignell and others to try to make it someone else’s problem).

What kept me up till I got a mock up made was figuring out how to take advantage of @TomNicholas work with [Datatrees](https://xarray-datatree.readthedocs.io/en/latest/index.html).

```plaintext
dt = xarray_fmrc.from_model_runs([ds0, ds1])
dt
DataTree('None', parent=None)
│ Dimensions: (forecast_reference_time: 2,
│ constant_forecast: 242, constant_offset: 121)
│ Coordinates:
│ * forecast_reference_time (forecast_reference_time) datetime64[ns] 2022-12...
│ * constant_forecast (constant_forecast) datetime64[ns] 2022-12-02 .....
│ * constant_offset (constant_offset) timedelta64[ns] 06:00:00 ... 5...
│ Data variables:
│ model_run_path (forecast_reference_time) <U29 'model_run/2022-1...
└── DataTree('model_run')
    ├── DataTree('2022-12-01T18:00:00')
    │ Dimensions: (forecast_reference_time: 1, time: 121,
    │ latitude: 220, longitude: 215)
    │ Coordinates:
    │ * longitude (longitude) float64 -79.95 -79.86 ... -60.13 -60.04
    │ * latitude (latitude) float64 27.03 27.12 ... 47.32 47.41
    │ * time (time) datetime64[ns] 2022-12-02 ... 2022-12-07
    │ * forecast_reference_time (forecast_reference_time) datetime64[ns] 2022-12...
    │ Data variables:
    │ wind_speed (forecast_reference_time, time, latitude, longitude) float32 ...
    │ wind_from_direction (forecast_reference_time, time, latitude, longitude) float32 ...
    │ Attributes: (12/178)
    │ ...
    └── DataTree('2022-12-12T18:00:00')
            Dimensions: (forecast_reference_time: 1, time: 121,
                                          latitude: 220, longitude: 215)
            Coordinates:
              * longitude (longitude) float64 -79.95 -79.86 ... -60.13 -60.04
              * latitude (latitude) float64 27.03 27.12 ... 47.32 47.41
              * time (time) datetime64[ns] 2022-12-13 ... 2022-12-18
              * forecast_reference_time (forecast_reference_time) datetime64[ns] 2022-12...
            Data variables:
                wind_speed (forecast_reference_time, time, latitude, longitude) float32 ...
                wind_from_direction (forecast_reference_time, time, latitude, longitude) float32 ...
            Attributes: (12/178)
                ...

```

I think that by placing each model run in the datatree, then accessor methods could be used for the various views that users want from a collection of forecasts. Thus methods along the lines of `dt.fmrc.constant_offset("12H")` to get values that are 12 hours from the time that a forecast was generated or `dt.fmrc.best()` to give a dataset with the best forecast data (least amount of time between generation and a forecasted time).

In trying to go from a quickly slapped together mock up to a post with slightly easier to understand example I ended up writing enough code tinkering that the mock up has started to become a reality: [xarray\_fmrc](https://github.com/abkfenris/xarray_fmrc).

I know there has also been some exploration of using [`kerchunk.subchunk()`](https://fsspec.github.io/kerchunk/reference.html#kerchunk.utils.subchunk) and other ways to represent forecast collections, so I’d love to hear other thoughts.

---

<div class="post-metadata">

**Author:** ![rabernat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/rabernat/32/22_2.png) [@rabernat](https://discourse.pangeo.io/u/rabernat)\
**Post date:** [January 6, 2023, 4:00am UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/2 "2023-01-06T04:00:10Z")

</div>

This is very cool @abkfenris! Our community has wrestled with the question of how to best to organize forecast data for a long time…but now we have Datatree! 🙏 📈🌲 I think you’re heading in the right direction.

If you haven’t seen this Pangeo Forge thread, it might have some useful insights from @chiaral

> <https://github.com/pangeo-forge/staged-recipes/issues/17>
>
> \<!--
> This template is to describe a potential pipeline for Pangeo Forge to crea…te analysis-ready, cloud-optimized data from an upstream data repository.
> 
> A pipeline has three basic stages:
> 1. Download the source files from the upstream repository in whatever format they are stored.
> 2. Perform any transformations that are needed in order to make the data "analysis ready."
> 3. Write out a new dataset in a cloud optimized format
> \--\>
> 
> \## Source Dataset
> 
> 
> These are “reforecasts” of the new GEFSv12 system. They are retrospective forecasts spanning the period 2000-2019.
> These reforecasts are not as numerous as the real-time data; they were generated only once
> per day, from 00 UTC initial conditions, and only 5 members were provided, with the following
> exception. Once weekly, an 11-member reforecast was generated, and these extend in lead
> time to +35 days. 
> 
> 
> \- Link to \[documentation\](https://noaa-gefs-retrospective.s3.amazonaws.com/Description\_of\_reforecast\_data.pdf)
> \- Link to \[data AWS bucket\](https://noaa-gefs-retrospective.s3.amazonaws.com/index.html#GEFSv12/reforecast/)
> \- The file format is \*grib2\* - both the original and the cloud copy.
> \- The files have multiple dimensions because the are forecast data (besides times/x/y/z they have lead time, ensemble members), they are organized as (from the PDF file linked above):
> - the directory tree structure under GEFSv12/reforecast/ is by year;
> - there are separate subdirectories for each yyyymmddhh, thus 2000010100 to 2000123100 for the year 2000. 
> - Under each yyyymmddhh subdirectory, there are subdirectories c00, p01, p02, p03, p04 for the five individual member forecasts. Once per week, 11 reforecast members were computed, and the directories for those days extend
> through p10.
> - Individual grib files have file names such as “variable\_yyyymmddhh\_member.grib2”
> \- Currently I played with loading some part of the files (CONUS precip/temp) in the AWS pangeo deployment using xarray using the rasterio engine and concatenating them, but it's veeeeeery slow and I loose info from the grib files because I don't use a grib engine:
> \`\`\`
> import s3fs
> import xarray as xr
> s3 = s3fs.S3FileSystem(anon=True)
> 
> def preprocessing\_function\_ptc(path\_loop, im):
> prova = \['/'+ifg for ifg in path\_loop.split('/')\[1:\]\]
> ds = xr.open\_rasterio('https://noaa-gefs-retrospective.s3.amazonaws.com'+''.join(prova),
> chunks={'x':'200MB', 'band':-1})
> ds = ds.sel(y = slice(51,20), x = slice(229, 302))
> start\_time = pd.to\_datetime(str(path\_loop.split('/')\[4\]),format='%Y%m%d%H')
> ds.coords\['time'\] = start\_time
> ds = ds.expand\_dims('time')
> ds.coords\['member\_id'\] = im
> ds = ds.expand\_dims('member\_id')
> return ds
> \`\`\`
> \`\`\`
> path\_b = 's3://noaa-gefs-retrospective/GEFSv12/reforecast/\*'
> lglob = s3fs.S3FileSystem.glob(s3,path=path\_b)
> for ilyear in lglob\[0:4\]:
> print(ilyear)
> lglob2 = s3fs.S3FileSystem.glob(s3,path='s3://'+ilyear+'/\*')
> ds\_p\_all = \[\]
> ds\_t\_all = \[\]
>     
> for ilyear2 in lglob2\[0:4\]:
> print(ilyear2)
> ds\_p\_all1 = \[\]
> ds\_t\_all1 = \[\]
> for member in \['c00','p01','p02','p03','p04'\]:
> lglob3 = s3fs.S3FileSystem.glob(s3,path='s3://'+ilyear2+'/'+member+'/Days:1-10/\*')
> for ilg3 in lglob3:
> if 'apcp' in ilg3:
> ds\_p = preprocessing\_function\_ptc(ilg3, member)
> ds\_p\_all1.append(ds\_p)
> elif 'tmp\_2m' in ilg3:
> ds\_t = preprocessing\_function\_ptc(ilg3, member)
> ds\_t\_all1.append(ds\_t)                    
> temp\_p = xr.concat(ds\_p\_all1, dim='member\_id')
> temp\_t = xr.concat(ds\_t\_all1, dim='member\_id')
> ds\_p\_all.append(temp\_p)
> ds\_t\_all.append(temp\_t)
> precip = xr.concat(ds\_p\_all, dim='time')
> print(precip)
> tmp\_2m = xr.concat(ds\_t\_all, dim='time')
> \`\`\`
> \- There is no password.
> 
> \## Transformation / Alignment / Merging
> 
> How to combine forecast data in a zarr format is not entirely clear to me yet - meaning what type of chunking should be used or merging. Often times analysis are carried out both as a function of start times and lead times, so I don't think we can have one structure that makes everyone happy. But I will figure out what is the structure that makes more sense.
> \<!--
> Describe below how the files should be combined into one analysis-ready dataset.
> For example, "the files should be concatenated along the time dimension."
> Are there any other transformations or checks that should be performed to make the data more "analysis ready"?
> \--\>
> 
> 
> \## Output Dataset
> 
> I am interested in transforming them in zarr format. It might be necessary to first download them and them opening them using a grib engine and then push them to the cloud, because as of now I couldn't use a grib engine.
> \<!--
> How do you want the output of the pipeline to be stored?
> Cloud optimized formats such as zarr, tiledb, or parquet are recommended.
> If possible, provide details on how you would like the output to be structured
> (e.g. number of different output datasets, chunk / partition size, etc.)
> \--\>

---

<div class="post-metadata">

**Author:** ![TomNicholas](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/tomnicholas/32/3189_2.png) [@TomNicholas](https://discourse.pangeo.io/u/TomNicholas)\
**Post date:** [January 6, 2023, 4:44pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/3 "2023-01-06T16:44:44Z")

</div>

This is great to see @abkfenris! I would love to help out however I can.

If you have any thoughts on what kind of API improvements on datatree would help you then I’m all ears. There is [some discussion on API methods here](https://github.com/xarray-contrib/datatree/issues/79). I actually added a `.filter` method [just now](https://github.com/xarray-contrib/datatree/pull/185) which might be useful.

---

<div class="post-metadata">

**Author:** ![abkfenris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/abkfenris/32/507_2.png) [@abkfenris](https://discourse.pangeo.io/u/abkfenris)\
**Post date:** [January 8, 2023, 7:22pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/4 "2023-01-08T19:22:05Z")

</div>

I’ve gotten the accessor methods working at least for my few test cases, so `dt.fmrc.best()`, `dt.fmcr.constant_offset("6h")`, and `dt.fmrc.constant_forecast("2023-01-10")` all work now.

@TomNicholas I’d definitely love some help. I’m somewhat just using datatree as a pre-existing structure to jam model runs into and somewhat naively looping over them.

While filter might take care of this, it would also be nice to be able to filter nodes with broadcasted dataset methods.

Maybe an extra argument along the lines of `dt.sel(time=dt, _filter=True)` to catch `KeyError`s and instead only return a subset of the tree with the selection.

@rabernat sometimes it seems like the only commonality amongst forecasts is that you can at least expect that all the data generated at the same time will be in the same place and be able to be used together.

That’s why I’m treating each one as it’s own group when storing in the tree. Only when you try to use a specific accessor method does it actually go in and try to get them to play together. At that point you should have enough information to filter out the unwanted duplicate data that would otherwise make it ugly to jam into a dataset.

---

<div class="post-metadata">

**Author:** ![TomNicholas](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/tomnicholas/32/3189_2.png) [@TomNicholas](https://discourse.pangeo.io/u/TomNicholas)\
**Post date:** [January 9, 2023, 7:25am UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/5 "2023-01-09T07:25:02Z")

</div>

> [@abkfenris](#):
>
> I’ve gotten the accessor methods working at least for my few test cases, so `dt.fmrc.best()`, `dt.fmcr.constant_offset("6h")`, and `dt.fmrc.constant_forecast("2023-01-10")` all work now.

Amazing!

> [@abkfenris](#):
>
> Maybe an extra argument along the lines of `dt.sel(time=dt, _filter=True)` to catch `KeyError`s and instead only return a subset of the tree with the selection.

I can see how that would be useful. If I’m understanding correctly this is quite a similar suggestion to in [this datatree issue](https://github.com/xarray-contrib/datatree/issues/67), just you’re talking about ignoring errors from missing variables instead of ignoring errors from missing dimensions. ANother idea would be to add something like `dt.sel(time=dt, errors="ignore")` inspired by the option xarray’s `.drop_vars` already has.

@abkfenris in general any snippets of “I would like it if this code acted like this” is really helpful to me.

---

<div class="post-metadata">

**Author:** ![abkfenris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/abkfenris/32/507_2.png) [@abkfenris](https://discourse.pangeo.io/u/abkfenris)\
**Post date:** [January 9, 2023, 3:12pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/6 "2023-01-09T15:12:02Z")

</div>

I thought I had seen an issue along those lines, but I think I searched ‘coordinates’ and ‘exception’ (as technically I’m selecting a non-dimension coordinate in one case) so I didn’t find that one. `dt.sel(time=dt, errors="ignore")` would work great. I’ve got some more ideas that I’ll continue in the issue.

A method like `.concat_child_datasets` would also be helpful. Between those two it could simplify two of my methods from.

```plaintext
    def constant_offset(self, offset: Union[str, int, float, timedelta]) -> xr.Dataset:
        timedelta = pd.to_timedelta(offset)

        filtered_ds = []

        for child in self.datatree_obj["model_run"].children.values():
            try:
                selected = child.ds.sel(forecast_offset=timedelta)
                filtered_ds.append(selected)
            except KeyError:
                pass

        combined = xr.concat(filtered_ds, "time")
        combined = combined.sortby("time")
        return combined

```

to

```plaintext
    def constant_offset(self, offset: Union[str, int, float, timedelta]) -> xr.Dataset:
        return self.datatree_obj["model_run"].sel(forecast_offset=offset).concat_child_datasets("time").sortby("time")

```

I think it could allow other datatrees with similarly really closely related but non-alignable datasets to be quickly collapsed into aligned datasets with less boilerplate.

---

<div class="post-metadata">

**Author:** ![rabernat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/rabernat/32/22_2.png) [@rabernat](https://discourse.pangeo.io/u/rabernat)\
**Post date:** [January 9, 2023, 3:43pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/7 "2023-01-09T15:43:36Z")

</div>

`.concat_child_datasets` is analogous to the `.merge_child_datasets` idea we discussed at the LEAP workshop last week. It seems like collapsing datatrees into datasets could be a broadly useful feature.

---

<div class="post-metadata">

**Author:** ![TomNicholas](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/tomnicholas/32/3189_2.png) [@TomNicholas](https://discourse.pangeo.io/u/TomNicholas)\
**Post date:** [January 9, 2023, 4:07pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/8 "2023-01-09T16:07:42Z")

</div>

`concat_child_datasets` could be something like

```python
def concat_child_datasets(self, **concat_options):
    collapsed_ds = xr.concat([node.ds for node in self.subtree], **concat_options)

    self.ds = collapsed_ds
    
    # then delete / orphan all the children of this node

    return self

```

if you’re interested in submitted a PR?.. Else I’ll try to get to it soon.

---

<div class="post-metadata">

**Author:** ![abkfenris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/abkfenris/32/507_2.png) [@abkfenris](https://discourse.pangeo.io/u/abkfenris)\
**Post date:** [January 9, 2023, 9:30pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/9 "2023-01-09T21:30:44Z")

</div>

I was thinking that `concat_child_datasets` would return a dataset rather than a datatree.

How about `concat_children()` and `merge_children()` that return single nodes as then datasets are also directly accessible from `.ds`?

---

<div class="post-metadata">

**Author:** ![kdl0013](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/kdl0013/32/2023_2.png) [@kdl0013](https://discourse.pangeo.io/u/kdl0013)\
**Post date:** [January 9, 2023, 10:16pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/10 "2023-01-09T22:16:36Z")

</div>

You mentioned forecast data and python has a community for that [climpred: verification of weather and climate forecasts — climpred documentation](https://climpred.readthedocs.io/en/stable/) called climpred. It could be possible to collaborate with them to see how forecast data is currently processed and they might have an idea.

---

<div class="post-metadata">

**Author:** ![TomNicholas](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/tomnicholas/32/3189_2.png) [@TomNicholas](https://discourse.pangeo.io/u/TomNicholas)\
**Post date:** [January 9, 2023, 11:13pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/11 "2023-01-09T23:13:43Z")

</div>

> [@abkfenris](#):
>
> How about `concat_children()` and `merge_children()` that return single nodes as then datasets are also directly accessible from `.ds`?

That’s what I was thinking too, I just didn’t express it very clearly. [I’ve raised an issue to continue this discussion](https://github.com/xarray-contrib/datatree/issues/192)

---

<div class="post-metadata">

**Author:** ![abkfenris](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/abkfenris/32/507_2.png) [@abkfenris](https://discourse.pangeo.io/u/abkfenris)\
**Post date:** [January 11, 2023, 4:00pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/12 "2023-01-11T16:00:25Z")

</div>

Hi @kdl0013 thanks for the reminder about climpred! I had forgotten about that project and it was interesting to explore what’s happened since I last checked it out.

Climpred has made things work by reshaping forecasts into `init` (`forecast_reference_time`, the time the forecast was ‘made’) and `lead` (integer offset with units encoded as an attribute) dimensions so that they can be aggregated into a single dataset.

So far, I’ve been largely approaching this problem with my [NERACOOS](https://www.neracoos.org/) data manager hat on. Most of our users (fishermen, Coast Guard, local communities) really just want the best available data and aren’t worried about individual model runs, but I’ve also got a handful of scientists who are interested in those. We are currently running ERDDAP and THREDDS servers, but don’t have full THREDDS FMRC serving set up (instead I’ve got something manually hacked together with separate datasets for each model run).

For us (and I’m guessing a lot of other orgs), it probably doesn’t make too much sense to change our model storage to be `init` x `lead` shaped with our normal usage patterns.

Since there is enough information to get `init` x `lead` (even with some changes I’m pondering below), I do think it would be reasonable to provide a method that could return a tree as a climpred compatible dataset.

Right now I’m capturing `forecast_reference_time` (or `init`) as a dimension and `forecast_period` (or `lead`) as a non-dimension coordinate (`forecast_offset`) by adding them to each dataset within the tree. `forecast_reference_time` is also captured in the path within the tree (`model_run/{forecast_reference_time}`). I’m also making a dataset at the root with that info.

Looking at how Kerchunk wants to do aggregations, it probably makes the most sense to avoid doing any reshaping of datasets with extra data. `xarray-fmrc` could depend on just the `forecast_reference_time` in the path and derive all the related info as needed for the reshaping it’s doing (well, as long as there is a reasonable time dimension, but I think even that is manageable). That way Kerchunk can refer to existing model data as is and we can still access them in the various FMRC-ish and climpred ways.

---

<div class="post-metadata">

**Author:** ![rsignell](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/rsignell/32/447_2.png) [@rsignell](https://discourse.pangeo.io/u/rsignell)\
**Post date:** [May 26, 2023, 5:44am UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/13 "2023-05-26T05:44:44Z")

</div>

Slightly off topic, but just in case it’s useful here in some way, here is [my GEFS kerchunking notebook](https://nbviewer.org/gist/rsignell-usgs/36742ffb32d8fe61938aa1fb116d9147).

---

<div class="post-metadata">

**Author:** ![rbrady](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/rbrady/32/2164_2.png) [@rbrady](https://discourse.pangeo.io/u/rbrady)\
**Post date:** [May 30, 2023, 4:35pm UTC](https://discourse.pangeo.io/t/xarray-and-collections-of-forecasts/3054/14 "2023-05-30T16:35:15Z")

</div>

Love the approach you’re taking with datatree @abkfenris. When we developed `climpred`, datatree wasn’t around. In hindsight, we effectively created something like it with our nested dictionary structure to hold forecasts, observations, perfect model runs, etc. We also hacked together our classes to apply the same function across all objects ([climpred/classes.py at main · pangeo-data/climpred · GitHub](https://github.com/pangeo-data/climpred/blob/main/climpred/classes.py) I think under ` __getattr__ `). It would be great to make it so `climpred` accepts datatree objects. Although I left for an industry job without much time/legality for open-source work, as did @aaronspring I believe. Would be great for others to pick up the torch if the package still provides value for the community.