2015-09-27 107 views
1

當我從netCDF文件中提取數據再分析(可變壓力(SLP),01/01/2014)數據是非常高分辨率(9公里網格),這使得產生的圖像相當嘈雜。我想將數據放入較低分辨率的網格(例如1度)。我試圖使用函數meshgrid和gridata,但經驗不足無法使它工作。有誰知道如何解決?謝謝。高分辨率再分析數據

from netCDF4 import Dataset 
    import numpy as np 
    from scipy.interpolate import griddata 
    file = Dataset('slp_2014_01_01.nc', 'r') 
    # Printing variables 
    print ' ' 
    print ' ' 
    print '----------------------------------------------------------' 
    for i,variable in enumerate(file.variables): 
     print ' '+str(i),variable 
     if i == 2: 
      current_variable = variable 
    print ' ' 
    print 'Variable: ', current_variable.upper() 
    print 'File name: ', file_name 
    lat = file.variables['lat'][:] 
    lon = file.variables['lon'][:] 
    slp = file.variables['slp'][:] 

    lon_i = np.linspace(lon[0], lon[len(REANALYSIS_lon)-1], num=len(lon)*2, endpoint=True, retstep=False)  
    lat_i = np.linspace(lat[0], lat[len(lat)-1], num=len(lat)*2, endpoint=True, retstep=False) 

    lon_grid, lat_grid = np.meshgrid(lon_i,lat_i) 

    temp_slp = np.asarray(slp).squeeze() 
    new_slp = temp_slp.reshape(temp_slp.size) 

    slp_grid = griddata((lon, lat), new_slp, (lon_grid, lat_grid),method='cubic') 

正如我所說,我試圖用meshgrid和DataGrid的功能,但將產生以下錯誤:

回溯(最近通話最後一個): 文件 「REANALYSIS_LOCAL.py」,346行,在
LON,經緯度,時間,VAR,variavel_atual = netCDF_builder_local(caminho_netcdf_local,nome_arquivo,DT) 文件 「REANALYSIS_LOCAL.py」,管線143,在netCDF_builder_local
slp_grid =的GridData((LON,LAT),new_slp,(lon_grid ,lat_grid),method ='cubic')
File 「/home/carlos/anaconda/lib/python2.7/site-packages/scipy/interpolate/ndgriddata.py」,第182行,在griddata中 points = _ndim_coords_from_arrays(points)
文件「interpnd.pyx」,行176 ,在scipy.interpolate.interpnd._ndim_coords_from_arrays(scipy/interpolate/interpnd.c:4064)
文件「/home/carlos/anaconda/lib/python2.7/site-packages/numpy/lib/stride_tricks.py」,第101行,在broadcast_arrays 「軸%r上的不兼容維度」。 %(軸))
ValueError異常:形狀不匹配:兩個或更多個陣列具有在軸線上不相容的維度0

變量的尺寸爲:)(144,
緯度::
LON(73)
lon_i:(288,)
lat_i:(146,)
lon_grid:(146,288)
lat_grid:(146,288)
new_slp:(10512)

在new_slp的值是: new_slp:[102485. 102485. 102485. ...,100710. 100710. 100710.]

因爲再分析分辨率較高,所以目的是增加變量(lon,lat和slp)中的值。然後,該決議可能是最詳細的(更多點)。

例如:變量LAT有幾點:

原始尺寸可變緯度:(73)
緯度:[90 87.5 82.5 85 77.5 80 72.5 75 67.5 70 62.5 65 60. 57.5 55.52.5 50. 47.5 45. 42.5 40. 37.5 35. 32.5 30. 27.5 25. 22.5 20. 17.5 15. 12.5 10. 7.5 5. 2.5 0. -2.5 -5。 -7.5 -10。 -12.5 -15。 -17.5 -20。 -22.5 -25。 -27.5 -30。 -32.5 -35。 -37.5 -40。 -42.5 -45。 -47.5 -50。 -52.5 -55。 -57.5 -60。 -62.5 -65。 -67.5 -70。 -72.5 -75。 -77.5 -80。 -82.5 -85。 -87.5 -90。當我定義代碼行時:lat_i = np.linspace(lat [0],lat [len(lat)-1],num = len(lat)* 2,endpoint = True,retstep = False)[0124]我將lat變量la_i的值加倍(146)

lat_i:[90. 88.75862069 87.51724138 86.27586207 85.03448276 83.79310345 82.55172414 81.31034483 80.06896552 78。82758621 77.5862069
...
-78.82758621 -80.06896552 -81.31034483 -82.55172414 -83.79310345 -85.03448276 -86.27586207 -87.51724138 -88.75862069 -90。 ]

我需要的想法是相同的,在此代碼中可用,其中x是lon,y是lat,slp是z。
從scipy.interpolate進口的GridData 進口numpy的爲NP 進口matplotlib.pyplot如PLT

x=np.linspace(1.,10.,20) 
y=np.linspace(1.,10.,20) 
z=z = np.random.random(20) 
xi=np.linspace(1.,10.,40) 
yi=np.linspace(1.,10.,40) 

X,Y= np.meshgrid(xi,yi) 

Z = griddata((x, y), z, (X, Y),method='nearest') 

plt.contourf(X,Y,Z) 

回答

3

取決於你的最終目的,你可以使用CDO來REGRID整個文件

cdo remapbil,r360x180 infile outfile 

或就像這樣從原始文件中繪製每秒或第三個值:

plt.pcolormesh(lon[::2,::2],lat[::2,::2],var1[::2,::2]) 

您顯示的錯誤消息只是說維度不多,只是在出現錯誤之前打印變量的形狀並嘗試使其工作。

爲什麼你的代碼不起作用? 您選擇的方法要求輸入座標爲數據點的lon,lat對,而不是網格座標。如果您的數據點的形狀爲10000,則您的座標必須爲形狀(10000,2),而不是(100,100)。 但是,由於griddata是爲非結構化數據,它不會有效爲您的目的,我建議使用類似scipy.interpolate.RegularGridInterpolator

但無論如何,如果您需要使用插值數據不止一次,我建議用cdo創建新的netCDF文件並處理它們,而不是每次運行腳本時插入數據。

+0

感謝您的回答......但我添加了更多關於我的目的/我需要做什麼的細節。 –

+0

好吧,先說你想要1度分辨率而不是9公里?但是,然後你創建新的座標(lon_i,lat_i),它是初始數組的兩倍,從而獲得「更精細」的分辨率?!你仍然可以用cdo來做到這一點,它正是爲了這個目的而開發的。 – kakk11

+0

我編輯了我的第一個答案,以解釋爲什麼你的代碼看起來不起作用。 – kakk11

0

感謝您的幫助。真的,我的問題是關於維度。我正在學習使用海洋學數據。所以,我解決了這個代碼的問題。

lonbounds = [25,59] 
latbounds = [-10,-33] 

#longitude lower and upper index 
lonli = np.argmin(np.abs(lon - lonbounds[0])) 
lonui = np.argmin(np.abs(lon - lonbounds[1])) 

#latitude lower and upper index 
latli = np.argmin(np.abs(lat - latbounds[0])) 
latui = np.argmin(np.abs(lat - latbounds[1])) 

#limiting of the interest region/data 
lon_f = file.variables['lon'][lonli:lonui] 
lat_f = file.variables['lat'][latli:latui] 
slp_f = file.variables['slp'][0,latli:latui,lonli:lonui] 

#creating a matrix with the filtered data (area to be searched) for use in gridData function of python 
lon_f_grid, lat_f_grid = np.meshgrid(lon_f,lat_f) 

#adjusting the data (size 1) for use in gridData function of python 
lon_f1 = lon_f_grid.reshape(lon_f_grid.size) 
lat_f1 = lat_f_grid.reshape(lat_f_grid.size) 
slp_f1 = slp_f.reshape(slp_f.size) 

#increasing the resolution of data (1000 points) of longitude and latitude for the data to be more refined 
lon_r = np.linspace(lon_f[0], lon_f[len(lon_f)-1], num=1000, endpoint=True, retstep=False) 
lat_r = np.linspace(lat_f[0], lat_f[len(lat_f)-1], num=1000, endpoint=True, retstep=False) 

#creating a matrix with the filtered data (area to be searched) and higher resolution for use in gridData function of python 
lon_r_grid, lat_r_grid = np.meshgrid(lon_r,lat_r) 

#applying gridata that can be generated since pressure (SLP) with higher resolution. 
slp_r = griddata((lon_f1,lat_f1),slp_f1,(lon_r_grid,lat_r_grid),method='cubic') 

擁抱, 卡洛斯。