統計と可視化の力



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


  • Global Moderator

    花粉の飛散量

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

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

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

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

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

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

    さて、試しに2018年の関東のデータを見てみます。各観測点ごとに、2~5月の4ヶ月間のデータがありました。
    0_1553939204958_pollen2018-kantou.png

    2018年関東の全ての測定点のデータを重ねるとこうなりました。
    0_1553942285467_pollen2018-kantou-overlay.png

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

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

  • Global Moderator

    花粉の飛散量 測定法

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

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

    追加(上記文献の理解)

    装置の構造と機能:

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

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

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

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

    残った問題:

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

  • Global Moderator

    花粉飛散量の時間変化

    1時間ごと、1日24時間の計数データがあるので、詳細な時間変化を見ることが可能です。日ごとに色をわけて部分表示すると以下のようになりました(例):
    0_1553990918224_pollen2018kantou7_1297_1800.png

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

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


  • Global Moderator

    花粉の飛散量 年の違い

    2017年関東のデータについて同じように見てみました:
    0_1553993307746_pollen2017-kantou.png

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



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


  • Global Moderator

    花粉の飛散量 長期変化(関東)

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

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


  • Global Moderator

    花粉の飛散量 データ処理

    再現方法です。測定データの記載されたエクセル形式のファイルを 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]
    

  • Global Moderator

    ニホンウナギの稚魚の採捕量(日本)

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

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

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

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

    などなどです。これを含むウナギの国内流通に関しては、以下の水産庁「うなぎに関する情報」のページが詳しいです。
    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年からは水産庁調べ(池入れ数量-輸入量))

    プロットすると以下のようになりました(水産省の「うなぎをめぐる現状と対策」にあるのと同等です):
    0_1554516464252_unagi-5.png
    縦線は、元になる統計データが前後で異なることを示します。

    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}]
    

    1985年ぐらいまでに大きく減少し、過去15年ほど低調が続いています。最近も減少傾向にあるように見えます。
    0_1554517564260_unagi-5-14yfit.png
    黒線:過去14年間のデータに最小二乗法で直線フィット。

    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と呼ぶそうです。


  • Global Moderator

    ウナギの流通(日本)

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

    資料でまとめてあるのが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 (1億)匹が池に入れられたことが分かります。

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

    @2017年

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

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

    「一般的に蒲焼に用いられるニホンウナギの体重は 200~250g/尾④」「養殖ウナギの適性出荷」重量は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


  • Global Moderator

    シラスウナギの流通(日本)

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

    このスライドにあるシラスウナギの流れ(重量ベース)を Sankey diagram にしてみたのが以下です:
    0_1554594716417_shirasu-flow.png

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

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


    ダイヤグラムの作成には、以下のページを用いました。
    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"}
     	]
    }
    

Log in to reply
 

Looks like your connection to パズルハウス was lost, please wait while we try to reconnect.