動作環境
Xeon E5-2620 v4 (8コア) x 2
32GB RAM
CentOS 6.8 (64bit)
openmpi-1.8.x86_64 とその-devel
mpich.x86_64 3.1-5.el6とその-devel
gcc version 4.4.7 (とgfortran)
NCAR Command Language Version 6.3.0
WRF v3.7.1を使用。
Python 2.6.6 (r266:84292, Aug 18 2016, 15:13:37)
Python 3.6.0 on virtualenv
関連 http://qiita.com/7of9/items/6fd5d52439b0adc92d95
LN-test_nc: WRF関連の処理結果ファイル
LN-latlon_nc: 緯度経度情報ファイル
2つのファイルを読込んで、(lat, lon, 対象値)の形式とするために以下を実装していた。
- proc_vallatlon_170315_exec: bash script
- read_lat_lon.ncl : NCAR Command Language script
- combine_3files_170315.py: Python script
NCLはよく分かっていなく、bashは苦手。Pythonもまだまだ学習中。
これら3つのファイルを使い続けるのはメンテナンス性が悪い。
PythonでnetCDF4パッケージを使えば簡単、と聞いてそのようにした。
参考 http://qiita.com/okadate/items/954574a95545b06ca257
code
read_netcdf_170323.py
import netCDF4
import sys
# on Python 2.6.6
# codingrule:PEP8
# LN-test_nc: netCDF file including processed data
# LN-latlon_nc: netCDF file including (lat, lon)
#
# 1. read data
nc = netCDF4.Dataset('LN-test_nc', 'r')
cnc = nc.variables['CONC'][:]
nc.close()
# 2. read lat lon
nc = netCDF4.Dataset('LN-latlon_nc', 'r')
lat = nc.variables['XLAT'][:]
lon = nc.variables['XLONG'][:]
nc.close()
# debug
print(cnc.shape) # (24, 1, 1, 3, 65, 82)
print(lat.shape) # (65, 82)
print(lon.shape) # (65, 82)
print(cnc[0, 0, 0, 0, 0, 81]) # (24, 1, 1, 3, 65, 82)
# 3. extract map
timidx = 2 # 2: arbitrary
cmap = cnc[timidx, 0, 0, 0, :, :] # (65, 82)
print(cmap.shape)
# 4. output
for alat, alon, acnc in zip(lat, lon, cmap):
for idx in range(len(alat)):
print('%.5f, %.5f, %.5f' % (alat[idx], alon[idx], acnc[idx]))
# debug
if idx > 5:
sys.exit() # debug
結果
$ python read_netcdf_170323.py
(24, 1, 1, 3, 65, 82)
(65, 82)
(65, 82)
0.0
(65, 82)
36.07000, 111.17000, 0.00000
36.07000, 111.29000, 0.00000
36.07000, 111.41000, 0.00000
36.07000, 111.53000, 0.00000
36.07000, 111.65000, 0.00000
36.07000, 111.77000, 0.00000
36.07000, 111.89000, 0.00000
これで1つのファイルだけメンテナンスすればよくなった。
NCL(NCAR Command Language script)は既存スクリプトを読めればいい程度。書くための学習コストは不要となった。
上記コードはvirtualenvのPython 3.6.0でも動作する。
(easy_install netCDF4した)
教えていただいた事項
@shiracamus さんのコメントにてenumerate()
とzip(*)
の使用例を教えていただきました。
情報感謝です。