TES水同位体データの解析

TESによる水同位体データの解析についてまとめる。 TESデータの読み込み、水同位体比の計算、そして描画までについて。

同位体NICAMの気候値と比較するのが目的であるため、TES L3の月データを使用。

Import Moduels

In [75]:
import h5py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as colors
from mpl_toolkits.basemap import Basemap

Read

In [99]:
f1 = h5py.File('./TES-Aura_L3-H2O-M2007m07_F01_11.he5', 'r') 
f2 = h5py.File('./TES-Aura_L3-HDO-M2007m07_F01_11.he5', 'r')

読み込んだデータは階層構造になっている。 データの階層構造がどのようになっているかは、下記のように確認できる。

In [3]:
for k in f:
    print (k)

HDFEOS
HDFEOS INFORMATION
In [119]:
for k in f1['HDFEOS']['GRIDS']['NadirGrid']['Data Fields']:
          print (k)

H2O
H2ODataCount
H2OMaximum
H2OMinimum
H2OStdDeviation
Latitude
Longitude
Pressure
TotColDensDataCount
TotColDensMaximum
TotColDensMinimum
TotColDensStdDeviation
TotalColumnDensity
In [120]:
for k in f2['HDFEOS']['GRIDS']['NadirGrid']['Data Fields']:
          print (k)

HDO
HDODataCount
HDOMaximum
HDOMinimum
HDOStdDeviation
Latitude
Longitude
Pressure
TotColDensDataCount
TotColDensMaximum
TotColDensMinimum
TotColDensStdDeviation
TotalColumnDensity

H2OやHDOのデータは下記のように取得できる。

In [122]:
h2o = f1['HDFEOS']['GRIDS']['NadirGrid']['Data Fields']['H2O'].value
hdo = f2['HDFEOS']['GRIDS']['NadirGrid']['Data Fields']['HDO'].value

取得したデータの構造は下記の通り。 左から高さ、緯度、経度の順番になっている。欠損値は-999.0

In [123]:
print h2o.shape
print np.max(h2o), np.min(hdo)

(15, 83, 90)
0.0535723 -999.0

緯度、経度、気圧のデータも格納されているの取得可能。 それぞれのデータについても表示する。

In [124]:
lon = f1['HDFEOS']['GRIDS']['NadirGrid']['Data Fields']['Longitude'].value
lat = f1['HDFEOS']['GRIDS']['NadirGrid']['Data Fields']['Latitude'].value
pres = f1['HDFEOS']['GRIDS']['NadirGrid']['Data Fields']['Pressure'].value
In [125]:
print (pres)
print lon
print lat

[ 825.40197754  681.29101562  464.16000366  316.22698975  215.44400024
  146.77900696  100.           68.12950134   46.41579819   31.62290001
   21.54430008   14.67800045   10.            6.81291008    4.64160013]
[-180. -176. -172. -168. -164. -160. -156. -152. -148. -144. -140. -136.
 -132. -128. -124. -120. -116. -112. -108. -104. -100.  -96.  -92.  -88.
  -84.  -80.  -76.  -72.  -68.  -64.  -60.  -56.  -52.  -48.  -44.  -40.
  -36.  -32.  -28.  -24.  -20.  -16.  -12.   -8.   -4.    0.    4.    8.
   12.   16.   20.   24.   28.   32.   36.   40.   44.   48.   52.   56.
   60.   64.   68.   72.   76.   80.   84.   88.   92.   96.  100.  104.
  108.  112.  116.  120.  124.  128.  132.  136.  140.  144.  148.  152.
  156.  160.  164.  168.  172.  176.]
[-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.]

Calculate vapor isotope

TES観測による水素同位体比を計算する。 SMOWの値が必要で、Wordan et al. (2007). Natureによると、ここでの値はVSMOW=3.1110-4とのこと。 基本的にはδ=（R/VSMOW-1）1000で算定できる。

In [126]:
# 欠損値の処理
h2o[np.where(h2o==-999.)] = np.nan
hdo[np.where(hdo==-999.)] = np.nan
In [127]:
# Wordan2007では、850-550hPaにおける値を使用していたので、ここでも同様に算定。
ahdo = np.nanmean(hdo[0:2,:,:],axis=0)
ah2o = np.nanmean(h2o[0:2,:,:],axis=0)
iso1=ahdo/ah2o
iso1 = (iso1/(3.11*0.0001)-1)*1000
iso1[np.where(ah2o==0.)]=np.nan
print np.nanmax(iso1), np.nanmin(iso1)

247.598 -358.836

/Users/masatano/.pyenv/versions/anaconda3-5.0.1/envs/py2/lib/python2.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: Mean of empty slice
  
/Users/masatano/.pyenv/versions/anaconda3-5.0.1/envs/py2/lib/python2.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: Mean of empty slice
  This is separate from the ipykernel package so we can avoid doing imports until

Draw map

In [129]:
LonMin, LonMax = np.min(lon), np.max(lon)
LatMin, LatMax = np.min(lat), np.max(lat)

plt.rcParams['figure.figsize']  =  (12.0, 12.0)
fig = plt.figure(figsize = (8, 6))

setcmap = cm.rainbow
setcmap.set_under('#aaaaaa')
setcmap.set_over('#ff0000')
bounds = np.array([-250, -200, -150, -100, -50, 0, 50])

m=Basemap(projection='cyl',resolution='l',llcrnrlon=LonMin,llcrnrlat=LatMin,urcrnrlon=LonMax,urcrnrlat=LatMax)
m.drawcoastlines(color='k',linewidth=0.4)
cs = m.imshow(iso1,interpolation='nearest',cmap=setcmap,norm=colors.BoundaryNorm(bounds, 256))

cax = fig.add_axes([0.1,0.05,0.8,0.025])
cb = plt.colorbar(cs,cax,orientation='horizontal',extend='both')
cb.ax.tick_params(labelsize=14)

plt.show()
In [118]:
dat=iso1.reshape(-1)
plt.hist(iso1,normed=True)
plt.show()