夏の異常気象をオープン・データで確認：実践！ Rで学ぶ統計解析の基礎（6）（2/2 ページ）
今夏は異様に暑かったですよね。でも、どのぐらい暑かったと思いますか？ オープン・データを利用し、Rで可視化してみましょう。
夏の異常気象：NASA GISSによるグローバル表面温度データ
さて、ここで視点をグローバルに変えて地球全体を考えてみましょう。今年の夏は地球全体をで見た場合どこでどのくらい「アノーマリー」だったのか、ということを、世界地図に空間プロットして確認してみましょう。
Rを用いてデータを地図上に空間プロットするには、Geographic Infomation System（GIS）関係の機能を実装した地図表示パッケージや空間分析パッケージが必要になります。空間プロットについては、地図情報が絡むことで普通のプロットより複雑になるせいか、まだデ・ファクト・スタンダードといえるパッケージは出ていないようです。そこで、今回はいろいろなパッケージを組み合わせて空間プロットを実現します。具体的には、spパッケージを中心に、世界地図のポリゴンはmaptoolsパッケージにあるものを利用し、実際の空間プロットにはfieldsにあるas.image.plotを利用したいと思います。
地球の地表面温度と海表面温度のアノーマリーを定期的に解析しているのが、先ほどもでてきたコロンビア大学にあるNASAのGoddard Institute for Space Studies（GISS）のSurface Temperature Analysisチームです。ちなみに、このNASAのGISSは、地球温暖化理論の父といわれるJames Hansenが長く所長を勤めていることからも分かるように、地球温暖化研究の中心的役割を果たしている機関の1つです。もともとは金星や水星の惑星探査データの解析や惑星大気のモデリングというNASAの王道的な業務をやっていたのですが、1970年代後半に入り惑星探査計画が下火になってリストラされる運命だったところを、地球温暖化研究という「ビジネスモデル」の転換をしてサバイバルしたそうです。
地球の地表面温度と海表面温度のアノーマリーデータはGISSのSurface Temperature Analysisの以下のサイトから手に入れることができます。
今回は“Mean Period”を6、7、8月の3カ月データで取り、“Smoothing Radius”を250kmにしたデータを利用します。データはここに上げておきました。
念のために記述しますが、このデータがどのようにして加工されたか、という詳細な過程を知りたい方は、データ加工用のソースコード（Fortran、C、Pythonおよびシェルスクリプト）と説明がオープンになっていますので、それを参照ください。
さて、必要なパッケージをインストールし、データを読み込むところまで話を進めましょう。パッケージのインストールに指定するCRANサイトは筑波大学を選ぶとよいでしょう。
> install.packages("sp")
> install.packages("maptools")
> install.packages("fields")
> gtemp <- read.csv("http://spreadsheets.google.com/pub?key=0AlBuJgqcP5f3dHVCYm5SZUY4cDN1VFFMeWNKQ3IxbFE&hl=en&output=csv", header = T, skip = 1, sep = ',')
> names(gtemp) <- c("i", "j", "lon", "lat", "anom")
> head(gtemp)
i j lon lat anom
1 1 1 -179 -89 0.29
2 2 1 -177 -89 0.29
3 3 1 -175 -89 0.29
4 4 1 -173 -89 0.29
5 5 1 -171 -89 0.29
6 6 1 -169 -89 0.29
> gtemp$anom[grep(9999.00, gtemp$anom, fixed=T)] <- NA
> summary(gtemp)
i j lon lat anom
Min. : 1.00 Min. : 1.0 Min. :-179.0 Min. :-89 Min. : -5.7000
1st Qu.: 45.75 1st Qu.:23.0 1st Qu.: -89.5 1st Qu.:-45 1st Qu.: 0.0200
Median : 90.50 Median :45.5 Median : 0.0 Median : 0 Median : 0.4600
Mean : 90.50 Mean :45.5 Mean : 0.0 Mean : 0 Mean : 0.5536
3rd Qu.:135.25 3rd Qu.:68.0 3rd Qu.: 89.5 3rd Qu.: 45 3rd Qu.: 0.9800
Max. :180.00 Max. :90.0 Max. : 179.0 Max. : 89 Max. : 6.0400
NA's :3729.0000
上記のコードにより、データはgtempという名前のデータフレームに格納されました。i、jはただのインデックスで、lonとlatがそれぞれ経度、緯度という地表面の座標を表します。anomが今年の夏の地表面温度／海表面温度のアノーマリーです。上記のコードの最後から2番目でやっていることには少し注意が必要です。GISSデータに含まれている“9999.00”というデータは、データが入手できないという意味ですが、このままだとRでは正しい数値として認識されてしまいます。そこで、これをRにおける“Not Avairable”の意であるNAに置き換えました。
ここまでデータを用意できたら、それを地図データに展開するべくspパッケージで空間プロット用のデータフレームに加工していきます。具体的には以下の手順です。
> library(sp)
> coordinates(gtemp) = c("lon", "lat")
> summary(gtemp)
Object of class SpatialPointsDataFrame
Coordinates:
min max
lon -179 179
lat -89 89
Is projected: NA
proj4string : [NA]
Number of points: 16200
Data attributes:
i j anom
Min. : 1.00 Min. : 1.0 Min. : -5.7000
1st Qu.: 45.75 1st Qu.:23.0 1st Qu.: 0.0200
Median : 90.50 Median :45.5 Median : 0.4600
Mean : 90.50 Mean :45.5 Mean : 0.5536
3rd Qu.:135.25 3rd Qu.:68.0 3rd Qu.: 0.9800
Max. :180.00 Max. :90.0 Max. : 6.0400
NA's :3729.0000
> gridded(gtemp) <- TRUE
> summary(gtemp)
Object of class SpatialPixelsDataFrame
Coordinates:
min max
lon -180 180
lat -90 90
Is projected: NA
proj4string : [NA]
Number of points: 16200
Grid attributes:
cellcentre.offset cellsize cells.dim
lon -179 2 180
lat -89 2 90
Data attributes:
i j anom
Min. : 1.00 Min. : 1.0 Min. : -5.7000
1st Qu.: 45.75 1st Qu.:23.0 1st Qu.: 0.0200
Median : 90.50 Median :45.5 Median : 0.4600
Mean : 90.50 Mean :45.5 Mean : 0.5536
3rd Qu.:135.25 3rd Qu.:68.0 3rd Qu.: 0.9800
Max. :180.00 Max. :90.0 Max. : 6.0400
NA's :3729.0000
上のコードですが、まず通常のデータフレームにcoodinatesで座標を指定することでSpatialPointsDataFrameに変換しました。次にgriddedを利用して全体がグリッドデータであることを明示してSpatialPixelsDataFrameにしました。
このSpatialPixelsDataFrameを実際のプロットイメージにするのが、as.image.SpatialGridDataFrameで、これを利用すると最後に経度、緯度、そして表面温度アノーマリーの値のリストになります。このリストがプロットイメージです。
> anomimage <- as.image.SpatialGridDataFrame(gtemp["anom"]) > summary(anomimage) Length Class Mode x 180 -none- numeric y 90 -none- numeric z 16200 -none- numeric > class(anomimage) [1] "list"
さて、いよいよこのプロットイメージを空間プロットします。空間プロットにはspパッケージのimage関数を利用し、世界地図はmaptoolsのwrld_simlを利用しましょう。
> image(anomimage) > library(maptools) > data(wrld_simpl) > plot(wrld_simpl, add = T)
これを見ると、確かにプロットはされていますが、同色のカラースキームが利用されているので、値の違いがこれではよく分かりませんし、カラーバーが表示されていないので値の大きさが把握できません。これを改善するために、fieldsパッケージのimg.plotを利用します。
> library(fields) > image.plot(anomimage, nlevel=24) > plot(wrld_simpl, add = T) > grid(col = "grey", lty = 1)
image.plotのnlevelは値を色に割り振る際の分解能を指定します。
上図なら、地球全体でこの夏に何が起きていたかがかなり分かります。この夏は緯度50度線の付近で温度の過度な上昇というアノーマリーが発生していました。ロシアやウクライナ、カザフスタン、モンゴルでは例年の平均気温より5度から6度も高いという異常さを表しています。また日本などの極東やヨーロッパ全域、そしてアメリカ東部も平均気温が2度以上も高くなっています。
結果の詳細はGISSのSurface Temperature Analysisチームの以下の報告にありますので、興味がある方は参照ください。
この中緯度付近で起きた世界的な熱波は、直接的には地球規模の温暖化とは関係がなく、偏西風のブロッキング現象により引き起こされたと説明する研究者もいます。
そして地球規模の温暖化が、そのブロッキング現象をより引き起こしやすくするという地球シミュレーターを利用した数値実験の結果もあるようです。
この夏の異常な温度上昇という現象がどのようなメカニズムで引き起こされたかということは、今後の気候学、気象学の研究を待つことになりますが、このように直近のオープン・データを利用し、少し手間をかけるだけで、気候や気象の状況をモニタリングするということさえ個人でできるようになりました。本当に良い時代になったものだと思います。
次回について
次回もまた、オープンソースのRを利用して、オープンデータにアクセスし、オープンアイディアの活用したいと思います。それでは、また2週間後に会いましょう。
Copyright © ITmedia, Inc. All Rights Reserved.
関連記事
いまさらアルゴリズムを学ぶ意味
コーディングに役立つ！ アルゴリズムの基本（1） コンピュータに「3の倍数と3の付く数字」を判断させるにはどうしたらいいか。発想力を鍛えよう- Zope 3の魅力に迫る
Zope 3とは何ぞや？（1） Pythonで書かれたWebアプリケーションフレームワーク「Zope 3」。ほかのソフトウェアとは一体何が違っているのか？ - 貧弱環境プログラミングのススメ
柴田 淳のコーディング天国 高性能なIT機器に囲まれた環境でコンピュータの動作原理に触れることは可能だろうか。貧弱なPC上にビットマップの直線をどうやって引く？ - Haskellプログラミングの楽しみ方
のんびりHaskell（1） 関数型言語に分類されるHaskell。C言語などの手続き型言語とまったく異なるプログラミングの世界に踏み出してみよう 
ちょっと変わったLisp入門
Gaucheでメタプログラミング（1） Lispの一種であるScheme。いくつかある処理系の中でも気軽にスクリプトを書けるGaucheでLispの世界を体験してみよう