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

# はじめに
デジタル庁が運営するオープンデータの情報ポータルサイト、[データカタログ](https://www.data.go.jp/?lang=ja)データカタログ</a>から2001~2021年の台風の位置データをPythonで作成したRPAでダウンロードし、Rで軽いデータ分析しようと思います。


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/74759b0a-8933-702a-ba77-413818f7fddd.png)

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

![Screenshot 2022-12-20 at 21.07.18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/42849dfa-ca24-a5bd-f463-e1a3bd53c015.png)

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/4ef20a05-3e01-c26c-8861-f0bf8d3d03d4.png)

# Pythonでスクレイピング！(RPA?)
それではPythonでそれぞれのcsvをダウンロードしていきます。今回はPythonライブラリの[Selenium](https://selenium-python.readthedocs.io/)を使って行きます。スクレイピングはSeleniumか<a [Beautiful Soup](https://www.crummy.com/software/BeautifulSoup/bs4/doc/)が王道ですね！Seleniumを使う時はChrome Driverのバージョンを確認しましょう！

```python:Set Up
import time #To give a rest 
from selenium import webdriver 
```
Chrome Dvier を起動
```python: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エレメントを取得して来ます。
```python: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つ開いてみて、ダウンロードへ繋がるリンクを確認します。

```python:Click First URL
driver.get(urls[0])
```
2021年の詳細ページが開かれたので、うまく取得できていそうです。さて、このページからさらにCSVをダウンロードする必要がある様です。。
`[ダウンロード]`ボタンか直書きのリンクでダウンロードできる様なのでこれらを押したいです。↓
![Screenshot 2022-12-20 at 21.51.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/99659a8e-3e88-4164-7238-e99417e6f9ca.png)

`[ダウンロード]`ボタンのHTMLを観察すると、`class` tagの値が`btn btn-primary resource-url-analytics resource-type-None`とユニークっぽくなっているので、これをつかいダウンロードボタンを指定します。
![Screenshot 2022-12-20 at 22.00.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/aef97d93-2329-87ae-579c-37dd8a56295f.png)

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

```python: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/f625f7ef-8554-e645-a68a-de185ac456cd.png)

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

```python: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/d9cfe2f3-9edf-511d-7297-6fb86c82e6fe.png)


# Rでデータ分析

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

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

```r: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')　#ついでに保存
```

### データフレームを確認 
```r
head(df)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/02f75054-8a16-d376-99ca-4379d37694ae.png)



各カラムのメタデータは[台風位置表（CSV形式）フォーマット](https://www.data.jma.go.jp/yoho/typhoon/position_table/format_csv.html#:~:text=50KT%E9%95%B7%E5%BE%84%E6%96%B9%E5%90%91%20%E6%9A%B4%E9%A2%A8%E5%9F%9F,%E4%BB%A5%E4%B8%8B%E3%81%AE%E9%80%9A%E3%82%8A%E8%A1%A8%E8%A8%98%E3%81%99%E3%82%8B%E3%80%82)by気象庁を参照。経度緯度の他に説明関数として使えそうなものがありそうです。
|列番号|内容|
|---|---|
|1|年（UTC）|
|2|月（UTC）|
|3|日（UTC）|
|4|時（UTC）|
|5|台風番号（西暦の下2桁＋その年の発生番号）|
|6|台風名|  
|7|階級   （1: 使用しない）<br> 2: 最大風速が17m/s（34ノット）未満の熱帯低気圧<br> 3: 最大風速が17m/s（34ノット）以上、25m/s（48ノット）未満の熱帯低気圧（台風<sup>※1</sup>）<br> 4: 最大風速が25m/s（48ノット）以上、33m/s（64ノット）未満の熱帯低気圧（台風<sup>※1</sup>）<br> 5: 最大風速が33m/s（64ノット）以上の熱帯低気圧（台風<sup>※1</sup>）<br> 6: 温帯低気圧<br> 7: 中心が北西太平洋（北緯0-60度、東経100-180度）の外にあり、6時間以内に中心が北西太平洋の中に入る熱帯低気圧<sup>※2</sup><br> ※1 台風とは熱帯低気圧のうち北西太平洋内にあり、最大風速が17m/s（34ノット）以上のものである。<br> ※2 北西太平洋の中から外に出た場合、外に出た最初の時刻まで台風位置表に記録するが、その時刻の階級は7ではなく2～6とする。|
|8|緯度（度）|  
|9|経度（度）|  
|10|中心気圧（hPa）|  
|11|最大風速（ノット）<br> ただし、34ノット未満は「0」とする。|  
|12|50KT長径方向<br> 暴風域（最大風速が50ノット以上の領域）は、台風中心から長径方向に（長径半径－短径半径）/2ずらした点を中心点とする半径（長径半径＋短径半径）/2 の円内で表し、暴風域の長径方向を8方位で以下の通り表記する。<br> 1: 北東、2: 東、3: 南東、4: 南、5: 南西、6: 西、7: 北西、8: 北<br> ただし、暴風域が台風の中心に対して同心円状である場合には「9」とし、13列目の長径、14列目の短径に同じ数値を表記する。|  
|13|50KT長径（海里）<br> ただし、暴風域がない（最大風速が50ノット未満の）場合は、「0」とする。|  
|14|50KT短径（海里）<br> ただし、暴風域がない（最大風速が50ノット未満の）場合は、「0」とする。|  
|15|30KT長径方向<br> 強風域（最大風速が30ノット以上の領域）は、台風中心から長径方向に（長径半径－短径半径）/2ずらした点を中心点とする半径（長径半径＋短径半径）/2 の円内で表し、強風域の長径方向を8方位で以下の通り表記する。<br> 1: 北東、2: 東、3: 南東、4: 南、5: 南西、6: 西、7: 北西、8: 北<br> ただし、強風域が台風の中心に対して同心円状である場合には「9」とし、16列目の長径、17列目の短径に同じ数値を表記する。|  
|16|30KT長径（海里）<br> ただし、強風域がない（最大風速が30ノット未満の）場合は、「0」とする。|
|17|30KT短径（海里）<br> ただし、強風域がない（最大風速が30ノット未満の）場合は、「0」とする。|
|18|上陸<br> 1: 上陸（再上陸を含む）または通過の直前、0: それ以外|


### 緯度と観測期間を可視化
観測期間が進むにつれて緯度が上がって行くのがわかります
```r
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/d61f4426-5f6a-13b8-a574-59a2ba9b925e.png)



### 経度と観測期間を可視化
緯度とは違い、観測期間が進むにつれて経度が下がって行くもと、30〜60回目(1日約4回観測)の観測付近グンと上がってあるものがあることがわかります。偏西風の影響かと思われます。
```r
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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/61ac2405-6944-398f-91e6-ba8ac27e58cf.png)




### [Leaflet](https://rstudio.github.io/leaflet/)を使い可視化
緯度軽度情報が取れている様なので、leafletを使い可視化してみます。

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

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/d794407c-7041-72f3-c988-3113458bfeaa.png)

### 2001年だけ
多すぎるので2001年だけ、、
```r: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/3dea161a-041c-7690-a1ac-507356d32fca.png)

### 最大風速で色付け
オレンジだけだと味気ないので最大風速で色付け（緑→低、赤→高）
```r: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/8c9bd351-8c97-2d3e-b174-42d6705b8a72.png)

### 全期間をプロット。最大風速を色と点の大きさで表現
```r: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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/145c9d82-70ba-21bf-3b6b-31868ba10e3c.png)

↓意外と本島にかかってない..?
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/c27a5219-b336-30c7-4cb4-ec0a83bd6387.png)

### 上陸箇所の可視化
やはり西日本が多いですね！
```r: 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](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/3c776566-9d4e-817b-66c8-fac720d8c157.png)


# 台風データを分析


# データを[Delika](https://delika.io/Shu25/thphoon2001-2021)にアップロード

21個のcsvをひとつにまとめたデータDelikaに公開します。公開したデータは[こちら](https://delika.io/Shu25/thphoon2001-2021)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1044781/2328f1f4-d763-0b61-47d6-9ded6f0117eb.png)



以上です！




# ソースコード

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

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

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