2017-09-01 5 views
0

私は大きな(5GB)温度netCDFファイルを持っています。ファイルには、時間、圧力レベル、緯度、経度の4つの次元があります。NetCDFファイルをサブセット化してタプルを返す

データセットには31のタイムポイントがあり、5つの圧力レベルにしか興味がありません。

私のパラメータは温度tです:

from netCDF4._netCDF4 import Dataset 
# Load the dataset 
dataset = Dataset(path) 
factor = dataset.variables['t'] 

中心セルの周りに私のfactor変数からの温度データの「キューブ」を抽出するには、私は単純にサブセット化を行うだろう、このような:

radius = 5 
# +1 because the subsetting does not include last index 
lats_bounds = [nearest_latitude_index-radius,nearest_latitude_index+radius + 1] 
lons_bounds = [nearest_longitude_index-radius,nearest_longitude_index+radius +1] 

#all timepoints 
times_bounds = [0, len(times)] 

#just the last 5 pressure levels 
pressure_level_bounds = [len(levels)-5, len(levels)] 

results = factor[times_bounds[0]:times_bounds[1],pressure_level_bounds[0]:pressure_level_bounds[1], lats_bounds[0]:lats_bounds[1],lons_bounds[0]:lons_bounds[1]] 

resultsは、タイプndarrayの形状が(31,5,11,11)で、サイズが18755(31 * 5 * 11 * 11)で、すべてのインデックスが単一の値を保持するという問題があります。

私はresultsの値が必要ですが、値ごとに対応するタイムポイント、圧力レベル、緯度と経度も必要です。

理想的には、...このような何か私が行ったようにサブセット化を行うにはしたいと思いますが、私の最終的な結果は、タプルの配列のようになります。どのように私はこれを達成することができます

# corresponding timestamp, pressure level, latitude, longitude 
# and the temperature value extracted. 
final = [ 
(2342342, 1000, 24.532, 53.531, 277), 
(2342342, 1000, 74.453, 26.123, 351), 
(2342342, 1000, 80.311, 56,345, 131), 
... 
] 

を?

答えて

-2

このタスクにはPandasを使用します。しかし、あなたは35倍と5倍の圧力しかないので、私は最初にあなたのアプローチを単純化し、単一の時間と圧力レベルと単一の緯度経度をどうやって行うのかを考え出します。次に、それらのインデックスをループしてタプルを得る方法を見つけます。次のようなものがあります:

for i in range(0, len(times)): 
    for j in range(0, len(levels): 
    print(results[i, j, nearest_lat_idx, nearest_lon_idx)) 

もちろん、緯度と経度のループも追加できますが、醜いです。

1

xarrayのiselを参照してください。

ds = xr.open_dataset(path) 
factor = ds['t'] 

# note that levels/lon/lat are the names of dimensions in your Dataset 
subset = factor.isel(levels=slice(-5, None), 
        lon=[1, 18, 48, 99], lat=[16, 28, 33, 35]) 
stacked = subset.stack(points=('time', 'levels', 'lon', 'lat')) 

# This subset can be converted to a `pandas.Series`: 
data = stacked.to_pandas() 

# or it can be converted to a list of tuples 
df = data.reset_index() 
final = [tuple(row[1].values) for row in df.iterrows()] 

Xarrayもラベルベースのインデクサーをサポートしています(すなわち、lat=[29.3, 42.3])が、そのために、あなたはiselの代わりにselメソッドを使用する必要があります。netCDF4から構文を翻訳することは、このようになります。

+0

ありがとうございました!私は、あなたがあなたのファクタに「時間」を含めていないことに気づいた。isel() 'それは意図的でしたか?すべてのタイムポイントを使用する必要がありますが、将来はサブセットが必要になることがあります。タイムポイントはUnixのタイムスタンプ形式です。 – pookie

+0

時間軸に沿ってスライスする場合は 'time = slice(start、end)'を追加することができますが、現在のユースケースではすべてのタイムステップが与えられます。 – jhamman

関連する問題