川をデータ分析してみた


## 概要
神通川について定量的に評価したいと思い、公開されているデータを用いてデータ分析してみた。

今回評価した項目は

- 河川延長(長さ)
- 流域面積
- 流量
- 水質

の4つ。

## 分析結果

- 河川延長では、全108水系中30位
- 河川流域面積では、全108水系中22位
- 豊水時の流量/渇水時の流量では、全108水系中4位
- 水質では、全108水系中12位

## 分析に使用した公開情報

- [一級河川の河川延長等調](https://www.mlit.go.jp/statistics/details/content/001477723.pdf)
	- 令和３年4月の河川延長に関わるデータがPDF内に表形式で記載されている
- [一級水系の流況](https://www.mlit.go.jp/river/toukei_chousa/kasen_db/pdf/2021/4-1-7.pdf)
	- 平成24年観測のデータだが、各水系ごとに豊水時・渇水時・最大・最小の流量がPDF内に表形式で記載されている
- [河川環境データベース](http://www.nilim.go.jp/lab/fbg/ksnkankyo/dl_00_index.html)
	- 1990年から2019年までの、河川ごとの生物調査の結果をCSV形式でまとめている
- [一級河川の主要地点水質](https://www.mlit.go.jp/river/toukei_chousa/kasen_db/pdf/2021/15-1-2.pdf)
	- 令和1年観測のBOD（生物学的酸素要求量）がPDF内に表形式で記載されている
- [国土数値情報　河川データ](https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-W05.html)
	- 平成21年の河川の流路のデータをGIS形式でまとめている

## 分析手法

### 河川延長、河川流域面積、豊水時の流量/渇水時の流量
PDFをWordで読み込んでから、表の部分をExcelにコピペして`RANK()`使っただけ

### 水質
まず、[一級河川の主要地点水質](https://www.mlit.go.jp/river/toukei_chousa/kasen_db/pdf/2021/15-1-2.pdf)を用いて分析してみた。しかし、BODが1mg以下(鮎が住めるとされる）の水系が多すぎたため、比較するための指標にしづらいと考え却下した。

次に、[河川環境データベース](http://www.nilim.go.jp/lab/fbg/ksnkankyo/dl_00_index.html)を用いて、指標生物と呼ばれる水底に住む生き物を基準に比較した。


```Python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

clearness_rank1 = ["アミカ","ナミウズムシ", "カワゲラ", "サワガニ","ナガレトビケラ","ヒラタカゲロウ","ブユ","ヘビトンボ","ヤマトビケラ","ヨコエビ"]
clearness_rank2 = ["イシマキガイ", "オオシマトビケラ", "カワニナ", "ゲンジボタル", "コオニヤンマ", "コガタシマトビケラ","ヒラタドロムシ", "ヤマトシジミ"]
clearness_rank3 = ["イソコツブムシ", "タニシ", "ニホンドロソコエビ", "シマイシビル", "ミズカマキリ", "ミズムシ"]
clearness_rank4 = ["アメリカザリガニ", "エラミミズ", "サカマキガイ","ユスリカ","チョウバエ"]

teisei_zenkoku = pd.read_csv('/content/drive/MyDrive/底生_全国.csv')

# 調査年度ごとの件数
teisei_zenkoku.groupby('調査年度').count()

# 生物種名を名寄せ
teisei_zenkoku['正規化種名'] = teisei_zenkoku['種名'].str.extract(f'({"|".join(clearness_rank1+clearness_rank2+clearness_rank3+clearness_rank4)})')
```

湿重量と個体数の関係を見ると、なんとなく、反比例の関係にありそう。

```Python
teisei_zenkoku[ teisei_zenkoku['湿重量(mg)'] != 0.0]['湿重量(mg)'].hist(range=[0,20],bins=40)

```
![Kobito.98rHZS.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/181667/edaf01f4-897d-72ed-cd43-789802531179.png "Kobito.98rHZS.png")

`個体数x湿重量`で種別の割合を見ると、ユスリカがかなり多い。

```Python
teisei_zenkoku[ '個体数x湿重量(mg)' ] = teisei_zenkoku['個体数'] * teisei_zenkoku['湿重量(mg)']
ax = teisei_zenkoku[(teisei_zenkoku['湿重量(mg)'] != 0.0) & ~(teisei_zenkoku['湿重量(mg)'].isnull())].groupby('正規化種名').agg({'個体数x湿重量(mg)':'sum'}).plot.pie(figsize=(8,6),subplots=True)
ax[0].get_legend().set_visible(False)
ax
```
![Kobito.JyHGnu.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/181667/54baa9a3-f403-5e64-9499-42540fb102cd.png "Kobito.JyHGnu.png")


階級4かつどの水質でも住むとされるユスリカがノイズになりそうなので除いて、再度`個体数x湿重量`で種別の割合を見る。



```Python
ax = teisei_zenkoku[ (teisei_zenkoku['正規化種名'] != 'ユスリカ') & (teisei_zenkoku['湿重量(mg)'] != 0.0) & ~(teisei_zenkoku['湿重量(mg)'].isnull()) ].groupby('正規化種名').agg({'個体数x湿重量(mg)':'sum'}).plot.pie(figsize=(8,6),subplots=True)
ax[0].get_legend().set_visible(False)
ax
```

![Kobito.luGJGs.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/181667/de3ebb9d-8f2e-d589-5e79-7e10c72660b2.png "Kobito.luGJGs.png")

さらに水系ごとの調査数（ユスリカを除く）を調べる。
どの水系も最低37回は調査されているため、これを元に分析して良さそう。

```Python
teisei_zenkoku[ (teisei_zenkoku['正規化種名'] != 'ユスリカ') & (teisei_zenkoku['湿重量(mg)'] != 0.0) & ~(teisei_zenkoku['湿重量(mg)'].isnull()) ].groupby('水系').count()['種名'].describe()
```

```
count     108.000000
mean      635.685185
std       549.618128
min        37.000000
25%       233.500000
50%       465.000000
75%       886.750000
max      2963.000000
```

水質階級Iの水底生物の割合(`個体数x湿重量`ベース）を調べた。


```
rivers = pd.DataFrame(index=['水系名'], columns=['水系名','ポイント', '水質階級1割合', '水質階級2割合', '水質階級3割合', '水質階級4割合'])
for river, r_group in teisei_zenkoku[ (teisei_zenkoku['正規化種名'] != 'ユスリカ') & (teisei_zenkoku['湿重量(mg)'] != 0.0) & ~(teisei_zenkoku['湿重量(mg)'].isnull()) ].groupby(['水系']):
  weights_rank = np.zeros(4)
  for name, n_group in r_group.groupby('正規化種名'):
    for weight_sum in n_group['個体数x湿重量(mg)']:
      if( name in clearness_rank1 ):
        weights_rank[0] += weight_sum
      elif( name in clearness_rank2 ):
        weights_rank[1] += weight_sum
      elif( name in clearness_rank3 ):
        weights_rank[2] += weight_sum
      else:
        weights_rank[3] += weight_sum
  weights_rank = weights_rank / weights_rank.sum()
  rivers = rivers.append({'水系名': river, '水質階級1割合': weights_rank[0], '水質階級2割合': weights_rank[1], '水質階級3割合': weights_rank[2], '水質階級4割合': weights_rank[3], ignore_index=True)
  
rivers.sort_values('水質階級1割合', ascending=False)[['水系名', '水質階級1割合']].head(15)
```

|index|水系名|水質階級1割合|
|---|---|---|
|62|渚滑川|0\.9935331109592532|
|108|黒部川|0\.9830947917949149|
|64|湧別川|0\.9791590189561669|
|30|姫川|0\.977703469767748|
|21|大井川|0\.9649509142743108|
|100|雲出川|0\.9583155849006285|
|32|安倍川|0\.9534003564087925|
|34|富士川|0\.9532341410833879|
|35|小丸川|0\.9444211296003793|
|94|釧路川|0\.9368800936376812|
|80|網走川|0\.9359008580244127|
|76|神通川|0\.9355345435771318|
|49|斐伊川|0\.9113224995624291|
|22|大分川|0\.8925777394241623|
|42|常呂川|0\.8888403795629263|


水質階級Iの割合をヒストグラムで見ると、広く分布しており指標として望ましいと考えた。

```Python
rivers['水質階級1割合'].hist()
```

![Kobito.mNt9Qg.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/181667/54da213b-6319-3e12-3dae-ee8fdc0dee4f.png "Kobito.mNt9Qg.png")



## 最後に
Pandasを使って分析するのに意外と手こずったが、中々面白い結果が出たのでヨシ！
欲を言うと[国土交通省 水文水質データベース](http://www1.river.go.jp/)を使って分析したかったが、スクレイピング禁止らしいので自重した。
スクレイピング禁止の割には検索UIがアレなので、数年毎のデータをまとめてCSV形式で公開してほしい。

統計的にNGな処理をしてるかもしれないので、指摘していただいたら嬉しいです。

## 参考にした文献

- [水道水の水温](https://www.waterworks.metro.tokyo.lg.jp/suigen/topic/03/)
- [水質調査項目の説明](http://www.thr.mlit.go.jp/bumon/b00037/k00290/river-hp/kasen/plaza/jiko/suisitu_top/suisitu/Yougo/yougo.htm#%E7%94%9F%E7%89%A9%E5%8C%96%E5%AD%A6%E7%9A%84%E9%85%B8%E7%B4%A0%E8%A6%81%E6%B1%82%E9%87%8F%EF%BC%88BOD%EF%BC%89)
- [全国水生生物調査のページ](https://water-pub.env.go.jp/water-pub/mizu-site/mizu/suisei/about/way/text1a.html)