Issue on page /notebooks/distributions.html

[remigandoin]

Hi team. Sorry for asking a rookie question here, but I can't seem to figure out how to subset ds_danra_sl using lat and lon.

Seems like lat and lon are 2d arrays, I can't subset these as slices or directly like I would do with other zarr repos such as https://earthdatahub.destine.eu/collections/era5/datasets/reanalysis-era5-single-levels.

How do you compute ds_pt_thor for instance?

Thank you
Rémi

[khintz]

I do see we haven't defined ds_pt_thor in the notebook. That is undoubtedly a mistake. Thanks for bringing this up. We will have a look at it as soon as possible.

[khintz]

A quick investigation revealed that we had forgotten to include the functions used to compute ds_pt_thor.
To get you going, I paste in the relevant code here. Then we will update the notebook asap. I will keep this issue open until that has happened.

def sel_nearest_to_latlon_pt(ds, pt, wrap_lon=True):
    """
    For a gridded Harmonie dataset (ds), find the nearest grid point to a given
    lat/lon point (pt) and return a selection of the dataset with only
    that grid point.

    NOTE: this function internally uses a KDTree for fast nearest neighbour
    lookup, this tree is cached based on the suite_name of the dataset. The
    distance is computed as Euclidean distance in the lat/lon space and so is
    not the actual distance on the earth's surface.

    Parameters
    ----------
    ds : xarray.Dataset
        Harmonie dataset
    pt : dict
        Dictionary with keys "lat" and "lon" specifying the lat/lon point to
        find the nearest grid point to.
    wrap_lon : bool, optional
        Whether to wrap the longitude values to the domain bounds (default: True)
    """
    identifier = ds.attrs.get("suite_name")
    if identifier is None:
        raise ValueError("dataset must have a suite_name attribute")

    if not isinstance(pt, dict):
        raise TypeError("pt_latlon must be a dict with keys 'lat' and 'lon'")

    lon, lat = pt["lon"], pt["lat"]

    @cache_to_pickle(identifier)
    def build_tree():
        logger.info(f"Building nearest-point to lat/lon lookup tree for {identifier}")
        values = list(zip(ds.lon.values.flatten(), ds.lat.values.flatten()))
        tree_kdtree = scipy.spatial.cKDTree(values)
        latlon_bounds = dict(
            lat=(ds.lat.min().values, ds.lat.max().values),
            lon=(ds.lon.min().values, ds.lon.max().values),
        )
        return dict(tree=tree_kdtree, latlon_bounds=latlon_bounds)

    if identifier in LOOKUP_INFO:
        lookup_info = LOOKUP_INFO[identifier]
    else:
        lookup_info = build_tree()
        LOOKUP_INFO[identifier] = lookup_info

    lookup_tree = lookup_info["tree"]
    latlon_bounds = lookup_info["latlon_bounds"]

    if wrap_lon:
        lon_wrapped, applied_lon_offset = _wrap_lon(lon, latlon_bounds["lon"])
    else:
        # check that the longitude is within the bounds
        if lon < latlon_bounds["lon"][0] or lon > latlon_bounds["lon"][1]:
            raise ValueError(
                f"Longitude value ({lon}) is not within the domain bounds ({latlon_bounds['lon']})"
            )

    # check that the latitude is within the bounds
    if lat < latlon_bounds["lat"][0] or lat > latlon_bounds["lat"][1]:
        raise ValueError(
            f"Latitude value ({lat}) is not within the domain bounds ({latlon_bounds['lat']})"
        )

    _, idx_nearest = lookup_tree.query((lon_wrapped, lat))
    nx, ny = ds.lon.shape
    i, j = idx_nearest % ny, idx_nearest // ny

    ds_point = ds.isel(x=i, y=j)

    if wrap_lon:
        ds_point["lon"] = ds_point["lon"] - applied_lon_offset

    return ds_point

pt_thor_windpark = dict(lat=56.330, lon=8.0)
ds_pt_thor = sel_nearest_to_latlon_pt(ds=ds_danra_sl, pt=pt_thor_windpark)

[remigandoin]

Thanks Kasper. I'm now getting issues with "@cache_to_pickle(identifier)" - maybe are there more dependencies, or this directory needs to be set in the main script. I asked chatGPT but it's struggling giving me specific instructions.

I'm not working on something urgent with these data, this can wait until the person in charge is back to the office.

[observingClouds]

Hi @remigandoin,
Please have a look at the code/notebook in the PR above. We need to fix some CI issues to run and cache this rather heavy computation in the CI, but the code should now be complete.

[remigandoin]

Hi @observingClouds, I'm running the following code for downloading data but it takes very long.

import xarray as xr
import numpy as np
import pandas as pd
from dask.diagnostics import ProgressBar
import os

# -----------------------------
# OPEN DATASET (LAZY)
# -----------------------------
ds = xr.open_zarr(
    "s3://dmi-danra-05/height_levels.zarr",
    consolidated=True,
    storage_options={"anon": True},
)

ds.attrs["suite_name"] = "danra"

# -----------------------------
# USER INPUT
# -----------------------------
target_lat = 56.15
target_lon = 10.20

out_dir = "monthly_wind_files"
os.makedirs(out_dir, exist_ok=True)

# -----------------------------
# FIND NEAREST GRID POINT
# -----------------------------
dist = (ds.lat - target_lat) ** 2 + (ds.lon - target_lon) ** 2
dist_stacked = dist.stack(z=("y", "x"))

flat_index = dist_stacked.argmin(dim="z").compute().item()
iy, ix = np.unravel_index(flat_index, (ds.dims["y"], ds.dims["x"]))

print(f"Nearest grid point: y={iy}, x={ix}")
print(
    "Grid lat/lon:",
    float(ds.lat.isel(y=iy, x=ix)),
    float(ds.lon.isel(y=iy, x=ix)),
)

point_ds = ds.isel(y=iy, x=ix)

# -----------------------------
# BUILD MONTHLY TIME SLICES
# -----------------------------
time_index = pd.to_datetime(point_ds.time.values)

months = pd.period_range(
    start=time_index.min(),
    end=time_index.max(),
    freq="M",
)

# -----------------------------
# NETCDF ENCODING
# -----------------------------
encoding = {
    "wind_speed": {"zlib": True, "complevel": 4},
    "u": {"zlib": True, "complevel": 4},
    "v": {"zlib": True, "complevel": 4},
}

# -----------------------------
# LOOP OVER MONTHS
# -----------------------------
for p in months:
    fname = f"wind_timeseries_{p.strftime('%Y-%m')}.nc"
    out_path = os.path.join(out_dir, fname)

    if os.path.exists(out_path):
        print(f"Skipping {p} (file exists)")
        continue

    print(f"Processing {p}...")

    month_ds = point_ds.sel(
        time=slice(p.start_time, p.end_time)
    )

    if month_ds.time.size == 0:
        continue

    wind_speed = np.sqrt(
        month_ds["u"] ** 2 + month_ds["v"] ** 2
    ).rename("wind_speed")

    out_ds = xr.Dataset(
        {
            "wind_speed": wind_speed,
            "u": month_ds["u"],
            "v": month_ds["v"],
        }
    )

    with ProgressBar():
        out_ds.to_netcdf(
            out_path,
            mode="w",
            format="NETCDF4",
            encoding=encoding,
            compute=True,
        )

print(f"\n✅ Monthly files saved in: {out_dir}")

[observingClouds]

Hi @remigandoin 👋
The script looks OK to me and running it I achieve similar speeds. When downloading single station data, you need to keep in mind that the dataset is chunked primarily for users retrieving spatial fields. Station data is included in these chunks and is extracted on the client side. The DANRA dataset bucket mentioned above is chunked along the time dimension, with each chunk containing about a month (256 timesteps @ 3 h) of output. As you request 2 variables and all levels, you request 9 altitude chunks @ 40MB x 2 variables = 720MB. Depending on your internet connection, e.g. 100mbps, this can result at roughly a minute per month.

[leifdenby]

The script looks OK to me and running it I achieve similar speeds. When downloading single station data, you need to keep in mind that the dataset is chunked primarily for users retrieving spatial fields. Station data is included in these chunks and is extracted on the client side. The DANRA dataset bucket mentioned above is chunked along the time dimension, with each chunk containing about a month (256 timesteps @ 3 h) of output. As you request 2 variables and all levels, you request 9 altitude chunks @ 40MB x 2 variables = 720MB. Depending on your internet connection, e.g. 100mbps, this can result at roughly a minute per month.

@SimonKamuk and I have been discussing whether we could ask to be allowed to host two mirrors of DANRA, one a) chunked for x/y (i.e. station-like) access with long chunks in time small in space, and b) one chunked with small chunks in time and large in x/y. What we have now tries to be a compromise with 3 by 2 chunks in space and then having the chunks size in time be so that every chunk is ~30MB. Happy to hear of you have any thoughts on this @observingClouds and @remigandoin :)

[remigandoin]

For most users in renewables and the energy sector, the primary use case of reanalysis dataset is using long-term time series at some selected locations - anything that that can help make this easier (such as https://earthdatahub.destine.eu/collections/era5/datasets/reanalysis-era5-single-levels) is welcome :)

Happy to discuss this further as part of the work we do within https://www.wemcouncil.org/wp/reanalysis-working-group/

[leifdenby]

For most users in renewables and the energy sector, the primary use case of reanalysis dataset is using long-term time series at some selected locations - anything that that can help make this easier (such as https://earthdatahub.destine.eu/collections/era5/datasets/reanalysis-era5-single-levels) is welcome :)

Yes, that was exactly what I was thinking too. This is why I suggested that we could maybe create mirror of this dataset which has longer chunks in time, and smaller in space. But, we do also have downstream users which need the opposite though (creating 2D spatial plots at a given altitude, for example windspeed over an area at hub hight). @observingClouds is also suggesting we could look into the sharding support in zarr3.

Interestingly the ERA5 copy you linked to actually has quite big chunks in space (64x64), but obviously not as large as ours (295x263), so you are probably ending up downloading ~20 times more data than you currently need (if we had similar chunking to that ERA5 dataset) On the other hand that would make it a lot worse for people wanting to do 2D spatial plots :)

<xarray.Dataset> Size: 404TB
Dimensions:     (valid_time: 753888, latitude: 721, longitude: 1440)
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 6MB 1940-01-01 ... 2025-12-31T23:...
  * latitude    (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 12kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
    number      int64 8B ...
    surface     float64 8B ...
Data variables: (12/129)
    alnid       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    alnip       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    aluvd       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    aluvp       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    anor        (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    asn         (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    ...          ...
    viiwn       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    vilwd       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    vilwe       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    vilwn       (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    z           (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
    zust        (valid_time, latitude, longitude) float32 3TB dask.array<chunksize=(4320, 64, 64), meta=np.ndarray>
Attributes:
    Conventions:             CF-1.7
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_edition:            1
    GRIB_subCentre:          0
    history:                 2025-12-04T16:25 GRIB to CDM+CF via cfgrib-0.9.1...
    institution:             European Centre for Medium-Range Weather Forecasts

Again, thanks for your feedback and we'll look into what we can do :)

[remigandoin]

Thanks. Some secondary use cases require mean wind speed maps, but this is typically for high level studies. Most renewable energy users will like to use time series, happy to discuss about their requirements and also help you at team DMI get in touch with these users. It's a very diverse community but you can reach many of them via https://groups.io/g/wrag/ for instance. In Denmark, DTU Wind and the team behind https://map.neweuropeanwindatlas.eu/ are also great contacts to have.

[khintz]

Hi @remigandoin, we (DMI) have regular contact with DTU Wind and also host visiting scientists from DTU Wind to work with us on a parameterization scheme for wind farms in our weather model. But that said, thanks a lot for the suggestions! :) I think it will be very valuable to have more onboard in a discussion on new chunking.