統計と可視化の力

ソム

３日間平均降水量

雨が降り続く効果を加味し、３日間平均降水量をプロット見てみましょう。時系列データについて、移動平均を行います。

dateMovingAverage[data_, r_Integer] := Module[{dates, values},
   {dates, values} = Transpose[data];
   Transpose[{dates[[r ;;]], MovingAverage[values, r]}]
   ];

DateListPlot[dateMovingAverage[precipitationLogKobe, 3], 
 PlotStyle -> Directive[Opacity[1]], AspectRatio -> 1/10, 
 PlotRange -> All, ImageSize -> 1000,
 FrameLabel -> {"西暦", "３日間平均 \n降水量 [mm]"}, 
 PlotLabel -> "気象庁データ　神戸気象台", BaseStyle -> {FontSize -> 16}]
Export["s:/precipitationLogKobe-m3.png", %]

DateListPlot[dateMovingAverage[precipitationLogKyoto, 3], 
 PlotStyle -> Directive[Opacity[1]], AspectRatio -> 1/10, 
 PlotRange -> All, ImageSize -> 1000,
 FrameLabel -> {"西暦", "３日間平均 \n降水量 [mm]"}, 
 PlotLabel -> "気象庁データ　京都気象台", BaseStyle -> {FontSize -> 16}]
Export["s:/precipitationLogKyoto-m3.png", %]

このプロットで見ても、過去に何度か経験したことのある大雨であることが分かります。例えば神戸では、歴史的に何度か水害がありました。http://www.city.kobe.lg.jp/life/town/river/suigaisonae/02kako_02.html

手塚治虫の「アドルフに告ぐ」でも描写されているように、昭和１３年（1938年）には阪神大水害がありました。このとは上のページによると、（7月3日：49.6mm、4日：141.8mm、5日：270.4mm）の降水量があったそうです。他にも、昭和３６年（1961年）（6月24日：76.8mm、25日：195.2mm、26日：127.7mm、27日：72.4mm）、昭和４２年（1967年）（7月9日：319.4mm）と、台風と梅雨前線の関係する水害が生じています。これらは上の一連のプロットでも見ることが出来ます。

ソム

降水量の偏り

これまでのようなプロットの欠点は、飛び値（他と較べて特に大きな値）を強調して見てしまうところかもしれません。実際には多くの日で降水量は小さいのですが、それらの値はプロット上では点として重なってしまい、見た目の効果が相対的に小さくなります。このようなときは、例えばヒストグラムを取ってやると、値の分布をより正確に理解できるでしょう。

Histogram[precipitationLogKobe[[All, 2]], "Log", Frame -> True, 
 FrameLabel -> {"一日あたり降水量 [mm]　対数軸", "日数"}, ImageSize -> 300, 
 BaseStyle -> {FontSize -> 16}, PlotLabel -> "気象庁データ　神戸気象台"]
Export["s:/histogram-precipitationLogKobe.png", %]

Histogram[precipitationLogKyoto[[All, 2]], "Log", Frame -> True, 
 FrameLabel -> {"一日あたり降水量 [mm]　対数軸", "日数"}, ImageSize -> 300, 
 BaseStyle -> {FontSize -> 16}, PlotLabel -> "気象庁データ　京都気象台"]
Export["s:/histogram-precipitationLogKyoto.png", %]

一日あたり100 mmを越える降水量はかなり稀であることが分かります。

riffraff

@ソムさん
ビジュアル化してみると、アナログなクリシェで見えないものがクリアに見えてきますね。
ちなみに当方の最大24時間雨量は166㎜、庭が池になっています。

ソム

花粉の飛散量

今年も花粉が飛ぶ季節になり、人によってはつらい時期になっています。さて、花粉の飛散はどのように変化しているのでしょうか？

２００３年ごろから、環境省が自動測定器を使って各地の花粉飛散量を測定しています。

環境省花粉観測システム（愛称：はなこさん）
http://kafun.taiki.go.jp

「花粉自動測定器の概要」によると、吸引ポンプで吸い込んだ大気にレーザーを照射し、散乱光の数から粒子の数を得るそうです。１時間ごとに、平均個数/m^3 が測定値として得られます。従来の、プレパラートを一日放置して、染色後に顕微鏡で花粉を計数するダーラム法に較べると、人手が必要無くほぼリアルタイムに測定できるという利点があります。一方、見て計数しているわけではないので、花粉の種類（スギ、ヒノキ、など）が分からない、花粉以外の粒子と区別つかない場合もある欠点があります。そのため、観測地点の近傍のアメダスやレーダーによる降雨・降雪などの測定環境の情報と合わせてデータを公開してます。公開されている花粉飛散量データには、降雪や黄砂などによる、データの欠測があります。おそらく、花粉と粒径が近い粒子がそれらの条件で多くなり、安定した測定ができなくなるのでしょう（※）。

※ただし文献によっては、３０μｍ前後の粒子以外はカットするので、黄砂（～４μm）と区別できるとしているものもあります。http://www.erc.pref.fukui.jp/center/publish/report/2007/3-3-6.pdf

ライブラリのページには、２００３年～２０１８年の測定データが公開されています。
http://kafun.taiki.go.jp/Library.html
２００３年時点では関東の限られた場所での測定しかありませんが、２０１８年では全国に広がり、観測点の増えています。

さて、試しに２０１８年の関東のデータを見てみます。各観測点ごとに、２～５月の４ヶ月間のデータがありました。

２０１８年関東の全ての測定点のデータを重ねるとこうなりました。

これらを見て分かったのは以下のことでした：

花粉の飛散量測定データは、測定時間によって大きく変動する。
花粉の飛散量測定データは、測定点によって、かなり違う。時期的な傾向も、飛散量の値も。

ソム

花粉の飛散量　測定法

先に言及した文献に、測定装置の開発に関する論文が引用されていました：

レーザー光学手法を用いた新しい花粉計測法とその成果、佐橋ら 2003
https://www.jstage.jst.go.jp/article/jriet1972/32/3/32_3_191/_pdf/-char/ja

追加（上記文献の理解）

装置の構造と機能：

大気を上部から吸引し、途中で砂を除くトラップを介して光学系へ。
大気は、内部を循環する清浄大気（フィルターを通した大気）をシースとする、シースフローとなる（花粉を含む測定したい大気の外側を鞘のように清浄大気が包み光学系を流れる）。
光学系では、シート状に成形したレーザー光により、粒子の通過時間、前方散乱光強度、側方散乱光強度を測定する。

細胞の解析で使われるフローサイトメーターの空気版というふうです。

従来のパーティクル・カウンターを花粉測定に用いるときの問題と解決：

花粉と砂や繊維状ダストの区別が付かない。
→光学系から得られたデータから、ほぼ球形に近い粒子を計数する。
花粉計測に適した粒径レンジがない。
→スギ・ヒノキを計数するために、28～35μmの粒子だけを計数する。
雨を直接吸い込むと機械の性能を維持できない。
→百葉箱や専用ラックに組み込む。

残った問題：

花粉以外にも大きさが同じくらいの空中浮遊粒子も測定してしまう。
実際の現場での予期せぬトラブル：雪や黄砂の影響で誤計数したことがある。

ソム

花粉飛散量の時間変化

１時間ごと、１日２４時間の計数データがあるので、詳細な時間変化を見ることが可能です。日ごとに色をわけて部分表示すると以下のようになりました（例）：

一日の中でも変化があるのが分かります。そのパターンは、山型に変化している場合もあれば、ピークが複数あるように見える日もあります。飛散量（測定値）の多い日もあれば、そうでない日もあるのが分かります。

補足
ここでの飛散量とは、スギやヒノキなどの花粉の飛散の時間変化に、気象状況（風）などの影響が加わり、さらに花粉以外の粒子のノイズも加わり、最終的に装置付近の大気中の粒子として計数された値と考えられます。すなわち、ここで見ているのは、生物学的な花粉の飛散の変化ではありません。一日ごとに変化のパターンがあるように見える理由には、植物の概日周期や、気象現象の日周期などが考えられると思います。

ソム

花粉の飛散量　年の違い

２０１７年関東のデータについて同じように見てみました：

まず、気づいたのは、相変わらず、測定地点による違いがかなりあります。それから、全体的に２０１８年の値に較べると低く見えるということです。このような長期的な変化を定量するにはどうすればよいでしょうか？

riffraff

@ソムさん
いろいろVoodooな技は思いつきますが決め手に欠ける上筋悪。
中でも筋悪な一例、スペクトル分析に持ち込み、補助パラメータ（風向、植生etc）てんこ盛りで、ビッグデータ多地点モデルを作る。だれかやってそうで怖いですが。コンピュータ神カルトです。!(^^)!

ソム

花粉の飛散量　長期変化（関東）

各測定地点で、代表値として各年の中央値を用い、長期変化が見えないか試みました。欠測は先に除いています。

同じ測定地点のデータを線で結んでいます。地点の増減や変更があります。グローバルな増減が見える年もありました。中央値で見る限り、平均気温やＣＯ２濃度のような、長期的な増減の傾向があるわけではなさそうです。

ソム

花粉の飛散量　データ処理

再現方法です。測定データの記載されたエクセル形式のファイルを http://kafun.taiki.go.jp/Library.html の下のリンクからダウンロードします。ただし、扱いにくい点が有るので注意です:

2003-2005年の書式と 2006年以降の書式が異なる。
2003-2014年はxls形式だが、2015年以降はxlsx形式である。
2015年以降と以前で欠測を示す数値コードに異なる点がある。

一括した扱いができないのが不便です。

それぞれの xls または、xlsx ファイルをMathematica上で以下の関数を使ってインポートします：

readType1Data[year_Integer, dataPath_String](*2003-2005*):= 
  Module[{data = Import[dataPath][[1]],
    points, dates,
    pollenDensities},
   points = DeleteCases[Drop[data[[1]], 1], ""];
   dates = 
    DateObject[{ToString[year] <> "年" <> #, {"Year", "年", "Month", 
         "月", "Day", "日", "Hour", "時"}}, "Hour", TimeZone -> 9] & /@ 
     data[[2 ;;]][[All, 1]];
   pollenDensities = Table[
     {points[[i]],
      Transpose[{dates, fix /@ data[[2 ;;]][[All, 1 + i]]}]},
     {i, Length@points}
     ]
   ];

readType2Data[dataPath_String](*2006-*):= 
  Module[{data = Import[dataPath][[1]],
    points, dates,
    pollenDensities},
   points = DeleteCases[Drop[data[[2]], 4], ""];
   dates = 
    DateObject[Round@#, "Hour", TimeZone -> 9] & /@ 
     data[[3 ;;]][[All, 1 ;; 4]];
   pollenDensities = Table[
     {points[[i]],
      Transpose[{dates, fix /@ data[[3 ;;]][[All, 4 + i]]}]},
     {i, Length@points}
     ]
   ];

(*欠測への対応*)
fix[x_] := Which[x == "", Missing[],
   NumericQ[x], If[Round@x < 0,
    Switch[Round@x,
     -9998, Missing["Snow"],
     -9997, Missing["Kousa"], -9991, Missing["Kousa"],
     -9996, Missing["Anomaly"], -9992, Missing["Anomaly"]]
    , Round@x]
   ];

readType1Dataは、2003-2005年のデータ、それ以降は readType2Data を使います。結果、欠測の処理も行った時系列データが得られます。このとき、日本語文字を含むファイル名のエクセルファイルではImport関数がエラーになる（Mathematica version 11の既知のバグ）ので、アルファベットや数字に置き換えておきます。pollen2018-kantou.xlsx などです。そして以下のように読み込み、プロットを行いました。

(*この前に SetDirectory[...] として、データのファイルのあるフォルダに移動しておく*)

pollen2018kantou = readType2Data["pollen2018-kantou.xlsx"];

Labeled[Grid[
  Partition[
   Table[DateListPlot[pollen2018kantou[[i, 2]], PlotRange -> All, 
     PlotLabel -> pollen2018kantou[[i, 1]],
     FrameLabel -> {Automatic, 
       "個/\!\(\*SuperscriptBox[\(m\), \(3\)]\)"}, ImageSize -> 250, 
     BaseStyle -> {FontSize -> 11}], {i, Length@pollen2018kantou}], 
   UpTo[4]]],
 "花粉飛散量 2018 関東\n環境省 花粉自動測定器による観測データ", Top]

ソム

ニホンウナギの稚魚の採捕量（日本）

ニホンウナギは稚魚（シラスウナギ）が河口付近で採捕され、室内外のプール（池）に入れ育てることで養殖されています。国内に流通するウナギの40%程度を占めるようです（養殖生産量：20,979t / 供給量：53,324t @2017年）。近年、シラスウナギの採捕量が減少しているとの報道がなされています。

ややこしいのは、ウナギには、いくつかの種があり、それぞれが世界に分布することです：

ニホンウナギ Anguilla japonica
マリアナ海溝で産卵、東アジアに分布
ＩＵＣＮ　絶滅危惧種　https://www.iucnredlist.org/species/166184/1117791
環境省レッドリスト　絶滅危惧　https://www.env.go.jp/press/16264.html
ヨーロッパウナギ Anguilla anguilla
サルガッソ海で産卵、ヨーロッパに分布
ＩＵＣＮ　絶滅危惧種　https://www.iucnredlist.org/species/60344/45833138
ワシントン条約で貿易取引が規制
アメリカウナギ

すなわち、いろいろなパターンがあります：

国内で稚魚が採捕され、国内で養殖されるニホンウナギ
国外で採捕され、国内に稚魚として輸入され、養殖されるニホンウナギ
国外で採捕され、養殖され、国内に輸入されるニホンウナギ
国外で採捕され、養殖され、国内に輸入されるヨーロッパウナギ

などなどです。これを含むウナギの国内流通に関しては、以下の水産庁「うなぎに関する情報」のページが詳しいです。
http://www.jfa.maff.go.jp/j/saibai/unagi.html

うなぎをめぐる現状と対策(2019年3月）
http://www.jfa.maff.go.jp/j/saibai/attach/pdf/unagi-109.pdf

これらの量的な流れを可視化するのにふさわしいのは、以下にあるようなフロー図（※）だと思います（ウナギとは関係ありません）：
https://flowcharts.llnl.gov/commodities/energy

しかし、ここではウナギの現状の指標として最も単純な、国内のシラスウナギの採捕量を見てみましょう。データは水産庁のページに公開されています：

http://www.jfa.maff.go.jp/j/saibai/attach/xls/unagi-5.xlsx
（2002年までは漁業・養殖業生産統計年報による・2003年からは水産庁調べ（池入れ数量－輸入量））

プロットすると以下のようになりました（水産省の「うなぎをめぐる現状と対策」にあるのと同等です）：

縦線は、元になる統計データが前後で異なることを示します。

unagi = Import["unagi-5.xlsx"];
DateListPlot[Transpose[{
   DateObject[{#}] & /@ Range[1957, 2018],
   unagi[[1, 2 ;;]][[All, 3]]}], Joined -> True,
 FrameLabel -> {"西暦", "シラスウナギ採捕量　　（トン）    "}, ImageSize -> 400,
 BaseStyle -> {FontSize -> 16}, 
 GridLines -> {{DateObject[{2002}]}, None}]

１９８５年ぐらいまでに大きく減少し、過去１５年ほど低調が続いています。最近も減少傾向にあるように見えます。

黒線：過去１４年間のデータに最小二乗法で直線フィット。

last14y = unagi[[1, 2 ;;]][[-14 ;; -1, 3]];

fit = FindFit[last14y, a t + b, {a, b}, t]
(* {a -> -0.542198, b -> 18.3165}*)

Show[
 ListPlot[last14y, PlotStyle -> PointSize@Medium, Frame -> True, 
  FrameLabel -> {"years", "シラスウナギ採捕量　　（トン）    "}, 
  PlotRangePadding -> Scaled[.05], BaseStyle -> {FontSize -> 12}], 
 Plot[a t + b /. fit, {t, 0, 14}, PlotStyle -> Black], 
 ImageSize -> 250]

※Sankey diagramと呼ぶそうです。

ソム

ウナギの流通（日本）

水産庁の資料から、ウナギの量的な流れを考えてみましょう。

資料でまとめてあるのが2017年までのデータなので、2017年を例に見ると、以下のようになります。

@2017年

ニホンウナギ池入れ数量のうち、国内漁獲量 15.5トン①
ニホンウナギ池入れ数量のうち、輸入 4.1 トン①
その他の種のウナギ池入れ数量 0.1トン（シラスウナギ換算）②

① http://www.jfa.maff.go.jp/j/saibai/attach/pdf/unagi-35.pdf
② http://www.jfa.maff.go.jp/j/saibai/attach/pdf/unagi-45.pdf

量的な流れを見るときに、トン数は保存量ではありません。ウナギは育つと重量が増すからです。保存されるのは匹数です。②の資料では、匹数と重量をシラスウナギで換算しています。このときの換算係数を計算すると、0.18g/シラスウナギ一匹としていることが分かります。この係数を用いて見積もると、2017年の国内の養殖ウナギは、(15.5+4.1+0.1)t/0.18 g = 1.0 x 10^8 （１億）匹が池に入れられたことが分かります。

一方、これらの養殖を含む最終的なウナギの供給量の統計は以下のようです。

@2017年

養殖生産量 20979トン③
漁業生産量 71トン③
輸入量　32274トン③

③ http://www.jfa.maff.go.jp/j/saibai/attach/pdf/unagi-109.pdf

「一般的に蒲焼に用いられるニホンウナギの体重は 200～250ｇ/尾④」「養殖ウナギの適性出荷」重量は200g⑤らしいので、ここでは、200g換算とすると、最終的な生産量は、53324 t/200 g = 2.4 x 10^8 匹となります。この換算係数の場合には、養殖生産が 9.5 x 10^7 匹になるので、池入れの95%が生産量になることになります。

④ http://www.pref.shizuoka.jp/sangyou/sa-130/documents/631.pdf
⑤ https://www.fra.affrc.go.jp/unagi/unagi_shigen.pdf

ソム

シラスウナギの流通（日本）

水産庁の「うなぎをめぐる現状と対策」は毎年更新されているようです。この資料には、前年のデータを元にした「シラスウナギの流通の実情」というスライドがあります。

このスライドにあるシラスウナギの流れ（重量ベース）を Sankey diagram にしてみたのが以下です：

「報告漏れ等」は、スライドによると以下のように説明されています：

都府県等からの聞き取りにより、次の原因が指摘されている。
① 採捕者が他人に自分の採捕数量を知られたくない（優良な採捕場所を秘密にしたい、大漁へのねたみを回避したい等）、報告するのが面倒などの理由で報告しない
② 採捕者が指定された出荷先以外へ、より高い価格で販売し、その分を報告しない
③ 無許可による採捕（いわゆる密漁）

ダイヤグラムの作成には、以下のページを用いました。
https://ricklupton.github.io/d3-sankey-diagram/

このページの "Try It!" という欄に、以下を張り付けると、ダイヤグラムが自動生成されます。

{
	"nodes" : [
		{"id":"shirasu-reported","title":"採捕\n報告"},
		{"id":"shirasu-unreported","title":"報告\n漏れ等"},
		{"id":"shirasu-captured-by-oneself","title":"自己採捕"},
		{"id":"shirasu-imported","title":"輸入"},
		{"id":"intermediary","title":"集荷業者"},
		{"id":"company","title":"会社・個人事業主"},
		{"id":"cooperative","title":"漁協"},
		{"id":"unprofit","title":"非営利法人"},
		{"id":"group","title":"任意団体"},
		{"id":"cultured","title":"池入れ"}
	],
 	"links" : [
 	{"source":"shirasu-reported","target":"intermediary","value":5.3,"type":"1"},
      	{"source":"shirasu-unreported","target":"intermediary","value":3.6,"type":"1"},
      	{"source":"shirasu-imported","target":"intermediary","value":5.2,"type":"1"},
      	{"source":"intermediary","target":"company","value":11.9,"type":"2"},
      	{"source":"intermediary","target":"cooperative","value":1.9,"type":"2"},
      	{"source":"intermediary","target":"unprofit","value":0.1,"type":"2"},
      	{"source":"intermediary","target":"group","value":0.2,"type":"2"},
      	{"source":"company","target":"cultured","value":11.9,"type":"3"},
      	{"source":"cooperative","target":"cultured","value":1.9,"type":"3"},
      	{"source":"unprofit","target":"cultured","value":0.1,"type":"3"},
      	{"source":"group","target":"cultured","value":0.2,"type":"3"},
        {"source":"shirasu-captured-by-oneself","target":"cultured","value":0.03,"type":"3"}
 	]
}