1

I'm using .nc files with spatially referenced timeseries data in xarray, i.e. inherently 3-dimensional data (lat, lon, time). My issue is, that the dataset's coordinates lat and lon are merged to a site_id, therefore my xarray has the coordinates site_id and time, and lat and lon coordinates are stored as data variables:

<xarray.Dataset>
Dimensions:      (time: 166560, site_id: 2751) 
Coordinates:
    time         (time) datetime64[ns] 2000-01-01 ... 2018-12-31T23:00:00
    site_id      (site_id) float64 20.0 21.0 22.0 ... 2.769e+03 2.77e+03
Data variables:
    data         (site_id, time) float64 0.5698 0.5349 0.4891 ... 0.5762 0.6608
    lat          (site_id) float64 36.77 37.21 38.52 38.96 ... 35.02 35.45 35.88
    lon          (site_id) float64 -9.008 -9.131 -9.512 ... 34.39 34.54 34.69

What I want to achieve is a restructured dataset, where time, lat and lon are set as coordinates, and the original index site_id is dropped:

<xarray.Dataset>
Dimensions:      (time: 166560, lat: 2751, lon: 2751) 
Coordinates:
    time         (time) datetime64[ns] 2000-01-01 ... 2018-12-31T23:00:00
    lat           (lat) float64 36.77 37.21 38.52 38.96 ... 35.02 35.45 35.88
    lon           (lon) float64 -9.008 -9.131 -9.512 ... 34.39 34.54 34.69
Data variables:
    data         (time, lat, lon) float64 0.5698 0.5349 0.4891 ... 0.5762 0.6608

Following a comparable thread (https://stackoverflow.com/questions/43015638/xarray-reshape-data-split-dimension), I tried to specify site_id as MultiIndex and unstacking it with ds.set_index(site_id=("lat","lon")).unstack("site_id") which did not work due to a MemoryError:

numpy.core._exceptions._ArrayMemoryError: Unable to allocate 9.17 TiB for an array with shape (166560, 2751, 2751) and data type float64

Is there a way to solve this problem without having 10TB of RAM?

4
  • Could you clarify why you want to create that restructured dataset? Because it will be a huge file thats mostly empty. Mar 14, 2023 at 11:34
  • Why do you think that the file will be mostly emtpy? I need to select the data by lat/lon [i.e. ds.sel(lat=slice(..), lon=slice(..))] and interpolate all timeseries within the set to match a predefined grid [i.e. ds.data.interp_like(grid, method="nearest")]. Once this is done, I can merge this dataset with other datasets also aligned to the grid.
    – Jonny32418
    Mar 14, 2023 at 13:17
  • For each time-slice: In your original dataset, you have 2751 data-points. In the transformed dataset, you have 2,751*2,751 = 7,568,001 data-points of which 7,565,250 points are empty or "missing data". Mar 14, 2023 at 14:01
  • You are right, good point. As mentioned above, I would like to project the data to a defined grid, e.g. 0.25° resolution over some lat/lon range. For simple point data, I would just define the grid an make a spatial join in GIS. But I have no idea what to do with timeseries data. I was thinking of defining the grid as e.g. GeoTIFF and interpolating all timeseries data to the grid's centroids. To do this, I first wanted to replace the site_ids by the lat lon coordinates (as my grid is also expressed in lat/lon coordinates) and perform ds.data.interp_like(grid, method="nearest").
    – Jonny32418
    Mar 14, 2023 at 15:17

1 Answer 1

1

First import modules and create some test data:

import xarray as xr
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.interpolate import griddata

# Determine dimension sizes
nt = 20
nsite_id = 200

# Create dummy coordinates
times = pd.date_range("2000-01-01", "2018-01-01", periods = nt)
site_ids = np.linspace(0, nsite_id-1, nsite_id)
data = np.random.rand(nsite_id, nt)
lats = np.linspace(50, 60, nsite_id)
lons = np.linspace(20, 10, nsite_id)
np.random.shuffle(lats)
np.random.shuffle(lons)

# Create dummy dataset
ds = xr.Dataset({"data": (["site_id", "time"],data), "lat": (["site_id"], lats), "lon": (["site_id"], lons)}, coords = {"time": times, "site_id": site_ids})

print(ds)

Gives the following dataset.

<xarray.Dataset>
Dimensions:  (site_id: 200, time: 20)
Coordinates:
  * time     (time) datetime64[ns] 2000-01-01 ... 2018-01-01
  * site_id  (site_id) float64 0.0 1.0 2.0 3.0 4.0 ... 196.0 197.0 198.0 199.0
Data variables:
    data     (site_id, time) float64 0.04678 0.2404 0.6671 ... 0.00551 0.3914
    lat      (site_id) float64 59.95 52.36 54.27 50.1 ... 55.08 50.35 58.79
    lon      (site_id) float64 19.4 15.28 13.62 16.98 ... 19.5 14.82 18.89 11.31

Then for your first question, you can select data by lat/lon "slices" like this:

# Select data by lat/lon
ds_selected = ds.where((ds.lat > 51) &
                       (ds.lat < 59) &
                       (ds.lon > 11) &
                       (ds.lon < 19), drop = True)

Just to get an idea of what happens when you do the .unstack thing:

test = ds_selected.set_index(site_id=("lat", "lon")).unstack("site_id")
test["data"].isel(time=0).plot()

You can see that there is a lot of missing data:

The data is sampled irregular and there are a lot of empty pixels.

To interpolate our selected data onto a regular grid, we can use scipy.interpolate.griddata, vectorize it using xr.apply_ufunc and allow for dask chunks like this:

# Define function to interpolate sampled data to grid.
def interp_to_grid(u, xc, yc, new_lats, new_lons):
    new_points = np.stack(np.meshgrid(new_lats, new_lons), axis = 2).reshape((new_lats.size * new_lons.size, 2))
    z = griddata((xc, yc), u, (new_points[:,1], new_points[:,0]), method = 'nearest', fill_value = np.nan)
    out = z.reshape((new_lats.size, new_lons.size), order = "F")
    return out 

# Create chunks (you'd probably want to write `ds_selected` to disk 
# first using `.to_netcdf()` and then load it from there as a `dask.array`).
ds_selected = ds_selected.chunk({"site_id": -1, "time": 10})

values = ds_selected.data
lons = ds_selected.lon
lats = ds_selected.lat

# Create some dummy grid on which to interpolate.
_new_lats = np.linspace(ds_selected.lat.min(), ds_selected.lat.max(), 120)
_new_lons = np.linspace(ds_selected.lon.min(), ds_selected.lon.max(), 90)
new_lons = xr.DataArray(_new_lons, dims = "lon", coords = {"lon": _new_lons})
new_lats = xr.DataArray(_new_lats, dims = "lat", coords = {"lat": _new_lats})

# Vectorize the `interp_to_grid` function.
gridded_ds = xr.apply_ufunc(interp_to_grid,
                     values, lons, lats, new_lats, new_lons,
                     vectorize = True,
                     dask = "parallelized",
                     input_core_dims = [['site_id'],['site_id'],['site_id'],["lat"],["lon"]],
                     output_core_dims = [['lat', 'lon']],
                     )

Now we have the data interpolated on a grid. Lets make a plot.

# Create a plot
fig = plt.figure()
ax = fig.gca()
gridded_ds.isel(time=0).plot(ax = ax)
ax.scatter(lons, lats, c = values.isel(time=0), edgecolors="white")

White circles show the sampled data, its seems to match nicely with the interpolated data.

2
  • Thanks Bert! Works super smooth! I converted the gridded_ds from DataArray back to a dataset using gridded_ds.to_dataset(name="data").
    – Jonny32418
    Mar 15, 2023 at 16:47
  • @Jonny32418 feel free to accept the answer and upvote ;) Mar 16, 2023 at 7:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.