注目キーワード
         

山口県のコロナ罹患者公表数でヒートマップを作ってみました

山口県のコロナ罹患者公表数でヒートマップを作ってみました
  1. HOME
  2. 雑記 日日是好日
  3. 山口県のコロナ罹患者公表数でヒートマップを作ってみました

オミクロン株は感染力が高いとかいう噂。週末、山に行くにあたり公衆衛生的に安全な場所を知りたかったので、コロナ罹患者公表数からヒートマップを作ってみました。

以下の画像から特設ページに飛べます。山口県の公表データを使っていますのでヒートマップは毎週更新されます。

【山口県】今週のコロナ罹患者公表数ヒートマップ
【山口県】今週のコロナ罹患者公表数ヒートマップ

以下、簡単にヒートマップの作り方もまとめておきます。

目次

ヒートマップの材料

材料は以下の３つです。

制作過程

ほとんどのプログラムは「なにか入力」→「うまいこと処理」→「欲しい形で出力」という流れで作られます。今回は「コロナの患者数」→「うまいこと処理」→「ヒートマップで出力」の流れとなりました。以下、作成過程を見てゆきますが、参考のため事前調査の段階から記しておきます。

入力：山口県のコロナ罹患者公表数

まずはデータ探しから。

フリーで使える新型コロナの患者データはあるのか？

そもそもフリーで使える新型コロナの患者データがなければ何も始まりません。まずは検索からスタートです。

  • “山口県 新型コロナ データ”で検索すると、「陽性者属性（令和4年1月26日時点） – CKAN」というページが見つかりました。山口県のオープンデータです。
  • データライセンスを見るとCC BY 4.0。営利目的を含めて共有も改変も大丈夫というフリーのライセンスです。のっけから良い感じ。
  • 同データセットには検査実施数、入院死亡数、相談件数等もありましたが、今回の目的には陽性者属性がマッチしているようです。
  • データセットの更新履歴（アクティビティストリーム）を見ると、山口県の厚生課さんが週に一度更新をされていました。

厚生課さんお疲れさまです。

陽性者属性データの中身は？

県が公開しているのですからこれ以上に正確なデータはありません。まずは手作業で陽性者属性データをダウンロードしてみました。

  • 陽性者属性は山口県の陽性者第1号から今週までの陽性者が並ぶ1万数千行のデータで、形式はCSVになっていました。
  • データ項目は[No,全国地方公共団体コード,都道府県名,市区町村名,陽性確定日,公表日,患者年代,患者性別,備考]で、緯度や経度は入っていません。

陽性者の位置情報なんてのは個人情報ですからフリーで公開されるわけもありません。そうなると使えそうなのは市区町村名と公表日、データ件数くらいだなーと扱えるデータにアタリを付けておきます。

データ更新日の取得方法は？

さて、週に1度しか更新されないのに、毎回1万数千行のデータをやり取りしていては互いのサーバに負担がかかります。大きなデータは更新された時にだけ取りに行きたいものです。

  • まずはphpのget_headers() 関数で陽性者属性の更新日付取得を試みましたが、戻ってきたデータにはLastModifiedが含まれていませんでした。
  • 困りましたが、陽性者属性のURLを見ているうちにCKAN（オープンソースのデータポータルプラットホーム）が使われていることに気が付きました。
  • 山口県もオープンソースを使うようになったんやねえというわけで、CKANのドキュメントを読んでみます。
  • CKANのAPIにリソースの状態を返してくれるものがありましたので、山口県独自サイトのURLに合わせて呼び出してみました。
  • ＯＫ。これで陽性者属性の更新日とURLをJSON形式で取得できるようになりました。

戻ってきたURLから陽性者属性を取得できることも確認できたので、入力の仕様についてはこれで大体OKです。
次は出力を見てみましょう。

出力：Esri JavaScript

地図上にデータを表示して直感的に分かるような何かを作りたいのですが、諸々検索をしてEsri JavaScriptにたどり着きました。

  • “Esri Javascript Sample”で検索するとSample Code | ArcGIS API for JavaScript 4.22がトップに表示されます。
    • このページにアクセスすると地図表示のサンプルコードをたくさん見ることができます。その他に[API Reference]からもサンプルコードのページを探れます。
  • コロナ陽性者数を地図表示する場合は、以下のようなバルーンが表示されるタイプが一般的かと思いますが、
  • 僕はヒートマップを作ってみたい。[Search]から”sample heatmap”を検索して以下のページを見つけました。
  • Sample Codeページでは[Explore in the sandbox]からJavaScriptのサンプルコードを見ることができます。ヒートマップのサンプルコードを眺めてみました。
    • urlでCSVファイルを取り込んでいます
    • rendererはコメントを見る限り単位ピクセル内のポイント数で色を変える処理っぽい
    • mapは背景地図の設定、viewは初期の中心座標とズームの設定
    • Legendは色味の説明
    • シンプルで分かりやすいです。
  • 試しにサンプルコードのCSVファイルをダウンロードしてみましたが、やはり緯度経度以外にrendererに使われていそうなものは見つかりませんでした。
  • となると、サンプルコードの入力処理を担うCSVLayerをもう少し掘り下げる必要があります。
    • [API Reference]からCSVLayerの説明に飛んでみました。
    • ほほー。CSVファイルのヘッダー行にlatitude, longitudeが書かれていると、その列が緯度、経度として認識されるようです。
    • サンプルコードのCSVには地震発生個所の緯度経度が格納されてるんで、その点の分布によって色が決まるってことね

なるほど、CSVさえ上手に作れば出力はサンプルのJavaScriptコードをほぼそのまま使えそうです。

処理：ざっくりとプログラム構成

入力と出力が分かったのでその間の処理を埋めていきます。メインとなる処理は「陽性者属性からEsri JavaScriptが読めるCSVを作ること」になります。ざっとプログラム構成はこんな感じですね。

【php】
// 入力
陽性者属性データが更新されているかを確認｛
　更新されていれば陽性者属性データを取得
　// メイン処理
　市町の中心の緯度経度を定義
　陽性者属性データの直近１週間分のみを対象として
　市町の中心の緯度経度＋分布のデータを作成
　CSVに吐き出す
｝
【Esri JavaScript】
// 出力
CSVを読み込んでヒートマップを表示する

分布について

分布については、円の中心ほど密になるデータを作るのが簡単なので以下の式を使いました。/1000とCORRECT35を除けばありがちな式です。

$r=rand(0,100) / 1000;
$theta=deg2rad(rand(0,360));
$x=$r * cos($theta);
$y=$r * sin($theta) * CORRECT35;

最終的なx,yは10進緯度経度単位です。半径rは市町の中心からのコロナの影響範囲ということになりますが、最大で0.1度、アバウトに10km程度としておきました。
CORRECT35は地図上で円を作るための補正係数です。北緯35度での緯度経度1度あたりの長さ比から予め計算しておきます。
latitude+y,longitude+xとすれば、市町の中心から半径rの範囲内のデータ点の分布ができます。

市町の中心（円の中心）をどこに置くか

山口県には19市町村があり陽性者属性にも市区町村名が入っています。さて、各市町の中心をどこに置くべきなんでしょう？

最初は各役場の緯度経度を充てていたのですが、どうにも海よりです。もう少し意味のあるデータはないものかと、”市町村 重心”等で検索していたところ「人口重心」という言葉を発見しました。地理的な重心よりもそれっぽいので円の中心に置くのにふさわしそうです。

我が国の人口重心－平成27年国勢調査結果から – 総務省統計局
https://www.stat.go.jp/data/kokusei/topics/topi102.html

上記サイトに「各都道府県及び市区町村の人口重心（エクセル：2,318KB）」があります。平成27年とちょっとデータが古いのですが、令和2年の国勢調査の人口重心がまだ上がってきていないので、まあ、役場の位置データよりましだろうということでこちらを採用することにしました。人口重心を用いることで各市町の中心が少しだけ海から内陸に入りました。

というわけで完成したのが以下のサイトです。

【山口県】今週のコロナ罹患者公表数ヒートマップ
【山口県】今週のコロナ罹患者公表数ヒートマップ

これで、どのあたりの山を歩けば安全なのか「見える化」ができました。

今後の課題

個々の罹患者の緯度経度が公開されない（公開されるわけがない）以上、ヒートマップはランダム点で作るしかありません。つまりヒートマップがどういう形になるかは分布の計算しだいです。ここでもう一度、分布の計算式を見てみることにします。

$r=rand(0,100) / 1000;
$theta=deg2rad(rand(0,360));
$x=$r * cos($theta);
$y=$r * sin($theta) * CORRECT35;

上記式のrは人口重心からのコロナの影響範囲ということになりますが、例えば人口の多い市町または罹患者の多い市町についてはrを長くするなど、rを可変にすることで他市との関係が変わりヒートマップに変化が起きて面白いのではないかと考えています。ただし中間点が濃くなる可能性もあり、これを避けるために恣意的なパラメータが必要になるかもしれません。そうなると泥沼です。山口県人の平均的な通勤距離など、rにはなんらかの意味を持たせた方が簡単かもしれません。

もう一点。今はthetaを0-360度の間でランダムに発生させています。rが海に及ばない山口市・美祢市などはこれでよいのですが、山口県の多くの市町は沿岸部に人口が集中しており0-360では円の半分が海にかかってしまうのです。これを避けるために、地形と中心点を考慮して市町ごとにrandの範囲を変え内陸のみを対象にするというのも、また一つ面白い結果を生むでしょう。

ちなみに本章を「今後の課題」としましたが、自分の用途としては現状で十分だと思っていますので、僕自身はこれ以上の改修を行うつもりはありません。より面白いものを作りたいという方は、ぜひ自作のヒートマップに挑戦してみてください。地図は面白いよ。

謝辞

最後に材料を提供いただいた皆様に謝辞を。

山口県様
CC BY 4.0のオープンデータ展開、ありがとうございます。本記事は一つの例ですが、緯度経度付きのデータを公開いただければいっそう県民の見える化が進むことと存じます。どうぞこれからもよろしくお願いいたします。
また厚生課様、毎週のデータ集計・更新は大変な作業かと存じます。ご担当も健康に気を付けられ滞りなく公開が続けられますよう、ご健勝をお祈りいたします。

Esri様
平易かつ格好いいたくさんのサンプルの提示をありがとうございます。ヒートマップのサンプルを勝手に使って良かったのか分かりませんが、習作として見逃していただけると嬉しいです。もし利用料金が発生するようでしたら速攻で削除しますのでご連絡ください。皆が地図に興味を持ちより高度な処理のできるArcGISが爆売れしますよう、貴社の益々のご発展をお祈り申し上げます。

総務省統計局様
国勢調査から都道府県・市区町村毎の人口重心なんてーものを計算していることを、今回初めて知りました。昔は手作業で計算されていたとか。お疲れ様です。そうした歴史ある人口重心とは思いますが、Excel表の度分秒に10進を併記いただけるとより平易なGISの利用促進につながるのではないかと存じます。（いや、変換のひと手間が面倒なだけなんですけどね） 令和２年度の国勢調査結果を楽しみにしております。

ありがとうございました。

       雑記 日日是好日前後記事