世界の気象観測データGSODを使ってみる

#世界の気象観測データを探す

データ分析で世界の気象観測データを使いたいことがあり，探してみたというお話です．

世界の気象観測データといえば，SYNOP（地上実況気象通報式）やMETAR（定時飛行場実況気象通報式）があります．これらはオープンデータとして利用でき，中には10年以上前のデータまで公開しているWEBサイトもあります．
(e.g. [OGIMET](https://www.ogimet.com/home.phtml.en))

ただ利用するにはちょっと不便な点があって，電文形式で公開していることが多いんですよね．例えばこんな感じ．

```
[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年分程度しかありませんが，入手可能です．

[気象庁｜過去の気象データ・ダウンロード](https://www.data.jma.go.jp/gmd/risk/obsdl/index.php)
[気象庁｜世界の天候データツール](http://www.data.jma.go.jp/gmd/cpd/monitor/dailyview/index.php)

今回3年以上前の世界の気象観測データが必要で，NOAA（アメリカ海洋大気庁）のデータベースを中心に調べていたところ，ちょうど良いデータを見つけたのでご紹介します．それがGSODです．


#GSODデータについて
[GSOD (Global Surface Summary of the Day)](https://catalog.data.gov/dataset/global-surface-summary-of-the-day-gsod) とは，世界の気象観測データについて，時間単位のデータを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ファイルでダウンロードできるところが見つかりました](https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/)．年ごとにディレクトリが分かれています．とりあえず2019年を見てみます．

すると…なにこれ？って感じですね…．なんだかよくわからない数字11桁がファイル名となったCSVファイルがズラ〜っと並んでいます．数字は，少なくとも最初の5桁はWMO(世界気象機関)の国際地点番号だと思うんですが，残り6桁がよくわからない．

とりあえず 47662 を検索してみます．検索すると 47662099999.csv というファイル1つだけが引っ掛かります．これは大手町の気象庁本庁の地点番号で，実際にファイルを開いてみると TOKYO という文字が見られます．
このデータをダウンロードしてExcel等で開いてみると，なんだか気温の値が変ですね．1月なのに気温が45度って…これ華氏ですね(笑)．他にも露点温度はあるけど湿度がなく，例えば不快指数などを計算するには一手間必要になってきます．

ということで，このデータを使うにはまずデータが欲しい地点の国際地点番号を探す必要があり，仮にそれが見つかっても単位の変換とか要素の変換とかが必要で，有益なんだけど使い勝手が悪い，そんなデータです．

これは大変だ…．そう思っていたら，何と偶然！GSODへのアクセスと取り扱いを簡単にするRパッケージを見つけてしまいました！


#R package : GSODR
GSODについて検索してて偶然，たまたま見つけたRパッケージ，[GSODR](https://cran.r-project.org/web/packages/GSODR/vignettes/GSODR.html)を紹介します．これを使うと簡単にGSODを入手して利用することができます．また気温の単位は自動で華氏から摂氏に変換されており，露点温度から湿度も計算されているのでとても便利です．

使い方は[本家サイト](https://cran.r-project.org/web/packages/GSODR/vignettes/GSODR.html)にまとめられていますが，やはり気象特有の国際地点番号を検索するところなど，慣れないとややこしいところもあります．以下で簡単にご紹介いたします．

まずはパッケージのインストールと呼び出しからです．

```
# 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
```

各カラムの説明については，[本家サイト](https://cran.r-project.org/web/packages/GSODR/vignettes/GSODR.html)の下の方に説明がありますのでご参照下さい．

カラムのうち STNID 〜 END は地点情報なので，不要でしょう．
また 要素名_ATTRIBUTES も品質管理情報みたいなものなので，使わなくて良いと思います．

気象要素としては，TEMP(平均気温)，DEWP(露点温度)，MAX(最高気温)，MIN(最低気温)，RH(平均湿度) あたりが使える要素かなと思います．先に述べた通り，気温の単位は摂氏，また湿度はパッケージが自動で計算してくれています．

なおデータの取得方法として，緯度経度を指定して最寄りのデータを取得する，という方法もあるようです．[本家サイト](https://cran.r-project.org/web/packages/GSODR/vignettes/GSODR.html)に例が載っていますので，ご参照下さい．

これで世界の気象観測データGSODを使えるようになりました．