はじめに
デジタル庁が運営するオープンデータの情報ポータルサイト、データカタログデータカタログから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のバージョンを確認しましょう！
import time #To give a rest
from selenium import webdriver
Chrome Dvier を起動
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エレメントを取得して来ます。
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つ開いてみて、ダウンロードへ繋がるリンクを確認します。
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()ファンクションでクリックする処理を行います。
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のリストに向け走らせます。
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が形成されました。
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を使い可視化してみます。
m <- leaflet(year2001) %>%
addTiles() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addCircleMarkers(lng=~経度,lat=~緯度,radius=1,color='orange',weight=2)
2001年だけ
多すぎるので2001年だけ、、
year2001 = df |>
filter(年==2001)
m <- leaflet(year2001) %>%
addTiles() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
addCircleMarkers(lng=~経度,lat=~緯度,radius=1,color='orange',weight=2)
最大風速で色付け
オレンジだけだと味気ないので最大風速で色付け（緑→低、赤→高）
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)
全期間をプロット。最大風速を色と点の大きさで表現
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 = 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)