LoginSignup
3

More than 1 year has passed since last update.

はじめに

デジタル庁が運営するオープンデータの情報ポータルサイト、データカタログデータカタログから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
  3. You can use dark theme
What you can do with signing up
3