データ解析を行う際に、データファイルを読み込むというのは最も基本的な操作になります。
地球物理学の分野でメジャーなデータファイルの形式は、以下の3つだと思います。
Pythonを用いれば、この3つのどの形式のデータファイルも簡単に読み込むことが可能です。
以下でそれぞれ順を追って説明していきますが、ご自身の必要に応じて2.、3.から読み進める等していただいても構いません。
import numpy as np
data=np.loadtxt('NINO.3.csv',comments='year',delimiter=',',dtype='float')
data
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 >" ;
: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モジュールを用いて取り出してみることにします。
import netCDF4
nc=netCDF4.Dataset('ersst.v4.201702.nc','r')
data=nc.variables['sst'][:]
たったこれだけです。
2行目で、netCDFファイルの中身の情報を持った「netCDF4オブジェクト」を生成し、
3行目で、そのオブジェクトの一要素(今回の場合'sst'という変数名で格納されたデータの中身)を取り出しているわけです。
print(type(nc))
print(data)
ここでは、まず
リトルエンディアンのヘッダなし4バイト浮動小数点バイナリ(俗に言うGrADS形式と呼ばれるもの)
の中身を読み込んでみることにします。
!cat write_binary.f90
まず、fortranプログラムを実行して架空のデータを作成。
!gfortran write_binary.f90
!./a.out
続いてpythonでファイルの読み込みです。
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)])
data
最終行は、
for j in range(0,M):
data[j,:]=chunk[j]['data']
を1行で書き換えたものです。
chunk[k-1]が、fortranで言うところのレコード番号kのデータに対応しています。
なので、例えばレコード番号6のデータだけを取り出したい場合には、最終行を
data=chunk[5]['data']
に置き換えれば良いです。
data