2017-05-09 61 views
2

我有兩種格式(時間,lats,lons)的液態水當量厚度月度全球網格數據集。兩者具有相同的空間和時間分辨率。我想關聯它們,但numpy.corrcoef()只適用於二維數組,而不適用於三維。 所以我想關聯整個時間序列的兩個變量的同一個網格點(x,y)。事實上,我想要一個具有相關係數網格的新文件。關聯python中的網格數據集

import numpy as np 
from netCDF4 import Dataset 

wdir = '.../Data/' 

# read GRACE NCs 
GRACE_GFZ = Dataset(wdir+'GRACE/GRCTellus.GFZ.200204_201607.nc','r') 
GRACE_JPL = Dataset(wdir+'GRACE/GRCTellus.JPL.200204_201607.nc','r') 

兩個變量(gfz和jpl)在相同位置具有相同數量的缺失值。

GRACE_GFZ.variables['lwe_thickness'] 
    <type 'netCDF4._netCDF4.Variable'> 
    float32 lwe_thickness(time, lat, lon) 
     long_name: Liquid_Water_Equivalent_Thickness 
     units: cm 
     _FillValue: 32767.0 
     missing_value: 32767.0 
    unlimited dimensions: time 
    current shape = (155, 72, 144) 
    filling off 

GRACE_JPL.variables['lwe_thickness'] 
    <type 'netCDF4._netCDF4.Variable'> 
    float32 lwe_thickness(time, lat, lon) 
     long_name: Liquid_Water_Equivalent_Thickness 
     units: cm 
     _FillValue: 32767.0 
     missing_value: 32767.0 
    unlimited dimensions: time 
    current shape = (155, 72, 144) 
    filling off 

由於它們具有相同的時間和空間分辨率,因此它們的時間,經度和緯度都可以用於兩者。

time = GRACE_GFZ.variables['time'][:] 
lons = GRACE_GFZ.variables['lon'][:] 
lats = GRACE_GFZ.variables['lat'][:] 
gfz = GRACE_GFZ.variables['lwe_thickness'][:] 
jpl = GRACE_JPL.variables['lwe_thickness'][:] 

這裏我想穿過網格並將corrcoef放入數組中。這給了我一堆2x2矩陣。

test = [] 
for x in range(len(lats)): 
    for y in range(len(lons)): 
     print(np.corrcoef(gfz[:,x,y],jpl[:,x,y])) 

我怎樣才能把每個點的corrcoef放到一個新的陣列在正確的位置?

+1

你能用[最小,完整和可驗證的例子](http://stackoverflow.com/help/mcve)更新你的問題嗎? – Kasramvd

+0

這是你需要的信息@Kasramvd? –

回答

0

我克服我的問題有以下:

temp =[] 
corrcoefMatrix_gfzjpl = [[0 for i in range(len(lons))] for j in range(len(lats))] 
for x in range(len(lats)): 
    for y in range(len(lons)): 
     temp = np.corrcoef(gfz[:,x,y],jpl[:,x,y]) 
     corrcoefMatrix_gfzjpl[x][y] = temp[0,1] 

corrcoefMatrix_gfzjpl = np.squeeze(np.asarray(corrcoefMatrix_gfzjpl)) 

基本上我由含有零的矩陣,並與來自corrcoef矩陣的相關係數的值替換它們。我爲每個網格單元做了這樣的操作,通過lats和lons每個單元的for循環。 然後我創建了一個新的netcdf文件,定義了所有的尺寸和變量。

nc_data.createDimension('lat', len(lats)) 
nc_data.createDimension('lon', len(lons)) 
nc_data.createDimension('data', 1) 

latitudes = nc_data.createVariable('latitude', 'f4', ('lat',)) 
longitudes = nc_data.createVariable('longitude', 'f4', ('lon',)) 
corrcoef_gfzjpl = nc_data.createVariable('corrcoef_gfzjpl', 'f4', ('lat', 'lon', 'data'), fill_value=-999.0) 

latitudes.units = 'degree_north' 
longitudes.units = 'degree_east' 
latitudes[:] = np.arange(-90, 90, 180./len(lats)) 
longitudes[:] = np.arange(0, 360, 360./len(lons)) 
corrcoef_gfzjpl[:,:] = corrcoefMatrix_gfzjpl[:,:] 

改進建議非常受歡迎!

+0

或者,您可以將係數矩陣初始化爲'corrcoefMatrix_gfzjpl = np.zeros([len(lats),len(lons)])'',然後將數據存儲在double for-loop中,作爲corrcoefMatrix_gfzjpl [x,y] =溫度[0,1]'。 – N1B4