search
LoginSignup
2

posted at

updated at

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

はじめに

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

image.png

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

Screenshot 2022-12-20 at 21.07.18.png

HTMLを検証

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

image.png

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をダウンロードする必要がある様です。。
[ダウンロード]ボタンか直書きのリンクでダウンロードできる様なのでこれらを押したいです。↓
Screenshot 2022-12-20 at 21.51.21.png

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

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()

無事にダウンロードされた事を確認しました(左下)
Screenshot 2022-12-20 at 22.19.28.png

これを先ほど取得したすべての年の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を使い分析していきます。
image.png

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)

image.png

各カラムのメタデータは台風位置表(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)

Screenshot 2023-01-20 at 18.41.37.png

経度と観測期間を可視化

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

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

Screenshot 2023-01-20 at 18.43.13.png

Leafletを使い可視化

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

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

image.png

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)

image.png

最大風速で色付け

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

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)

image.png

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

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)

image.png

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

上陸箇所の可視化

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

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)

image.png

台風データを分析

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

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

以上です!

ソースコード

スクレイピング/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)

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
2