はじめに
TechTrainのアドベントカレンダー(https://qiita.com/advent-calendar/2019/techtrain) 11日目の担当で、大学１年生のPython学び始めて2ヶ月半の私が今の実力でオープンデータを使って軽く分析してみたものをまとめます。
内容
各都道府県の花粉症有病率のデータを見て、そんなにスギが生えているイメージがない東京を含む関東圏の広範囲で花粉症有病率が高いことがわかりました。そこでスギの植生と花粉飛散量の間に何か要因があるのだろうと、、、それは風ではないかと予想して、それらを可視化することで予想が正しいか、またどの程度影響しているかを見てみようという趣旨です。
やることが多いので今回はデータの前処理の部分のみまとめます。
使うデータ一覧
- 環境省花粉観測システム(はなこさん)
(http://kafun.taiki.go.jp/DownLoad1.aspx)
期間を2019/2/15〜2019/4/30に指定して北海道と沖縄を除く全国の観測所で観測された花粉飛散数とアメダスのデータ - おまるくんの技術ノート 市区町村コード表(緯度経度付き)
(http://note.omarukun.com/notes/4/)
最終的に日本地図にプロットしたいので観測地点の座標が欲しいため - 風向コード16方位対応表 (気象庁の資料を基に自作)
アメダスのデータの風向が方角ではなく番号で記されていたのでわかりやすいようにするため - スギ花粉症有病率の地域差について
(https://www.jstage.jst.go.jp/article/arerugi/59/1/59_KJ00006063684/_pdf)
都道府県別の花粉症有病率のデータ - 林野庁都道府県別森林率・人工林率
(http://www.rinya.maff.go.jp/j/keikaku/genkyou/h29/1.html)
都道府県別のスギの植生面積のデータ
データ読み込みとデータクレンジング
花粉飛散の地域と量と天気に関するデータ
import pandas as pd
import numpy as np
#データフレームの表示数を設定する
#8行分を表示する
pd.set_option('display.max_rows', 10)
各地方ごとに出力されたCSVファイル読み込んでいきます
# 東北
touhoku = pd.read_csv('./pollen_touhoku.csv', encoding='shift_jis', header=None)
touhoku
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 50210100 | 31312 | 20190215 | 1 | 青森市中央卸売市場 | 1 | 2 | 青森県 | 2201 | 青森市 | 0 | 10 | 3 | -3.2 | 0 | 0 |
| 50210100 | 31312 | 20190215 | 2 | 青森市中央卸売市場 | 1 | 2 | 青森県 | 2201 | 青森市 | 323 | 9 | 5 | -3.4 | 0 | 0 |
| 50210100 | 31312 | 20190215 | 3 | 青森市中央卸売市場 | 1 | 2 | 青森県 | 2201 | 青森市 | 69 | 9 | 4 | -3.1 | 0 | 0 |
| 50210100 | NaN | 20190215 | 4 | 青森市中央卸売市場 | 1 | 2 | 青森県 | 2201 | 青森市 | 86 | NaN | NaN | NaN | NaN | 0 |
| ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ |
| 50620100 | 35141 | 20190308 | 1 | 庄内総合支庁分庁舎 | 2 | 6 | 山形県 | 6426 | 東田川郡三川町 | 0 | 15 | 3 | 3.3 | 0 | 1 |
| 50620100 | 35141 | 20190308 | 2 | 庄内総合支庁分庁舎 | 2 | 6 | 山形県 | 6426 | 東田川郡三川町 | 45 | 16 | 3 | 1.1 | 1 | 0 |
| 50620100 | 35141 | 20190308 | 3 | 庄内総合支庁分庁舎 | 2 | 6 | 山形県 | 6426 | 東田川郡三川町 | 8 | 14 | 3 | 2.9 | 0 | 0 |
| 50620100 | 35141 | 20190308 | 4 | 庄内総合支庁分庁舎 | 2 | 6 | 山形県 | 6426 | 東田川郡三川町 | 4 | 13 | 3 | 3.1 | 0 | 0 |
23400 rows × 16 columns
これを残りの関東(kanto)、中部(chubu)、関西(kansai)、中国・四国(chugoku_shikoku)、九州(kyusyu)と読み込んでいきます。
同じこと且つ長いのでコードと結果は割愛します。
余分な列を削除する
東北地方のデータフレームを見てわかるように列のラベルが不明なので一つずつ調べて要らない列を削除します
touhoku.drop([0, 1, 5, 6, 13, 14, 15], axis=1, inplace=True)
kanto.drop([0, 1, 5, 6, 13, 14, 15], axis=1, inplace=True)
chubu.drop([0, 1, 5, 6, 13, 14, 15], axis=1, inplace=True)
kansai.drop([0, 1, 5, 6, 13, 14, 15], axis=1, inplace=True)
chugoku_shikoku.drop([0, 1, 5, 6, 13, 14, 15], axis=1, inplace=True)
kyushu.drop([0, 1, 5, 6, 13, 14, 15], axis=1, inplace=True)
欠損値の処理
まず欠損値の有無の確認
print(touhoku.isnull().sum())
print(kanto.isnull().sum())
print(chubu.isnull().sum())
print(kansai.isnull().sum())
print(chugoku_shikoku.isnull().sum())
print(kyushu.isnull().sum())
欠損値のある行の削除
touhoku.dropna(how='any', inplace=True)
kanto.dropna(how='any', inplace=True)
chubu.dropna(how='any', inplace=True)
kansai.dropna(how='any', inplace=True)
chugoku_shikoku.dropna(how='any', inplace=True)
kyushu.dropna(how='any', inplace=True)
削除後のデータを表示する
touhoku
| 2 | 3 | 4 | 7 | 8 | 9 | 10 | 11 | 12 |
|---|---|---|---|---|---|---|---|---|
| 20190215 | 1 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 0 | 10 | 3 |
| 20190215 | 2 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 323 | 9 | 5 |
| 20190215 | 3 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 69 | 9 | 4 |
| 20190215 | 5 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 16 | 9 | 5 |
| ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ |
| 20190430 | 21 | 株式会社江東微生物研究所 東北中央研究所 | 福島県 | 7204 | いわき市 | 0 | 8 | 2 |
| 20190430 | 22 | 株式会社江東微生物研究所 東北中央研究所 | 福島県 | 7204 | いわき市 | 0 | 9 | 3 |
| 20190430 | 23 | 株式会社江東微生物研究所 東北中央研究所 | 福島県 | 7204 | いわき市 | 0 | 8 | 3 |
| 20190430 | 24 | 株式会社江東微生物研究所 東北中央研究所 | 福島県 | 7204 | いわき市 | 0 | 9 | 2 |
23310 rows × 9 columns
欠損値のある行を削除することで東北地方に関しては90行減りました
これを先ほどと同様に各地方見ていきます
各データにどの市区町村が含まれるのかを確認する
それぞれのデータにたくさんの観測地点がありますが、一体いくつあるのかを観測地点がある市区町村の数でカウントしてみます
# 北海道と沖縄を除く各地方における観測地点の市区町村の数を確認
print(len(touhoku.groupby(9).count()))
print(len(kanto.groupby(9).count()))
print(len(chubu.groupby(9).count()))
print(len(kansai.groupby(9).count()))
print(len(chugoku_shikoku.groupby(9).count()))
print(len(kyushu.groupby(9).count()))
# 北海道と沖縄を除く全国の観測地点の市区町村の数の合計を確認
print(len(touhoku.groupby(9)) + len(kanto.groupby(9)) + len(chubu.groupby(9)) + len(kansai.groupby(9)) + len(chugoku_shikoku.groupby(9)) + len(kyushu.groupby(9)))
東北： 13
関東： 20
中部： 23
関西： 19
中国・四国： 21
九州： 20
合計： 116
各地方ごとに分かれたデータフレームを１つに結合
結合したら約21万行のデータフレームになりました
alljapan = pd.concat([touhoku, kanto, chubu, kansai, chugoku_shikoku, kyushu])
alljapan
| 2 | 3 | 4 | 7 | 8 | 9 | 10 | 11 | 12 |
|---|---|---|---|---|---|---|---|---|
| 20190215 | 1 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 0 | 10 | 3 |
| 20190215 | 2 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 323 | 9 | 5 |
| 20190215 | 3 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 69 | 9 | 4 |
| 20190215 | 5 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 16 | 9 | 5 |
| ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ |
| 20190430 | 21 | 鹿児島県大隅地域振興局本庁舎 | 鹿児島県 | 46203 | 鹿屋市 | 8 | 11 | 1 |
| 20190430 | 22 | 鹿児島県大隅地域振興局本庁舎 | 鹿児島県 | 46203 | 鹿屋市 | 4 | 10 | 1 |
| 20190430 | 23 | 鹿児島県大隅地域振興局本庁舎 | 鹿児島県 | 46203 | 鹿屋市 | 0 | 11 | 1 |
| 20190430 | 24 | 鹿児島県大隅地域振興局本庁舎 | 鹿児島県 | 46203 | 鹿屋市 | 0 | 10 | 1 |
208273 rows × 9 columns
列のラベルがこのままではわかりにくいので列名を編集します
alljapan.columns = ['年月日', '時間', '観測地点', '都道府県', '市区町村コード', '市区町村', '花粉飛散数(個/m^3)', '風向コード', '風速(m/s)']
風向のデータをわかりやすくしたい
風向コードと16方位で表される風向の対応表を用意する
wind_direction = pd.read_csv('./wind_direction_code.csv', encoding='UTF-8', header=0)
wind_direction
| 風向 | 風向コード |
|---|---|
| 北北東 | 1 |
| 北東 | 2 |
| 東北東 | 3 |
| 東 | 4 |
| 東南東 | 5 |
| 南東 | 6 |
| 南南東 | 7 |
| 南 | 8 |
| 南南西 | 9 |
| 南西 | 10 |
| 西南西 | 11 |
| 西 | 12 |
| 西北西 | 13 |
| 北西 | 14 |
| 北北西 | 15 |
| 北 | 16 |
この風向コードの対応表と観測データを結合します
結合にあたりデータ型を揃えるためにalljapanの風向列のデータ型をfloatからintに変換
alljapan['風向コード'] = alljapan['風向コード'].astype(np.int64)
風向コードを軸にして2つのデータフレームを結合する
alljapan2 = pd.merge(
alljapan, wind_direction,
on = '風向コード',
how = 'left'
)
alljapan2
| 年月日 | 時間 | 観測地点 | 都道府県 | 市区町村コード | 市区町村 | 花粉飛散数(個/m^3) | 風向コード | 風速(m/s) | 風向 |
|---|---|---|---|---|---|---|---|---|---|
| 20190215 | 1 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 0 | 10 | 3 | 南西 |
| 20190215 | 2 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 323 | 9 | 5 | 南南西 |
| 20190215 | 3 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 69 | 9 | 4 | 南南西 |
| 20190215 | 5 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 16 | 9 | 5 | 南南西 |
| ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ |
| 20190430 | 21 | 鹿児島県大隅地域振興局本庁舎 | 鹿児島県 | 46203 | 鹿屋市 | 8 | 11 | 1 | 西南西 |
| 20190430 | 22 | 鹿児島県大隅地域振興局本庁舎 | 鹿児島県 | 46203 | 鹿屋市 | 4 | 10 | 1 | 南西 |
| 20190430 | 23 | 鹿児島県大隅地域振興局本庁舎 | 鹿児島県 | 46203 | 鹿屋市 | 0 | 11 | 1 | 西南西 |
| 20190430 | 24 | 鹿児島県大隅地域振興局本庁舎 | 鹿児島県 | 46203 | 鹿屋市 | 0 | 10 | 1 | 南西 |
風向コード0番目が含まれていたので削除する
風向コード0は「静穏」を示すが今回は不要なため
alljapan3 = alljapan2[alljapan2['風向コード'] != 0]
座標に関するデータを結合したい
市区町村コードと座標の対応表を読み込む
データフレームをマークダウン用に出力したので緯度経度の小数点以下が2桁しか表示されていませんが実際にはもっと細かい値が出ています。
city_lon_lat = pd.read_csv('./citycode_coordinate_correspondence_table.csv', encoding='UTF-8', header=0)
city_lon_lat
| 市区町村コード | 都道府県 | 市区町村 | Unnamed: 3 | 緯度 | 経度 |
|---|---|---|---|---|---|
| 11011 | 北海道 | 札幌市 | 中央区 | 43.06 | 141.3 |
| 11029 | 北海道 | 札幌市 | 北区 | 43.09 | 141.3 |
| 11037 | 北海道 | 札幌市 | 東区 | 43.08 | 141.4 |
| ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ |
| 473758 | 沖縄県 | 宮古郡 | 多良間村 | 24.67 | 124.7 |
| 473812 | 沖縄県 | 八重山郡 | 竹富町 | 24.29 | 123.9 |
| 473821 | 沖縄県 | 八重山郡 | 与那国町 | 24.47 | 123.0 |
1896 rows × 6 columns
余分な列を削除する
city_lon_lat.drop(['都道府県', '市区町村', 'Unnamed: 3'], axis=1, inplace=True)
市区町村コードの下一桁はチェックデジットなので左の5桁が一般的に使われている市区町村コードになるので下一桁を取り除く
# city_lon_latの市区町村コード列のデータ型をintからstrに変換
city_lon_lat['市区町村コード'] = city_lon_lat['市区町村コード'].astype(str)
# city_lon_latの市区町村コード列の下一桁を削除して再代入する
city_lon_lat['市区町村コード'] = city_lon_lat['市区町村コード'].str[:-1]
# city_lon_latの市区町村コード列のデータ型をstrからintに変換
city_lon_lat['市区町村コード'] = city_lon_lat['市区町村コード'].astype(int)
alljapan3とcity_lon_latを「市区町村コード」列を軸にしてマージする
join_alljapan = pd.merge(
alljapan3, city_lon_lat,
on='市区町村コード',
how='left'
)
join_alljapan
| 年月日 | 時間 | 観測地点 | 都道府県 | 市区町村コード | 市区町村 | 花粉飛散数(個/m^3) | 風向コード | 風速(m/s) | 風向 | 緯度 | 経度 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 20190215 | 1 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 0 | 10 | 3 | 南西 | 40.82 | 140.7 |
| 20190215 | 2 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 323 | 9 | 5 | 南南西 | 40.82 | 140.7 |
| 20190215 | 3 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 69 | 9 | 4 | 南南西 | 40.82 | 140.7 |
| ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ | ・・・ |
| 20190430 | 22 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 0 | 13 | 1 | 西北西 | 40.82 | 140.7 |
| 20190430 | 23 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 4 | 11 | 1 | 西南西 | 40.82 | 140.7 |
| 20190430 | 24 | 青森市中央卸売市場 | 青森県 | 2201 | 青森市 | 0 | 14 | 1 | 北西 | 40.82 | 140.7 |
データ型と各列のデータ数を確認
join_alljapan.info()
年月日 201422 non-null int64
時間 201422 non-null int64
観測地点 201422 non-null object
都道府県 201422 non-null object
市区町村コード 201422 non-null int64
市区町村 201422 non-null object
花粉飛散数(個/m^3) 201422 non-null float64
風向コード 201422 non-null int64
風速(m/s) 201422 non-null float64
風向 201422 non-null object
緯度 196486 non-null float64
経度 196486 non-null float64
緯度・経度に欠損値があるようであるので具体的に確認すると3つの観測地点の市区町村の市区町村コードに誤りがあったようです。
join_alljapan[join_alljapan['緯度'].isnull()]
# 静岡県庁の所在地、静岡市葵区の市区町村コードが22201ではなく正しくは22101のようであるので修正する
alljapan3.市区町村コード[alljapan3.市区町村コード == 22201] = 22101
# 熊本市医師会ヘルスケアセンターの所在地、熊本市中央区の市区町村コードが43201ではなく正しくは43101のようであるので修正する
alljapan3.市区町村コード[alljapan3.市区町村コード == 43201] = 43101
# 鹿児島県姶良・伊佐地域振興局伊佐庁舎の所在地、鹿児島県伊佐市の市区町村コードが46209ではなく正しくは46224のようであるので修正する
alljapan3.市区町村コード[alljapan3.市区町村コード == 46209] = 46224
改めてalljapan3とcity_lon_latを市区町村コード列を軸にしてマージする
join_alljapan = pd.merge(
alljapan3, city_lon_lat,
on='市区町村コード',
how='left'
)
データ型と各列のデータ数を確認
今度はデータ数が全て一致し、欠損値なくマージできたようです
join_alljapan.info()
年月日 201422 non-null int64
時間 201422 non-null int64
観測地点 201422 non-null object
都道府県 201422 non-null object
市区町村コード 201422 non-null int64
市区町村 201422 non-null object
花粉飛散数(個/m^3) 201422 non-null float64
風向コード 201422 non-null int64
風速(m/s) 201422 non-null float64
風向 201422 non-null object
緯度 201422 non-null float64
経度 201422 non-null float64
観測地点の数の検算
print(len(join_alljapan.groupby('市区町村コード')))
116
スギの植生のデータ
sugi_area = pd.read_csv('./sugi_h29.csv', encoding='UTF-8', header=4)
sugi_area
52 rows × 24 columns
sugi_area.drop(['Unnamed: 0', 'Unnamed: 1', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '20＋'], axis=1, inplace=True)
sugi_area.columns = ['都道府県', 'スギ植生面積(m^2)']
# 欠損値のある行の削除
sugi_area.dropna(how='any', inplace=True)
sugi_area = sugi_area.drop(47)
sugi_area
| 都道府県 | スギ植生面積(m^2) |
|---|---|
| 北海道 | 32,128 |
| 青森県 | 198,690 |
| 岩手県 | 200,525 |
| ・・・ | ・・・ |
| 宮崎県 | 224,828 |
| 鹿児島県 | 153,298 |
| 沖縄県 | 237 |
都道府県別花粉症有病率のデータ
pollen_prevalence = pd.read_csv('./pollen_prevalence.csv', encoding='UTF-8', header=0)
pollen_prevalence = pollen_prevalence.drop(47)
pollen_prevalence
| 都道府県 | 有病率 |
|---|---|
| 北海道 | 2.2 |
| 青森 | 12.5 |
| 岩手 | 12.1 |
| ・・・ | ・・・ |
| 宮崎 | 8.2 |
| 鹿児島 | 12.1 |
| 沖縄 | 6.0 |
CSVファイルで取り込んだデータのデータクレンジングは以上です。
集計
次はクレンジング済みのデータを集計して必要なデータを求めていきます。
各市町村における花粉飛散数の平均
pollen_amount = pd.DataFrame(join_alljapan.groupby(['都道府県']).mean()['花粉飛散数(個/m^3)'])
# pollen_amount.insert(0, '市区町村', join_alljapan['市区町村'])
pollen_amount
各観測地点の風向の集計
各観測地点で毎時間観測されている風向を観測した期間において集計することで、その観測地点における大まかな風向を求めたい、、、ということで集計していきます。
# 都道府県・市区町村・風向でグルーピングして各方位ごとの観測回数をまとめてデータフレームを新しく作成する
count_wind = pd.DataFrame(join_alljapan.groupby(['都道府県', '市区町村', '風向']).size())
count_wind.columns = ['風向の各方位における観測回数']
# 「風向の全方位における合計観測回数」を新しい列として追加する
count_wind['風向の全方位における合計観測回数'] = join_alljapan.groupby(['都道府県', '市区町村']).size()
# 「観測された風向の各方角の全体に占める割合(%)」を新しい列として追加する
count_wind['観測された風向の各方角の全体に占める割合(%)'] = count_wind['風向の各方位における観測回数'] / count_wind['風向の全方位における合計観測回数'] * 100
count_wind
こんな具合です
花粉飛散数と風向の集計のテーブルを作成しましたが、列の左の方にデータフレームとして列になっていない状態です。
そこに都道府県の列を追加できたら可視化の準備が完了します。
ただ、そこに手間取っているので今回はここまでにします。。。
最後に
今回はやりたかった分析の一連の流れのうちデータ前処理のみまとめました。
あと少し改良してデータを可視化できる状態にして、実際に可視化するのが自分の宿題です。。。
まとめられたらQiita更新します。
関連
オープンデータから花粉飛散量とスギの植生と風の影響を分析する(可視化編)
https://qiita.com/TakuTaku36/items/e9fa09b38f86be188070