2017-08-04 45 views
0

雖然看起來應該是相當簡單的,但我對這個主題有困難。Xarray:沒有尺寸的切片座標

我想使用一組緯度和經度座標對xarray數據集進行切片。

這裏是我的數據集的樣子:

In [31]: data = xr.open_mfdataset(open_file, decode_cf=True) 

In [32]: data 
Out[32]: 
<xarray.Dataset> 
Dimensions: (time: 108120, x: 349, y: 277) 
Coordinates: 
    lons  (y, x) float64 -145.5 -145.3 -145.1 -144.9 -144.8 -144.6 -144.4 ... 
    lats  (y, x) float64 1.0 1.104 1.208 1.312 1.416 1.519 1.621 1.724 ... 
    * time  (time) datetime64[ns] 1980-01-01 1980-01-01T03:00:00 ... 
Dimensions without coordinates: x, y 
Data variables: 
    stp  (time, y, x) float64 0.1235 0.0867 0.07183 0.05389 0.05901 ... 

這是我做切片:

In [48]: lat_bnd = [25,30] 
    ...: lon_bnd = [-80,-75] 

In [49]: r = data.sel(y=slice(*lat_bnd),x=slice(*lon_bnd)) 

,一切似乎很大:

In [50]: r 
Out[50]: 
    <xarray.Dataset> 
    Dimensions: (time: 108120, x: 5, y: 5) 
    Coordinates: 
     lons  (y, x) float64 -82.52 -82.28 -82.05 -81.81 -81.57 -82.44 -82.2 ... 
     lats  (y, x) float64 13.54 13.46 13.38 13.3 13.22 13.77 13.69 13.61 ... 
     * time  (time) datetime64[ns] 1980-01-01 1980-01-01T03:00:00 ... 
    Dimensions without coordinates: x, y 
    Data variables: 
     stp  (time, y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 

但我的緯度/ lon值不正確:

In [53]: r.lats.values 
Out[53]: 
array([[ 13.53542397, 13.45647916, 13.37686013, 13.296571 , 
     13.21561592], 
     [ 13.76719053, 13.6878189 , 13.60776989, 13.52704767, 
     13.44565641], 
     [ 13.99938176, 13.91958109, 13.83909988, 13.75794233, 
     13.67611265], 
     [ 14.2319952 , 14.15176326, 14.07084762, 13.98925249, 
     13.90698214], 
     [ 14.46502833, 14.3843629 , 14.30301059, 14.22097564, 
     14.13826236]]) 

In [54]: r.lons.values 
Out[54]: 
array([[-82.52229969, -82.28438922, -82.0469968 , -81.8101255 , 
     -81.57377834], 
     [-82.44118948, -82.20260881, -81.96455096, -81.72701901, -81.490016 ], 
     [-82.3595596 , -82.12030558, -81.8815792 , -81.64338357, 
     -81.40572174], 
     [-82.27740522, -82.03747469, -81.79807668, -81.55921433, 
     -81.32089068], 
     [-82.19472148, -81.95411126, -81.71403851, -81.47450637, -81.2355179 ]]) 

當然,如果我嘗試使用lats/lons座標進行切片,由於尺寸不匹配,所以出現錯誤。

In [55]: r = data.sel(lats=slice(*lat_bnd),lons=slice(*lon_bnd)) 
--------------------------------------------------------------------------- 
ValueError        Traceback (most recent call last) 
<ipython-input-55-7c6237be5f22> in <module>() 
----> 1 r = data.sel(lats=slice(*lat_bnd),lons=slice(*lon_bnd)) 

/lib/anaconda2/lib/python2.7/site-packages/xarray/core/dataset.pyc in sel(self, method, tolerance, drop, **indexers) 
    1204   """ 
    1205   pos_indexers, new_indexes = indexing.remap_label_indexers(
-> 1206    self, indexers, method=method, tolerance=tolerance 
    1207  ) 
    1208   result = self.isel(drop=drop, **pos_indexers) 

/lib/anaconda2/lib/python2.7/site-packages/xarray/core/indexing.pyc in remap_label_indexers(data_obj, indexers, method, tolerance) 
    275  new_indexes = {} 
    276 
--> 277  dim_indexers = get_dim_indexers(data_obj, indexers) 
    278  for dim, label in iteritems(dim_indexers): 
    279   try: 

/lib/anaconda2/lib/python2.7/site-packages/xarray/core/indexing.pyc in get_dim_indexers(data_obj, indexers) 
    243  if invalid: 
    244   raise ValueError("dimensions or multi-index levels %r do not exist" 
--> 245       % invalid) 
    246 
    247  level_indexers = defaultdict(dict) 

ValueError: dimensions or multi-index levels ['lons', 'lats'] do not exist 

有沒有什麼我在我的理解中缺少這是一個NARR數據集?

回答

0

在您的第一個示例中,您不是按緯度/經度,而是按xy的數字索引進行索引。也就是說,您正在分割第25和第30個y和第80和第75個x的值。這解釋了爲什麼緯度/經度值在輸出中沒有意義。

您可以通過使用xr.Dataset.set_index()像這樣與你的尺寸的關聯座標:

data.set_index(y='lats', inplace=True) 
data.set_index(x='lons', inplace=True) 
+0

不幸的是,我得到以下錯誤:NotImplementedError:> 1個NDIM直言不支持在這個時候。我有xarray版本0.9.6。我的問題是NARR的lats和lons有兩個維度(x和y)。有沒有其他的見解? –

+0

@MariaMolina我對NARR數據並不熟悉,也沒有類似的數據集來測試,但是您可以嘗試將列表傳遞給'set_index':沿着'data.set_index(y = ['lats', 'lons'],inplace = True)' – Dan