#世界の気象観測データを探す
データ分析で世界の気象観測データを使いたいことがあり，探してみたというお話です．
世界の気象観測データといえば，SYNOP（地上実況気象通報式）やMETAR（定時飛行場実況気象通報式）があります．これらはオープンデータとして利用でき，中には10年以上前のデータまで公開しているWEBサイトもあります．
(e.g. OGIMET)
ただ利用するにはちょっと不便な点があって，電文形式で公開していることが多いんですよね．例えばこんな感じ．
[SYNOP]
202007080300 AAXX 08034 47662 41/70 71709 10253 20218 30009 40036
53008 78082 8327/==
[METAR]
RJTT 080300Z 21019KT 9999 -RA FEW015 SCT050 BKN100 27/21 Q1003
TEMPO 22024G34KT=
私は航空気象に親しんでいるのでMETARは読めますが，SYNOPはほとんど暗号みたいなもので，人がそらで読むもんじゃないです．それにデータ分析で活用するには，CSV形式などテーブルデータとして整理されたものが欲しいですよね．
実は日本の気象庁WEBサイトはしっかりしていて，日本国内のデータであればCSVダウンロードできます．また海外のデータは過去3年分程度しかありませんが，入手可能です．
気象庁｜過去の気象データ・ダウンロード
気象庁｜世界の天候データツール
今回3年以上前の世界の気象観測データが必要で，NOAA（アメリカ海洋大気庁）のデータベースを中心に調べていたところ，ちょうど良いデータを見つけたのでご紹介します．それがGSODです．
#GSODデータについて
GSOD (Global Surface Summary of the Day) とは，世界の気象観測データについて，時間単位のデータをNOAAが独自に日単位データとして集計したものです．具体的な要素としては，最高気温と最低気温，日平均の気温・露点温度・気圧・風などがあります．
観測地点によっては数十年以上の蓄積があり，日単位のデータ分析に使うのであれば，ても有益なデータです．
ただ注意点しなければならないのが、1日の定義が世界標準時(GMT)で00:00〜23:59となっていることです！そうするとローカルタイムとの日時のズレを気にしないといけません．
例えばアメリカ・カリフォルニア州ならば，現地時間とGMTとの差は，サマータイムも考慮するとGMT-7:00 〜 GMT-8:00 となります．
あるT日0:00〜23:59(GMT)の日単位集計は，ローカルタイムにするとGMT-7:00の場合は T-1日17:00〜T日16:59，GMT-8:00の場合は T-1日16:00〜T日15:59 となります．
ここで最高・最低気温を考えると，一般には最高気温は14時ごろ，最低気温は6時ごろに出ます．そうするとGSODデータの集計方法でも，実際の現地時間におけるT日の最高・最低気温は多くの場合正しく集計されていることが期待できます．
このようにローカルタイムとGMTとの差を意識しながら利用する必要があり，完璧ではありませんが，ある程度有用なデータとして使うことができそうです．
#GSODデータをダウンロードする
では早速データを入手してみましょう．GSODサイト内を探索していくと，CSVファイルでダウンロードできるところが見つかりました．年ごとにディレクトリが分かれています．とりあえず2019年を見てみます．
すると…なにこれ？って感じですね…．なんだかよくわからない数字11桁がファイル名となったCSVファイルがズラ〜っと並んでいます．数字は，少なくとも最初の5桁はWMO(世界気象機関)の国際地点番号だと思うんですが，残り6桁がよくわからない．
とりあえず 47662 を検索してみます．検索すると 47662099999.csv というファイル1つだけが引っ掛かります．これは大手町の気象庁本庁の地点番号で，実際にファイルを開いてみると TOKYO という文字が見られます．
このデータをダウンロードしてExcel等で開いてみると，なんだか気温の値が変ですね．1月なのに気温が45度って…これ華氏ですね(笑)．他にも露点温度はあるけど湿度がなく，例えば不快指数などを計算するには一手間必要になってきます．
ということで，このデータを使うにはまずデータが欲しい地点の国際地点番号を探す必要があり，仮にそれが見つかっても単位の変換とか要素の変換とかが必要で，有益なんだけど使い勝手が悪い，そんなデータです．
これは大変だ…．そう思っていたら，何と偶然！GSODへのアクセスと取り扱いを簡単にするRパッケージを見つけてしまいました！
#R package : GSODR
GSODについて検索してて偶然，たまたま見つけたRパッケージ，GSODRを紹介します．これを使うと簡単にGSODを入手して利用することができます．また気温の単位は自動で華氏から摂氏に変換されており，露点温度から湿度も計算されているのでとても便利です．
使い方は本家サイトにまとめられていますが，やはり気象特有の国際地点番号を検索するところなど，慣れないとややこしいところもあります．以下で簡単にご紹介いたします．
まずはパッケージのインストールと呼び出しからです．
# Install and import
install.packages("GSODR")
library(GSODR)
GSODRパッケージでは地点リストも提供してくれます．ここからデータが欲しい地点の国際地点番号を検索する必要があります．ここではアメリカ・カリフォルニア州を例にとってご紹介します．
# See locations list of US, CA
load(system.file("extdata", "isd_history.rda", package = "GSODR"))
US_locatoin <- subset(isd_history, CTRY == "US")
CA_location <- subset(US_locatoin, STATE == "CA")
print( dim(CA_location) )
print( head(CA_location) )
# 出力
454 11
STNID NAME LAT LON CTRY STATE
1: 690020-93218 JOLON HUNTER LIGGETT MIL RES 36.000 -121.233 US CA
2: 690020-99999 JOLON HUNTER LIGGETT MIL RES 36.000 -121.233 US CA
3: 690070-93217 FRITZSCHE AAF 36.683 -121.767 US CA
4: 690140-93101 EL TORO MCAS 33.667 -117.733 US CA
5: 690150-93121 TWENTY NINE PALMS 34.300 -116.167 US CA
6: 690150-99999 TWENTYNINE PALMS 34.300 -116.167 US CA
BEGIN END COUNTRY_NAME ISO2C ISO3C
1: 19640715 19970401 UNITED STATES US USA
2: 20030702 20030801 UNITED STATES US USA
3: 19600404 19930831 UNITED STATES US USA
4: 19890101 19990630 UNITED STATES US USA
5: 19900102 20200407 UNITED STATES US USA
6: 19891115 19891229 UNITED STATES US USA
カリフォルニアだけでも454地点もあります！
ここからさらに地点を絞ります．NAME が都市名っぽいので，試しにサンフランシスコを検索してみます．
# See locations of "SAN FRANCISCO" within its NAME
print( subset(CA_location, grepl("SAN FRANCISCO", NAME)) )
# 出力
STNID NAME LAT LON CTRY STATE
1: 724940-23234 SAN FRANCISCO INTERNATIONAL AIRPORT 37.620 -122.365 US CA
2: 745097-99999 SAN FRANCISCO PBS 37.750 -122.700 US CA
3: 994016-99999 SAN FRANCISCO 37.807 -122.465 US CA
4: 998011-99999 SAN FRANCISCO BAY RESERVE 38.217 -122.017 US CA
5: 998479-99999 SAN FRANCISCO (PIER 1) 37.798 -122.393 US CA
6: 999999-23234 SAN FRANCISCO INTL AP 37.620 -122.365 US CA
BEGIN END COUNTRY_NAME ISO2C ISO3C
1: 19730101 20200407 UNITED STATES US USA
2: 19760130 19960701 UNITED STATES US USA
3: 20080101 20200406 UNITED STATES US USA
4: 20080101 20200406 UNITED STATES US USA
5: 20130101 20200406 UNITED STATES US USA
6: 19480101 19721231 UNITED STATES US USA
SAN FRANCISCOだけでも6件ヒットしました.
次に見たいのは，観測データが存在する期間を表す BEGIN と END です．欲しいデータの期間を含む場所を選びます．2行目と6行目の地点は古いデータしかないようなので，それ以外から1つ選んでみます．1行目はサンフランシスコ国際空港ですが、STNID の最初の5桁 72494 はWMO国際地点番号だし，アメリカでは空港での観測が地方の気象観測として重要に機能していると聞いたことがあるので，ここを選ぶのが良さそうに思います（個人の主観がだいぶ入ってますが）．
ではサンフランシスコ国際空港のGSODデータを実際に入手してみましょう．
# Use get_GSOD function
# Arguments
# years : the year you want to get data
# station : STNID of the location you want to get data
data <- get_GSOD(years = 2011:2016, station = "724940-23234")
print( head(data) )
# 出力
STNID NAME CTRY STATE LATITUDE
1: 724940-23234 SAN FRANCISCO INTERNATIONAL AIRPORT US CA 37.62
2: 724940-23234 SAN FRANCISCO INTERNATIONAL AIRPORT US CA 37.62
3: 724940-23234 SAN FRANCISCO INTERNATIONAL AIRPORT US CA 37.62
4: 724940-23234 SAN FRANCISCO INTERNATIONAL AIRPORT US CA 37.62
5: 724940-23234 SAN FRANCISCO INTERNATIONAL AIRPORT US CA 37.62
6: 724940-23234 SAN FRANCISCO INTERNATIONAL AIRPORT US CA 37.62
LONGITUDE ELEVATION BEGIN END YEARMODA YEAR MONTH DAY YDAY TEMP
1: -122.365 2.4 19730101 20200407 2011-01-01 2011 1 1 1 9.3
2: -122.365 2.4 19730101 20200407 2011-01-02 2011 1 2 2 10.2
3: -122.365 2.4 19730101 20200407 2011-01-03 2011 1 3 3 9.4
4: -122.365 2.4 19730101 20200407 2011-01-04 2011 1 4 4 8.4
5: -122.365 2.4 19730101 20200407 2011-01-05 2011 1 5 5 8.0
6: -122.365 2.4 19730101 20200407 2011-01-06 2011 1 6 6 7.6
TEMP_ATTRIBUTES DEWP DEWP_ATTRIBUTES SLP SLP_ATTRIBUTES STP
1: 24 5.6 24 1017.3 22 16.5
2: 24 7.6 24 1014.5 19 13.8
3: 24 6.4 24 1017.1 19 16.3
4: 24 5.4 24 1022.9 23 22.5
5: 24 4.6 24 1026.0 24 25.4
6: 24 3.8 24 1024.4 24 23.9
STP_ATTRIBUTES VISIB VISIB_ATTRIBUTES WDSP WDSP_ATTRIBUTES MXSPD GUST MAX
1: 24 15.6 24 3.3 24 7.2 NA 16.0
2: 24 15.6 24 4.2 24 6.7 NA 11.1
3: 24 15.0 24 2.3 24 6.7 NA 11.7
4: 24 15.3 24 1.3 24 3.6 NA 11.7
5: 24 15.9 24 0.9 24 3.1 NA 12.8
6: 24 15.8 24 2.2 24 5.1 NA 11.7
MAX_ATTRIBUTES MIN MIN_ATTRIBUTES PRCP PRCP_ATTRIBUTES SNDP I_FOG
1: * 8.3 * 0.8 G NA 1
2: <NA> 8.3 <NA> 2.8 G NA 1
3: <NA> 7.2 <NA> 7.6 G NA 1
4: <NA> 3.9 <NA> 0.0 G NA 1
5: <NA> 3.3 <NA> 0.0 G NA 1
6: <NA> 2.8 <NA> 0.0 G NA 1
I_RAIN_DRIZZLE I_SNOW_ICE I_HAIL I_THUNDER I_TORNADO_FUNNEL EA ES RH
1: NA NA NA NA NA 0.9 1.2 77.7
2: NA NA NA NA NA 1.0 1.2 83.9
3: NA NA NA NA NA 1.0 1.2 81.5
4: NA NA NA NA NA 0.9 1.1 81.4
5: NA NA NA NA NA 0.8 1.1 79.1
6: NA NA NA NA NA 0.8 1.0 76.9
各カラムの説明については，本家サイトの下の方に説明がありますのでご参照下さい．
カラムのうち STNID 〜 END は地点情報なので，不要でしょう．
また 要素名_ATTRIBUTES も品質管理情報みたいなものなので，使わなくて良いと思います．
気象要素としては，TEMP(平均気温)，DEWP(露点温度)，MAX(最高気温)，MIN(最低気温)，RH(平均湿度) あたりが使える要素かなと思います．先に述べた通り，気温の単位は摂氏，また湿度はパッケージが自動で計算してくれています．
なおデータの取得方法として，緯度経度を指定して最寄りのデータを取得する，という方法もあるようです．本家サイトに例が載っていますので，ご参照下さい．
これで世界の気象観測データGSODを使えるようになりました．