トップページに戻る

ファイル読み込み

目次

  • はじめに
  • テキスト(csv)ファイル読み込み
  • NetCDFファイル読み込み
  • Fortranバイナリファイル読み込み

はじめに

データ解析を行う際に、データファイルを読み込むというのは最も基本的な操作になります。
地球物理学の分野でメジャーなデータファイルの形式は、以下の3つだと思います。

  1. テキストファイル
  2. NetCDFファイル
  3. Fortranバイナリファイル

Pythonを用いれば、この3つのどの形式のデータファイルも簡単に読み込むことが可能です。
以下でそれぞれ順を追って説明していきますが、ご自身の必要に応じて2.、3.から読み進める等していただいても構いません。

テキスト(csv)ファイル読み込み

カンマで区切られているテキストファイルを読み込んでみましょう。
テキストファイル内の数値データは、numpyのloadtxt関数で読み込むことが出来ます。
なお、今回読み込むファイルは、気象庁のページから取得した、エルニーニョ指数データです。

In [1]:
import numpy as np
data=np.loadtxt('NINO.3.csv',comments='year',delimiter=',',dtype='float')
  • commentsで、読み飛ばす行の左端に存在する文字列を指定する。
  • delimiterで、区切り文字を指定する。もしスペースで区切ってあった場合には、delimiter=...という記述は必要ない。
  • dtypeで、データをどの形式で読み込むかを指定する。デフォルトはfloat(浮動小数点数)。整数として読み込みたい場合にはintとすればよい。
In [2]:
data
Out[2]:
array([[  2.00100000e+03,  -4.00000000e-01,  -3.00000000e-01,
         -2.00000000e-01,  -1.00000000e-01,   0.00000000e+00,
          0.00000000e+00,  -1.00000000e-01,  -2.00000000e-01,
         -3.00000000e-01,  -4.00000000e-01,  -5.00000000e-01,
         -4.00000000e-01],
       [  2.00200000e+03,  -3.00000000e-01,  -1.00000000e-01,
          0.00000000e+00,   3.00000000e-01,   4.00000000e-01,
          5.00000000e-01,   6.00000000e-01,   7.00000000e-01,
          8.00000000e-01,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00],
       [  2.00300000e+03,   8.00000000e-01,   5.00000000e-01,
          0.00000000e+00,  -2.00000000e-01,  -2.00000000e-01,
         -2.00000000e-01,  -1.00000000e-01,   1.00000000e-01,
          3.00000000e-01,   4.00000000e-01,   4.00000000e-01,
          4.00000000e-01],
       [  2.00400000e+03,   4.00000000e-01,   3.00000000e-01,
          1.00000000e-01,   0.00000000e+00,  -1.00000000e-01,
          0.00000000e+00,   0.00000000e+00,   2.00000000e-01,
          3.00000000e-01,   4.00000000e-01,   4.00000000e-01,
          3.00000000e-01],
       [  2.00500000e+03,   2.00000000e-01,   1.00000000e-01,
          1.00000000e-01,   1.00000000e-01,   2.00000000e-01,
          3.00000000e-01,   2.00000000e-01,   1.00000000e-01,
         -2.00000000e-01,  -5.00000000e-01,  -7.00000000e-01,
         -8.00000000e-01],
       [  2.00600000e+03,  -8.00000000e-01,  -7.00000000e-01,
         -5.00000000e-01,  -3.00000000e-01,  -2.00000000e-01,
          1.00000000e-01,   3.00000000e-01,   4.00000000e-01,
          6.00000000e-01,   8.00000000e-01,   9.00000000e-01,
          8.00000000e-01],
       [  2.00700000e+03,   5.00000000e-01,   2.00000000e-01,
         -2.00000000e-01,  -5.00000000e-01,  -6.00000000e-01,
         -7.00000000e-01,  -9.00000000e-01,  -1.10000000e+00,
         -1.30000000e+00,  -1.40000000e+00,  -1.50000000e+00,
         -1.50000000e+00],
       [  2.00800000e+03,  -1.40000000e+00,  -1.10000000e+00,
         -8.00000000e-01,  -5.00000000e-01,  -1.00000000e-01,
          1.00000000e-01,   2.00000000e-01,   2.00000000e-01,
          2.00000000e-01,   0.00000000e+00,  -2.00000000e-01,
         -4.00000000e-01],
       [  2.00900000e+03,  -5.00000000e-01,  -5.00000000e-01,
         -3.00000000e-01,   0.00000000e+00,   3.00000000e-01,
          5.00000000e-01,   7.00000000e-01,   8.00000000e-01,
          9.00000000e-01,   1.00000000e+00,   1.00000000e+00,
          1.10000000e+00],
       [  2.01000000e+03,   1.10000000e+00,   9.00000000e-01,
          7.00000000e-01,   3.00000000e-01,   0.00000000e+00,
         -4.00000000e-01,  -8.00000000e-01,  -1.10000000e+00,
         -1.30000000e+00,  -1.40000000e+00,  -1.50000000e+00,
         -1.40000000e+00],
       [  2.01100000e+03,  -1.20000000e+00,  -9.00000000e-01,
         -7.00000000e-01,  -4.00000000e-01,  -2.00000000e-01,
         -2.00000000e-01,  -2.00000000e-01,  -4.00000000e-01,
         -6.00000000e-01,  -8.00000000e-01,  -8.00000000e-01,
         -7.00000000e-01],
       [  2.01200000e+03,  -6.00000000e-01,  -3.00000000e-01,
         -1.00000000e-01,   1.00000000e-01,   3.00000000e-01,
          4.00000000e-01,   5.00000000e-01,   5.00000000e-01,
          4.00000000e-01,   2.00000000e-01,   0.00000000e+00,
         -2.00000000e-01],
       [  2.01300000e+03,  -2.00000000e-01,  -3.00000000e-01,
         -4.00000000e-01,  -4.00000000e-01,  -5.00000000e-01,
         -6.00000000e-01,  -6.00000000e-01,  -5.00000000e-01,
         -3.00000000e-01,  -2.00000000e-01,  -1.00000000e-01,
         -2.00000000e-01],
       [  2.01400000e+03,  -2.00000000e-01,  -1.00000000e-01,
          0.00000000e+00,   2.00000000e-01,   4.00000000e-01,
          5.00000000e-01,   5.00000000e-01,   5.00000000e-01,
          6.00000000e-01,   7.00000000e-01,   7.00000000e-01,
          6.00000000e-01],
       [  2.01500000e+03,   5.00000000e-01,   5.00000000e-01,
          6.00000000e-01,   8.00000000e-01,   1.20000000e+00,
          1.60000000e+00,   1.90000000e+00,   2.20000000e+00,
          2.50000000e+00,   2.70000000e+00,   2.80000000e+00,
          2.70000000e+00]])

NetCDFファイル読み込み

netCDFファイルの読み込みには、netCDF4モジュールを用います。
netCDF4モジュールが使えるようにする方法に関してはこちらをご参照ください。
試しに、NOAAのサイトからダウンロードした2017年2月の全球海面水温データを読み込んでみることにします。
まず、ターミナル上で
$ ncdump -c ersst.v4.201702.nc
と実行して、当該netCDFファイル内のデータの構造に目を通しておくことにします。

In [3]:
# Jupyter上では、!を先頭につけることでターミナルのコマンドを実行することが出来る。
!ncdump -c ersst.v4.201702.nc
netcdf ersst.v4.201702 {
dimensions:
	lat = 89 ;
	lev = 1 ;
	lon = 180 ;
	time = 1 ;

variables:
	double lat(lat) ;
		lat:units = "degrees_north" ;
		lat:long_name = "Latitude" ;
		lat:standard_name = "latitude" ;
		lat:axis = "Y" ;
		lat:bounds = "lat_bnds" ;
		lat:grids = "Uniform grid from -88 to 88 by 2" ;
	double lev(lev) ;
		lev:units = "meters" ;
		lev:long_name = "Depth of sea surface temperature measurements" ;
		lev:standard_name = "depth" ;
		lev:axis = "Z" ;
		lev:positive = "down" ;
		lev:_CoordinateAxisType = "Height" ;
		lev:comment = "Measurement depth of in situ sea surface temperature varies" ;
	double lon(lon) ;
		lon:units = "degrees_east" ;
		lon:long_name = "Longitude" ;
		lon:standard_name = "longitude" ;
		lon:axis = "X" ;
		lon:bounds = "lon_bnds" ;
		lon:grids = "Uniform grid from 0 to 358 by 2" ;
	float sst(time, lev, lat, lon) ;
		sst:_FillValue = -999.f ;
		sst:long_name = "Extended reconstructed sea surface temperature" ;
		sst:standard_name = "SST" ;
		sst:units = "degree_C" ;
		sst:add_offset = 0.f ;
		sst:scale_factor = 1.f ;
		sst:valid_min = -3.f ;
		sst:valid_max = 45.f ;
	double time(time) ;
		time:long_name = "Center time of the day" ;
		time:standard_name = "time" ;
		time:axis = "T" ;
		time:delta_t = "0000-01-00" ;
		time:avg_period = "0000-01-00" ;
		time:units = "minutes since 2017-02-01 00:00" ;
	float ssta(time, lev, lat, lon) ;
		ssta:_FillValue = -999.f ;
		ssta:long_name = "Extended reconstructed SST anomalies" ;
		ssta:standard_name = "SSTA" ;
		ssta:units = "degree_C" ;
		ssta:add_offset = 0.f ;
		ssta:scale_factor = 1.f ;
		ssta:valid_min = -12.f ;
		ssta:valid_max = 12.f ;

// global attributes:
		:Conventions = "CF-1.6" ;
		:Metadata_Conventions = "CF-1.6, Unidata Dataset Discovery v1.0" ;
		:metadata_link = "C00884" ;
		:id = "ersst.v4.201702" ;
		:naming_authority = "gov.noaa.ncdc" ;
		:title = "NOAA ERSSTv4 (in situ only)" ;
		:summary = "ERSST.v4 is developped based on v3b after revisions of 11 parameters using updated data sets and advanced knowledge of ERSST analysis" ;
		:institution = "NOAA/NESDIS/NCEI/CCOG" ;
		:creator_name = "Boyin Huang" ;
		:creator_email = "boyin.huang@noaa.gov" ;
		:date_modified = "2017/03/03" ;
		:production_version = "Version 5" ;
		:history = "Fri Mar  3 13:35:18 2017: ncatted -O -a _FillValue,ssta,o,f,-999.0 ssta.nc\n",
    "Version 5 based on Version 4" ;
		:publisher_name = "Boyin Huang" ;
		:publisher_email = "boyin.huang@noaa.gov" ;
		:creator_url = "http://www.ncdc.noaa.gov" ;
		:license = "No constraints on data access or use" ;
		:time_coverage_start = "2017-02-15T000000Z" ;
		:time_coverage_end = "2017-02-15T000000Z" ;
		:geospatial_lon_min = -1.f ;
		:geospatial_lon_max = 359.f ;
		:geospatial_lat_min = -89.f ;
		:geospatial_lat_max = 89.f ;
		:geospatial_lat_units = "degrees_north" ;
		:geospatial_lat_resolution = 2.f ;
		:geospatial_lon_units = "degrees_east" ;
		:geospatial_lon_resolution = 2.f ;
		:spatial_resolution = "2.0 degree grid" ;
		:cdm_data_type = "Grid" ;
		:processing_level = "L4" ;
		:standard_name_vocabulary = "CF Standard Name Table v27" ;
		:keywords = "Earth Science > Oceans > Ocean Temperature > Sea Sur face Temperature &gt" ;
		:keywords_vocabulary = "NASA Global Change Master Directory (GCMD) Science Keywords" ;
		:project = "NOAA Extended Reconstructed Sea Surface Temperature (ERSST)" ;
		:platform = "Ship and Buoy SSTs from ICOADS R2.5 and NCEP GTS" ;
		:instrument = "Conventional thermometers" ;
		:source = "ICOADS R2.5 SST, NCEP GTS SST, HadISST ice, NCEP ice" ;
		:comment = "SSTs were observed by conventional thermometers in Buckets (in sulated or un-insulated canvas and wooded buckets), Engine Room Intakers, or floats and drifters" ;
		:references = "Huang et al, 2015: Extended Reconstructed Sea Surface Temperatures Version 4 (ERSST.v4), Part I. Upgrades and Intercomparisons. Journal of Climate, DOI: 10.1175/JCLI-D-14-00006.1." ;
		:climatology = "Climatology is based on 1971-2000 SST, Xue, Y., T. M. Smith, and R. W. Reynolds, 2003: Interdecadal changes of 30-yr SST normals during 1871.2000. Journal of Climate, 16, 1601-1612." ;
		:description = "In situ data: ICOADS2.5 before 2007, NCEP in situ data from 2008 to present. Ice data: HadISST ice before 2015 and NCEP ice after 2015." ;

data:

 lat = -88, -86, -84, -82, -80, -78, -76, -74, -72, -70, -68, -66, -64, -62, 
    -60, -58, -56, -54, -52, -50, -48, -46, -44, -42, -40, -38, -36, -34, 
    -32, -30, -28, -26, -24, -22, -20, -18, -16, -14, -12, -10, -8, -6, -4, 
    -2, 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 
    36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 
    72, 74, 76, 78, 80, 82, 84, 86, 88 ;

 lev = 0 ;

 lon = 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 
    38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 
    74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100, 102, 104, 106, 
    108, 110, 112, 114, 116, 118, 120, 122, 124, 126, 128, 130, 132, 134, 
    136, 138, 140, 142, 144, 146, 148, 150, 152, 154, 156, 158, 160, 162, 
    164, 166, 168, 170, 172, 174, 176, 178, 180, 182, 184, 186, 188, 190, 
    192, 194, 196, 198, 200, 202, 204, 206, 208, 210, 212, 214, 216, 218, 
    220, 222, 224, 226, 228, 230, 232, 234, 236, 238, 240, 242, 244, 246, 
    248, 250, 252, 254, 256, 258, 260, 262, 264, 266, 268, 270, 272, 274, 
    276, 278, 280, 282, 284, 286, 288, 290, 292, 294, 296, 298, 300, 302, 
    304, 306, 308, 310, 312, 314, 316, 318, 320, 322, 324, 326, 328, 330, 
    332, 334, 336, 338, 340, 342, 344, 346, 348, 350, 352, 354, 356, 358 ;

 time = 0 ;
}

以上の通り、海面水温データは、'sst'という変数名で格納されているということがわかりました。
そこで、この海面水温データをnetCDF4モジュールを用いて取り出してみることにします。

In [4]:
import netCDF4
nc=netCDF4.Dataset('ersst.v4.201702.nc','r')
data=nc.variables['sst'][:]

たったこれだけです。
2行目で、netCDFファイルの中身の情報を持った「netCDF4オブジェクト」を生成し、
3行目で、そのオブジェクトの一要素(今回の場合'sst'という変数名で格納されたデータの中身)を取り出しているわけです。

In [5]:
print(type(nc))
print(data)
<class 'netCDF4._netCDF4.Dataset'>
[[[[-- -- -- ..., -- -- --]
   [-- -- -- ..., -- -- --]
   [-- -- -- ..., -- -- --]
   ..., 
   [-1.5788253545761108 -1.5387723445892334 -1.518031358718872 ...,
    -1.7999999523162842 -1.7702364921569824 -1.7328710556030273]
   [-1.5536582469940186 -1.5258724689483643 -1.5181310176849365 ...,
    -1.6532748937606812 -1.6176127195358276 -1.5812113285064697]
   [-1.7319483757019043 -1.7999999523162842 -1.7999999523162842 ...,
    -1.741258144378662 -1.7337182760238647 -1.729689598083496]]]]

Fortranバイナリファイル読み込み

ここでは、まず
リトルエンディアンのヘッダなし4バイト浮動小数点バイナリ(俗に言うGrADS形式と呼ばれるもの)
の中身を読み込んでみることにします。

In [6]:
!cat write_binary.f90 
program main

  implicit none

  integer,parameter::N=10,M=20
  integer::i,j
  real,dimension(1:N,1:M)::x

  open(10,file='filename.out',form='unformatted',access='direct',recl=N*4)

  do i = 1,N
     do j = 1,M
        x(i,j) = i+j*2
     end do
  end do

  do j = 1,M
     write(10,rec=j)(x(i,j),i=1,N)
  end do

  close(10)

end program main

まず、fortranプログラムを実行して架空のデータを作成。

In [7]:
!gfortran write_binary.f90
In [8]:
!./a.out

続いてpythonでファイルの読み込みです。

In [9]:
import numpy as np
N=10  #1レコード番号あたりに格納されているデータの数。
M=20  #レコードの総数。
f=open('filename.out','r')
dty=np.dtype([('data','<'+str(N)+'f')])
chunk=np.fromfile(f,dtype=dty,count=M)
data=np.array([chunk[j]['data'] for j in range(0,M)])
In [10]:
data
Out[10]:
array([[  3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,  12.],
       [  5.,   6.,   7.,   8.,   9.,  10.,  11.,  12.,  13.,  14.],
       [  7.,   8.,   9.,  10.,  11.,  12.,  13.,  14.,  15.,  16.],
       [  9.,  10.,  11.,  12.,  13.,  14.,  15.,  16.,  17.,  18.],
       [ 11.,  12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.],
       [ 13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,  22.],
       [ 15.,  16.,  17.,  18.,  19.,  20.,  21.,  22.,  23.,  24.],
       [ 17.,  18.,  19.,  20.,  21.,  22.,  23.,  24.,  25.,  26.],
       [ 19.,  20.,  21.,  22.,  23.,  24.,  25.,  26.,  27.,  28.],
       [ 21.,  22.,  23.,  24.,  25.,  26.,  27.,  28.,  29.,  30.],
       [ 23.,  24.,  25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.],
       [ 25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.,  33.,  34.],
       [ 27.,  28.,  29.,  30.,  31.,  32.,  33.,  34.,  35.,  36.],
       [ 29.,  30.,  31.,  32.,  33.,  34.,  35.,  36.,  37.,  38.],
       [ 31.,  32.,  33.,  34.,  35.,  36.,  37.,  38.,  39.,  40.],
       [ 33.,  34.,  35.,  36.,  37.,  38.,  39.,  40.,  41.,  42.],
       [ 35.,  36.,  37.,  38.,  39.,  40.,  41.,  42.,  43.,  44.],
       [ 37.,  38.,  39.,  40.,  41.,  42.,  43.,  44.,  45.,  46.],
       [ 39.,  40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,  48.],
       [ 41.,  42.,  43.,  44.,  45.,  46.,  47.,  48.,  49.,  50.]], dtype=float32)

最終行は、

In [11]:
for j in range(0,M):  
     data[j,:]=chunk[j]['data']

を1行で書き換えたものです。
chunk[k-1]が、fortranで言うところのレコード番号kのデータに対応しています。
なので、例えばレコード番号6のデータだけを取り出したい場合には、最終行を

In [12]:
data=chunk[5]['data']  

に置き換えれば良いです。

In [13]:
data
Out[13]:
array([ 13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,  22.], dtype=float32)

誤字やおかしい点などがあったら @zawawahoge (Twitter) にお気軽にご連絡ください。

トップページに戻る