I have a 6d DataArray named data and want to use six 2d variables to interpolate the data.
The example netcdf file is here.
Data variables:
interp_var1 (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan
interp_var2 (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan
interp_var3 (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan
interp_var4 (y, x) float32 nan nan nan nan nan nan ... nan nan nan nan nan
interp_lon (y, x) float32 122.9 122.9 123.0 123.1 ... 95.07 94.95 94.82
interp_lat (y, x) float32 82.24 82.27 82.3 82.33 ... 86.48 86.53 86.58
data (var1, var2, var3, var4, lon, lat) float64 ...
Interp()
First, I tried the interp() function:
import xarray as xr
import matplotlib.pyplot as plt
ds_test = xr.open_dataset('test_multi_interp.nc')
interp_dict = {'var1': ds_test['interp_var1'],
'var2': ds_test['interp_var2'],
'var3': ds_test['interp_var3'],
'var4': ds_test['interp_var4'],
'lon': ds_test['interp_lon'],
'lat': ds_test['interp_lat'],
}
ds_test['data'].interp(interp_dict).plot(vmin=0)
There're many missing data as shown on the left panel. So, I tried to use extrapolation.
Interp with extrapolation
ds_test['data'].interp(interp_dict, kwargs={"fill_value": None}).plot(vmin=0)
Many missing data are filled, but the large white strips are still there.
Cause
So, I tried to find out the cause of the white stripes by neglecting one dimension of interpolation at each time.
When I drop any var (var1~var4) or lon from the interp_dict, the white stripes always exist:
novar4_dict = {i:interp_dict[i] for i in interp_dict if i!='var4'}
ds_test['data'].interp(novar4_dict, kwargs={"fill_value": None}).isel(var4=slice(0,6)).plot(col='var4', col_wrap=3, vmin=0)
Finally, the root is the lat dimension:
nolat_dict = {i:interp_dict[i] for i in interp_dict if i!='lat'}
ds_test['data'].interp(nolat_dict, kwargs={"fill_value": None}).isel(lat=slice(16,22)).plot(col='lat', col_wrap=3, vmin=0)
Let's take the first white stripe as an example.
The lat is from 83 to 84 in that region (as shown below) and the extrapolation should work well because there're valid data values between 81.5 and 83.5 (as shown above).
How to use the interpolation to fill the stripes correctly?
Attempt 1
I tried to apply interp() on lat dimension after the multi-dimension interpolation:
ds_test['data'].interp(nolat_dict, kwargs={"fill_value": None}).interp(lat=ds_test['interp_lat'])
But, it failed:
File ~/miniconda3/lib/python3.9/site-packages/xarray/core/utils.py:459, in Frozen.__getitem__(self, key)
458 def __getitem__(self, key: K) -> V:
--> 459 return self.mapping[key]
KeyError: 'y'
Attempt 2
Iterating the interpolation result without lat and applying interp1d manually gives similar results.
import numpy as np
from scipy.interpolate import interp1d
da_interp = ds_test['data'].interp(nolat_dict, kwargs={"fill_value": None})
result = np.empty(da_interp.shape[:-1])
for i in range(da_interp.sizes['y']):
for j in range(da_interp.sizes['x']):
f = interp1d(da_interp.lat, da_interp.isel(y=i, x=j), axis=-1, kind='linear', fill_value='extrapolate')
result[i, j] = f(ds_test['interp_lat'].sel(y=i, x=j))
xr.DataArray(result).rename({'dim_0': 'y', 'dim_1': 'x'}).plot(vmin=0)
Let's pick one grid to investigate why it fails:
da_interp.isel(y=30, x=40).plot(marker='o')
plt.axvline(ds_test['interp_lat'].isel(y=30, x=40), c='red', label='interp_lat at y=30,x=40')
plt.legend()
plt.title('da_interp at y=30,x=40')
i = 30; j = 40
f = interp1d(da_interp.lat, da_interp.isel(y=i, x=j), axis=-1, kind='linear', fill_value='extrapolate')
print(da_interp.isel(y=i, x=j))
print(f(ds_test['interp_lat'].sel(y=i, x=j)))
Although there're three valid values before the interp_lat, interp_1d returns nan.
<xarray.DataArray 'data' (lat: 25)>
array([ nan, nan, nan, nan,
nan, nan, nan, nan,
nan, nan, nan, nan,
nan, -2.35975262e-06, nan, nan,
6.67684124e-06, 6.35597345e-06, 3.44029654e-06, nan,
nan, nan, nan, nan,
nan])
Even if I interpolate data along lat dim and then apply interp1d:
i = 30; j = 40
f = interp1d(da_interp.lat, da_interp.isel(y=i, x=j).interpolate_na(dim='lat'), axis=-1, kind='linear', fill_value='extrapolate')
print(f(ds_test['interp_lat'].sel(y=i, x=j)))
print(da_interp.isel(y=i, x=j).interpolate_na(dim='lat'))
It's still nan:
nan
<xarray.DataArray 'data' (lat: 25)>
array([ nan, nan, nan, nan,
nan, nan, nan, nan,
nan, nan, nan, nan,
nan, -2.35975262e-06, 6.52445336e-07, 3.66464329e-06,
6.67684124e-06, 6.35597345e-06, 3.44029654e-06, nan,
nan, nan, nan, nan,
nan])