【RとPythonの合わせ技】過去20年間の台風データを取得して可視化してみる

はじめに

デジタル庁が運営するオープンデータの情報ポータルサイト、データカタログデータカタログから2001~2021年の台風の位置データをPythonで作成したRPAでダウンロードし、Rで軽いデータ分析しようと思います。

「台風位置表のCSVデータ」というのが欲しいデータです。2001~2021年なでのデータがある様なので、これらのcsvをRPAでダウンロードします。21ファイルあるのでダウンロードが大変＋毎年新しいデータが出た際に再度ダウンロードするのが面倒なので、RPAを作り自動ですべてのファイルをダウンロードします。

HTMLを検証

ディベロッパーツールを機動し、tagを観察。。。
titletagの"台風位置表のCSVデータ（20XX年）" というフィルターで持って来れそうです。そこから作成したリストからそれぞれの年のデータのダウンロードページに飛ぶURLを取得しFor Loop でダウンロードしていきます。

Pythonでスクレイピング！(RPA?)

それではPythonでそれぞれのcsvをダウンロードしていきます。今回はPythonライブラリのSeleniumを使って行きます。スクレイピングはSeleniumか<a Beautiful Soupが王道ですね！Seleniumを使う時はChrome Driverのバージョンを確認しましょう！

Set Up

import time #To give a rest 
from selenium import webdriver 

Chrome Dvier を起動

Activate Chrome Driver

main_directory = "/PATH/TO/YOUR/DIRECTORY" 
driver = webdriver.Chrome('/PATH/TO/YOUR/chromedriver 3') #Activate Chrome Driver
driver.get("https://www.data.go.jp/data/dataset/mlit_20140919_0747") # Access to the website

データ一覧ページから、'台風位置表のCSVデータ(20XX年)'のHTMLエレメントを搾取しました。
2行目でそのHTMLエレメントの中からhrefアトリビュートを搾取します。hrefアトリビュートはURLを格納しているので、この手順でダウンロードページに飛ぶURLを集めます。
集めたHTMLエレメントすべてからのURLを集めるためにFor Loopでリスト(urls)を作成しています。
ちなみに'台風位置表のCSVデータ'という部分を任意の文字列にすることで、その文字列を含んだHTMLエレメントを取得して来ます。

Get URLs

elems = driver.find_elements_by_xpath("//a[contains(@title,'台風位置表のCSVデータ')]") #Get elements title tags that contain '台風位置表のCSVデータ'
urls = [elem.get_attribute("href") for elem in elems] #Get all URLs from the list 

うまく'台風位置表のCSVデータ'をclassに含むHTMLエレメントが取得できたか確認したいので試しに1つ開いてみて、ダウンロードへ繋がるリンクを確認します。

Click First URL

driver.get(urls[0])

2021年の詳細ページが開かれたので、うまく取得できていそうです。さて、このページからさらにCSVをダウンロードする必要がある様です。。
[ダウンロード]ボタンか直書きのリンクでダウンロードできる様なのでこれらを押したいです。↓

[ダウンロード]ボタンのHTMLを観察すると、class tagの値がbtn btn-primary resource-url-analytics resource-type-Noneとユニークっぽくなっているので、これをつかいダウンロードボタンを指定します。

Pythonに戻りclass tagの値が"btn btn-primary resource-url-analytics resource-type-None"のエレメントをサーチし、click()ファンクションでクリックする処理を行います。

Try Clicking The Download URL

driver.get(urls[0])
button = driver.find_elements_by_xpath("//a[contains(@class,'btn btn-primary resource-url-analytics resource-type-None')]")
button[0].click()

無事にダウンロードされた事を確認しました（左下）

これを先ほど取得したすべての年のURLのリストに向け走らせます。

Run Download URLs

for url in urls:
    driver.get(url)
    button = driver.find_elements_by_xpath("//a[contains(@class,'btn btn-primary resource-url-analytics resource-type-None')]")
    button[0].click()
    time.sleep(3)

じゃんじゃんダウンロードされて...
綺麗に2001~2021年のデータが取れました！次はこのデータをRを使い分析していきます。

Rでデータ分析

ただデータをとってもつまらないので少し分析をしてみます。このデータに多いで一番に解決したいであろう課題は「台風の進路」でしょう。観測された日を1日目とし、経度緯度の推移を可視化してみましょう。

まずはダウンロードしたcsvを読み取り、結合します。19,014行x18列のdfが形成されました。

Read&Marge CSVs

df = data.frame()
file_list <- list.files("/PATH/TO/YOUR/typhonon_data")
for (csv in file_list){
  temp = read.csv(paste('typhonon_data',csv,sep = '/'),header = T,fileEncoding = "cp932")
  df = rbind(df,temp)
}
write.csv(df,'all_typhonon_data.csv')　#ついでに保存

データフレームを確認

head(df)

各カラムのメタデータは台風位置表（CSV形式）フォーマットby気象庁を参照。経度緯度の他に説明関数として使えそうなものがありそうです。

列番号	内容
1	年（UTC）
2	月（UTC）
3	日（UTC）
4	時（UTC）
5	台風番号（西暦の下2桁＋その年の発生番号）
6	台風名
7	階級 （1: 使用しない） 2: 最大風速が17m/s（34ノット）未満の熱帯低気圧 3: 最大風速が17m/s（34ノット）以上、25m/s（48ノット）未満の熱帯低気圧（台風※1） 4: 最大風速が25m/s（48ノット）以上、33m/s（64ノット）未満の熱帯低気圧（台風※1） 5: 最大風速が33m/s（64ノット）以上の熱帯低気圧（台風※1） 6: 温帯低気圧 7: 中心が北西太平洋（北緯0-60度、東経100-180度）の外にあり、6時間以内に中心が北西太平洋の中に入る熱帯低気圧※2 ※1 台風とは熱帯低気圧のうち北西太平洋内にあり、最大風速が17m/s（34ノット）以上のものである。 ※2 北西太平洋の中から外に出た場合、外に出た最初の時刻まで台風位置表に記録するが、その時刻の階級は7ではなく2～6とする。
8	緯度（度）
9	経度（度）
10	中心気圧（hPa）
11	最大風速（ノット） ただし、34ノット未満は「0」とする。
12	50KT長径方向 暴風域（最大風速が50ノット以上の領域）は、台風中心から長径方向に（長径半径－短径半径）/2ずらした点を中心点とする半径（長径半径＋短径半径）/2 の円内で表し、暴風域の長径方向を8方位で以下の通り表記する。 1: 北東、2: 東、3: 南東、4: 南、5: 南西、6: 西、7: 北西、8: 北 ただし、暴風域が台風の中心に対して同心円状である場合には「9」とし、13列目の長径、14列目の短径に同じ数値を表記する。
13	50KT長径（海里） ただし、暴風域がない（最大風速が50ノット未満の）場合は、「0」とする。
14	50KT短径（海里） ただし、暴風域がない（最大風速が50ノット未満の）場合は、「0」とする。
15	30KT長径方向 強風域（最大風速が30ノット以上の領域）は、台風中心から長径方向に（長径半径－短径半径）/2ずらした点を中心点とする半径（長径半径＋短径半径）/2 の円内で表し、強風域の長径方向を8方位で以下の通り表記する。 1: 北東、2: 東、3: 南東、4: 南、5: 南西、6: 西、7: 北西、8: 北 ただし、強風域が台風の中心に対して同心円状である場合には「9」とし、16列目の長径、17列目の短径に同じ数値を表記する。
16	30KT長径（海里） ただし、強風域がない（最大風速が30ノット未満の）場合は、「0」とする。
17	30KT短径（海里） ただし、強風域がない（最大風速が30ノット未満の）場合は、「0」とする。
18	上陸 1: 上陸（再上陸を含む）または通過の直前、0: それ以外

緯度と観測期間を可視化

観測期間が進むにつれて緯度が上がって行くのがわかります

ggplot(df, aes(time, lat,color = factor(typhoon_number))) + 
  geom_line(aes(group = factor(typhoon_number)),show.legend = FALSE)

経度と観測期間を可視化

緯度とは違い、観測期間が進むにつれて経度が下がって行くもと、30〜60回目(1日約4回観測)の観測付近グンと上がってあるものがあることがわかります。偏西風の影響かと思われます。

ggplot(df, aes(time, long,color = factor(typhoon_number))) + 
  geom_line(aes(group = factor(typhoon_number)),show.legend = FALSE)

Leafletを使い可視化

緯度軽度情報が取れている様なので、leafletを使い可視化してみます。

Plot with Leaflet

m <- leaflet(year2001) %>% 
  addTiles() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(lng=~経度,lat=~緯度,radius=1,color='orange',weight=2)

2001年だけ

多すぎるので2001年だけ、、

Plot with Leaflet for Year 2001

year2001 =  df |> 
  filter(年==2001) 

m <- leaflet(year2001) %>% 
  addTiles() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(lng=~経度,lat=~緯度,radius=1,color='orange',weight=2)

最大風速で色付け

オレンジだけだと味気ないので最大風速で色付け（緑→低、赤→高）

Plot with Leaflet for Year 2001 with Color Based on 最大風速

year2001 =  df |> 
  filter(年==2001) 

pal <- colorNumeric(
  palette = colorRampPalette(c('green','red'))(length(year2001$最大風速)), 
  domain = year2001$最大風速)

m <- leaflet(year2001) %>% 
  addTiles() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(lng=~経度,lat=~緯度,radius=1,color=pal(year2001$最大風速),weight=2)

全期間をプロット。最大風速を色と点の大きさで表現

Plot with Leaflet for Year 2001-2022 with Color Based on 最大風速

pal <- colorNumeric(
  palette = colorRampPalette(c('green','red'))(length(df$最大風速)), 
  domain = df$最大風速)

m <- leaflet(df) %>% 
  addTiles() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(lng=~経度,lat=~緯度,radius=df$最大風速/10,color=pal(df$最大風速),opacity = 0.0001)

↓意外と本島にかかってない..?

上陸箇所の可視化

やはり西日本が多いですね！

onLand == 1

onLand = df |> 
  filter(上陸 == 1)
pal <- colorNumeric(
  palette = colorRampPalette(c('green','red'))(length(onLand)), 
  domain = onLand$最大風速)

m <- leaflet(onLand) %>% 
  addTiles() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(lng=~経度,lat=~緯度,radius=onLand$最大風速/10,color=pal(onLand$最大風速),opacity = 0.001)

台風データを分析

データをDelikaにアップロード

21個のcsvをひとつにまとめたデータDelikaに公開します。公開したデータはこちら

以上です！

ソースコード

スクレイピング/RPA(Python)

ソースコード

# Set Up
import time #To give a rest 
from selenium import webdriver 

#Activate Chrome Driver
main_directory = "/PATH/TO/YOUR/DIRECTORY" 
driver = webdriver.Chrome('/PATH/TO/YOUR/chromedriver 3') #Activate Chrome Driver
driver.get("https://www.data.go.jp/data/dataset/mlit_20140919_0747") # Access to the website

#Get URLs
elems = driver.find_elements_by_xpath("//a[contains(@title,'台風位置表のCSVデータ')]") #Get elements title tags that contain '台風位置表のCSVデータ'
urls = [elem.get_attribute("href") for elem in elems] #Get all URLs from the list 

#Try Clicking The Download URL
driver.get(urls[0])
button = driver.find_elements_by_xpath("//a[contains(@class,'btn btn-primary resource-url-analytics resource-type-None')]")
button[0].click()

#Run Download URLs 
for url in urls:
    driver.get(url)
    button = driver.find_elements_by_xpath("//a[contains(@class,'btn btn-primary resource-url-analytics resource-type-None')]")
    button[0].click()
    time.sleep(3)

driver.close() # Close Chrome Driver

可視化(R)

ソースコード

library(dplyr)
library(leaflet)
library(ggplot2)

#Read CSVs
df = data.frame()
file_list <- list.files("/PATH/TO/YOUR/typhonon_data")
for (csv in file_list){
  temp = read.csv(paste('typhonon_data',csv,sep = '/'),header = T,fileEncoding = "cp932")
  df = rbind(df,temp)
}
write.csv(df,'all_typhonon_data.csv')　#ついでに保存

#Plot with Leaflet
m <- leaflet(year2001) %>% 
  addTiles() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(lng=~経度,lat=~緯度,radius=1,color='orange',weight=2)

# Plot with Leaflet for Year 2001
year2001 =  df |> 
  filter(年==2001) 

m <- leaflet(year2001) %>% 
  addTiles() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(lng=~経度,lat=~緯度,radius=1,color='orange',weight=2)

#Plot with Leaflet for Year 2001 with Color Based on 最大風速
year2001 =  df |> 
  filter(年==2001) 

pal <- colorNumeric(
  palette = colorRampPalette(c('green','red'))(length(year2001$最大風速)), 
  domain = year2001$最大風速)

m <- leaflet(year2001) %>% 
  addTiles() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(lng=~経度,lat=~緯度,radius=1,color=pal(year2001$最大風速),weight=2)

#Plot with Leaflet for Year 2001-2022 with Color Based on 最大風速
pal <- colorNumeric(
  palette = colorRampPalette(c('green','red'))(length(df$最大風速)), 
  domain = df$最大風速)

m <- leaflet(df) %>% 
  addTiles() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(lng=~経度,lat=~緯度,radius=df$最大風速/10,color=pal(df$最大風速),opacity = 0.0001)