将xarray对象重新采样以降低空间分辨率。

15

使用 xarray 将数据重采样为较低的空间分辨率

我想将我的 xarray 对象重采样为较低的空间分辨率(即更少的像素)。

import pandas as pd
import numpy as np
import xarray as xr

time = pd.date_range(np.datetime64('1998-01-02T00:00:00.000000000'), np.datetime64('2005-12-28T00:00:00.000000000'), freq='8D')
x = np.arange(1200)
y = np.arange(1200)

latitude = np.linspace(40,50,1200)
longitude = np.linspace(0,15.5572382,1200)

latitude, longitude = np.meshgrid(latitude, longitude)

BHR_SW = np.ones((365, 1200, 1200))

output_da = xr.DataArray(BHR_SW, coords=[time, y, x])
latitude_da = xr.DataArray(latitude, coords=[y, x])
longitude_da = xr.DataArray(longitude, coords=[y, x])

output_da = output_da.rename({'dim_0':'time','dim_1':'y','dim_2':'x'})
latitude_da = latitude_da.rename({'dim_0':'y','dim_1':'x'})
longitude_da = longitude_da.rename({'dim_0':'y','dim_1':'x'})

output_ds = output_da.to_dataset(name='BHR_SW')
output_ds = output_ds.assign({'latitude':latitude_da, 'longitude':longitude_da})

print(output_ds)

<xarray.Dataset>
Dimensions:    (time: 365, x: 1200, y: 1200)
Coordinates:
  * time       (time) datetime64[ns] 1998-01-02 1998-01-10 ... 2005-12-23
  * y          (y) int64 0 1 2 3 4 5 6 7 ... 1193 1194 1195 1196 1197 1198 1199
  * x          (x) int64 0 1 2 3 4 5 6 7 ... 1193 1194 1195 1196 1197 1198 1199
Data variables:
    BHR_SW     (time, y, x) float64 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    latitude   (y, x) float64 40.0 40.01 40.02 40.03 ... 49.97 49.98 49.99 50.0
    longitude  (y, x) float64 0.0 0.0 0.0 0.0 0.0 ... 15.56 15.56 15.56 15.56
```

我的问题是,如何按照x,y坐标重新采样成200x200网格?

这是对变量空间分辨率的降低。

我尝试过以下方法:

output_ds.resample(x=200).mean()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-54-10fbdf855a5d> in <module>()
----> 1 output_ds.resample(x=200).mean()

/home/mpim/m300690/miniconda3/envs/holaps/lib/python2.7/site-packages/xarray/core/common.pyc in resample(self, indexer, skipna, closed, label, base, keep_attrs, **indexer_kwargs)
    701         group = DataArray(dim_coord, coords=dim_coord.coords,
    702                           dims=dim_coord.dims, name=RESAMPLE_DIM)
--> 703         grouper = pd.Grouper(freq=freq, closed=closed, label=label, base=base)
    704         resampler = self._resample_cls(self, group=group, dim=dim_name,
    705                                        grouper=grouper,

/home/mpim/m300690/miniconda3/envs/holaps/lib/python2.7/site-packages/pandas/core/resample.pyc in __init__(self, freq, closed, label, how, axis, fill_method, limit, loffset, kind, convention, base, **kwargs)
   1198                              .format(convention))
   1199
-> 1200         freq = to_offset(freq)
   1201
   1202         end_types = set(['M', 'A', 'Q', 'BM', 'BA', 'BQ', 'W'])

/home/mpim/m300690/miniconda3/envs/holaps/lib/python2.7/site-packages/pandas/tseries/frequencies.pyc in to_offset(freq)
    174                     delta = delta + offset
    175         except Exception:
--> 176             raise ValueError(libfreqs._INVALID_FREQ_ERROR.format(freq))
    177
    178     if delta is None:

ValueError: Invalid frequency: 200

但是我遇到了所示的错误。

我该如何完成x和y的空间重采样?

理想情况下,我想要这样做:

output_ds.resample(x=200, y=200).mean()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-55-e0bfce19e037> in <module>()
----> 1 output_ds.resample(x=200, y=200).mean()

/home/mpim/m300690/miniconda3/envs/holaps/lib/python2.7/site-packages/xarray/core/common.pyc in resample(self, indexer, skipna, closed, label, base, keep_attrs, **indexer_kwargs)
    679         if len(indexer) != 1:
    680             raise ValueError(
--> 681                 "Resampling only supported along single dimensions."
    682             )
    683         dim, freq = indexer.popitem()

ValueError: Resampling only supported along single dimensions.

注意:真实数据具有不同的行为

这是我上面创建的测试数据。在从netcdf文件中读取的真实数据中会有所不同。

<xarray.Dataset>
Dimensions:    (time: 368, x: 1200, y: 1200)
Coordinates:
  * time       (time) datetime64[ns] 1998-01-02 1998-01-10 ... 2005-12-28
Dimensions without coordinates: x, y
Data variables:
    latitude   (y, x) float32 ...
    longitude  (y, x) float32 ...
    Data_Mask  (y, x) float32 ...
    BHR_SW     (time, y, x) float32 ...
Attributes:
    CDI:               Climate Data Interface version 1.9.5 (http://mpimet.mp...
    Conventions:       CF-1.4
    history:           Fri Dec 07 13:29:13 2018: cdo mergetime GLOBALBEDO/Glo...
    content:           extracted variabel BHR_SW of the original GlobAlbedo (...
    metadata_profile:  beam
    metadata_version:  0.5
    CDO:               Climate Data Operators version 1.9.5 (http://mpimet.mp...
```

我尝试过类似的事情:

ds.resample(x=200).mean()

/home/mpim/m300690/miniconda3/envs/holaps/lib/python2.7/site-packages/xarray/core/common.pyc in resample(self, indexer, skipna, closed, label, base, keep_attrs, **indexer_kwargs)
    686         dim_coord = self[dim]
    687
--> 688         if isinstance(self.indexes[dim_name], CFTimeIndex):
    689             raise NotImplementedError(
    690                 'Resample is currently not supported along a dimension '

/home/mpim/m300690/miniconda3/envs/holaps/lib/python2.7/site-packages/xarray/core/coordinates.pyc in __getitem__(self, key)
    309         if key not in self._sizes:
    310             raise KeyError(key)
--> 311         return self._variables[key].to_index()
    312
    313     def __unicode__(self):

KeyError: 'x'

非常感谢任何帮助。


1
请参考以下链接:https://dev59.com/NFgQ5IYBdhLWcg3wolao#42463491 - b-fg
4个回答

17
最近,xarray添加了coarsen方法。我认为这是空间降采样的最佳方法,即使不能使用它设置所需的最终分辨率并自动计算。coarsen将在非重叠窗口上执行操作(平均、最大、最小等),根据您设置的窗口大小,您将获得所需的最终分辨率。原始输入数据来自作者:
import pandas as pd
import numpy as np
import xarray as xr

​

time = pd.date_range(np.datetime64('1998-01-02T00:00:00.000000000'), np.datetime64('2005-12-28T00:00:00.000000000'), freq='8D')
x = np.arange(1200)
y = np.arange(1200)


latitude = np.linspace(40,50,1200)
longitude = np.linspace(0,15.5572382,1200)
latitude, longitude = np.meshgrid(latitude, longitude)

BHR_SW = np.ones((365, 1200, 1200))

output_da = xr.DataArray(BHR_SW, coords=[time, y, x])
latitude_da = xr.DataArray(latitude, coords=[y, x])
longitude_da = xr.DataArray(longitude, coords=[y, x])
output_da = output_da.rename({'dim_0':'time','dim_1':'y','dim_2':'x'})
latitude_da = latitude_da.rename({'dim_0':'y','dim_1':'x'})
longitude_da = longitude_da.rename({'dim_0':'y','dim_1':'x'})

output_ds = output_da.to_dataset(name='BHR_SW')
output_ds = output_ds.assign({'latitude':latitude_da, 'longitude':longitude_da})
print(output_ds)

​

<xarray.Dataset>
Dimensions:    (time: 365, x: 1200, y: 1200)
Coordinates:
  * time       (time) datetime64[ns] 1998-01-02 1998-01-10 ... 2005-12-23
  * y          (y) int64 0 1 2 3 4 5 6 7 ... 1193 1194 1195 1196 1197 1198 1199
  * x          (x) int64 0 1 2 3 4 5 6 7 ... 1193 1194 1195 1196 1197 1198 1199
Data variables:
    BHR_SW     (time, y, x) float64 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    latitude   (y, x) float64 40.0 40.01 40.02 40.03 ... 49.97 49.98 49.99 50.0
    longitude  (y, x) float64 0.0 0.0 0.0 0.0 0.0 ... 15.56 15.56 15.56 15.56

为了将空间分辨率从1200x1200降低到200x200,我们需要使用6x6的窗口来粗化数据。

output_ds.coarsen(x=6).mean().coarsen(y=6).mean()
# or output_ds.coarsen(x=6,y=6).mean()
<xarray.Dataset>
Dimensions:    (time: 365, x: 200, y: 200)
Coordinates:
  * time       (time) datetime64[ns] 1998-01-02 1998-01-10 ... 2005-12-23
  * y          (y) float64 2.5 8.5 14.5 20.5 ... 1.184e+03 1.19e+03 1.196e+03
  * x          (x) float64 2.5 8.5 14.5 20.5 ... 1.184e+03 1.19e+03 1.196e+03
Data variables:
    BHR_SW     (time, y, x) float64 1.0 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0 1.0
    latitude   (y, x) float64 40.02 40.07 40.12 40.17 ... 49.88 49.93 49.98
    longitude  (y, x) float64 0.03244 0.03244 0.03244 ... 15.52 15.52 15.52

14

更新

@clausmichele的回答使用coarsen现在是实现此目的的最佳方法。请注意,coarsen现在包括指定所需输出坐标的功能。

原始帖子

正如piman314所建议的,groupby是在xarray中执行此操作的唯一方法。Resample仅适用于日期时间坐标。

由于xarray当前不处理多维groupby,因此必须分两个阶段完成:

# this results in bin centers on 100, 300, ...
reduced = (
    output_ds
    .groupby(((output_ds.x//200) + 0.5) * 200)
    .mean(dim='x')
    .groupby(((output_ds.y//200) + 0.5) * 200)
    .mean(dim='y'))

如果您只想降采样数据,可以使用位置切片:

output_ds[:, ::200, ::200]

或者,使用命名维度:

output_ds[{'x': slice(None, None, 200), 'y': slice(None, None, 200)}]

最后,还有其他专门为与xarray兼容的快速重采样而设计的软件包。 xESMF 是其中一个不错的选择。


5
现在使用groupby_bins或甚至coarsen会使这变得稍微容易些。 - astrojuanlu
如何指定所需的输出坐标?我的理解是,您仍然需要为每个要粗化的维度传递一个窗口。 - HoosierDaddy
同样想知道,我该如何指定所需的输出坐标? - Rob

3

由于你正在使用已经通过CDO处理过的NetCDF文件,因此你也可以使用CDO的SAMPLEGRID函数或NCO的bilinear_interp函数:

SAMPLEGRID(https://code.mpimet.mpg.de/projects/cdo/embedded/cdo.pdf) 不进行插值,只是删除每个n-th网格点。

bilinear_interp(http://nco.sourceforge.net/nco.html#Bilinear-interpolation) 进行插值。

由于你可能需要平均值、最大值或其他反照率值,因此你可能更喜欢NCO的bilinear_interp。但是CDO的SAMPLEGRID可以为NOC的bilinear_interp提供所需的grid_out


2

使用xarray进行操作的最明显的方法是使用groupby_bins,但事实证明这样做非常缓慢。更有效的方法可能是使用numpy和超快速索引([:, :, frequency])。

nsamples = 200
bins = np.linspace(output_ds.x.min(),
                   output_ds.x.max(), nsamples).astype(int)
output_ds = output_ds.groupby_bins('x', bins).first()

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接