QGISで世帯数推計を行ってみた

世帯推計が必要になった背景

　令和5年7月中旬に豪雨災害があり、所属している課（福祉担当）が被災世帯支援業務を担当することになり「り災証明書が届き次第当日処理で給付手続きをとる」で対応するよう指示があった。
　通常業務に被災者支援の業務が上積みになるため、定数では対応できないことが想定され、どの程度の業務が増えるのか推定する必要があると考え、推定の基礎となる被災世帯の推定を行うこととした。
　その時、どのような作業を行ったか思い出しながら手順を解説する。
　気象変動による台風、豪雨のほか、大地震も頻発しているなか、被害を受けた範囲の確定は必要なスキルになりつつあると感じる。
　また、この手法はGISの基本的な操作だけで行え、応用の広い手法であり、通常業務でも応用が効くもの考える。

業務環境

　これほどの大規模災害は数十年ぶりであり、全く体制は整っていない。
　それどころか通常業務も人員不足の状況のため、体制整備のため業務量を見積もる必要があり世帯推計を行った。
　つまり市役所としてオフィシャルな推定では無く、担当部門の業務量推定のために行った。本来であれば初期段階で、担当部門が推定値を公表し、それに基づいて適切な体制整備が行われることが本来の姿である。
　ちなみに最初の推定で対象世帯数が2万件以上になりそうなため必要な人員や体制が必要と福祉保健部幹部に相談したが、２か月以上全く対策はなされなかった。この災害対応で私も含め２人が体調を壊し休職し、さらに他の1人が退職している。

引用・出典について

文中で使用している地図は国土地理院の地理院タイルを使用している。
　

世帯推計の方法

推計の考え方

　対象が多すぎるため全数調査を行うとしても、結果が出るまで待つことはできない。
　サンプル調査を行うとしても、その人員も無く時間も無かった。
　そこでフェルミ推定だが、目安を得られればよしとし、考えられる最大値を求めることとした。

採用した世帯推計の考え方

　あくまで根拠のある概算をするには、どうすれば良いかと考え次のようにした。

当初は街区・字単位で集計した

　災害発生時から、被害世帯から相談の電話が入ることが想定されていたため、相談記録を残すこととしていた。
　この相談記録は、MS Accessでデータベースを作成したが、その際、住所（街区、字単位）を庁内の標準コードで記録できるようにしていた。
　庁内では住所単位で人口推計を行っているので、二つの情報を組み合わせれば世帯数は出せる。
　次の地図で赤い太線が街区のサンプル。例えば秋田市山王一丁目や二丁目だったり楢山佐竹町など。

　そこで、当初段階の世帯数推計は、次のような考え方とした。

1. 被災の相談を住所コードに紐付ける。
2. 住所コード毎に相談件数をカウントする。
   3.相談件数が１以上（相談件数 > 0）の住所コードを抽出する。
3. 抽出された住所コードの世帯数を集計する。
   　この方法で災害発生後すぐに推計した被災世帯数は約2万9000世帯だった。
   　この方法の問題点は住所単位の集計のため、地理上の範囲は考慮されていない。
   　例えば範囲の広い街区や字では被災範囲が街区や字の一部しか無いとしても、全体が被災地域として集計することになり、過大に見積もることになる。

GISを用いて推定することに変更

　街区・字単位では、実態に即していない恐れがあり、より地理的な広がりに即した推定方法が望ましい。
　総務省統計局では国勢調査の集計結果を国土をメッシュ状に分割したデータを公開している( https://www.stat.go.jp/data/mesh/index.html )。
　最小のメッシュは250m四方となる。
　地図上に250mメッシュをかさねあわせると次のようになる。赤い線がメッシュを表す。

　集計単位を街区・字から250mメッシュに変更しすれば、より地理的な範囲で集計できる。
　この集計作業は手作業でもできるが、地図上にデータを位置づけ（マッピング）、被災したメッシュを集計するのは大変な作業になる。
　なお、災害から一ヶ月ほどたつと、り災証明書発行のための調査が始まり被災世帯データ（自称CSVファイル¹）が防災部門から届くようになった。

GISでの集計方法

　地理情報付きのデータがあれば、後はGISを準備すれば地図上へのマッピングから集計までほぼ自動化できる。
　GISとしてはOSSとして配布されているQGISを用いることとした。
　QGISはライセンスを守る限り自由に利用でき総務省が使用法を紹介するほか、国土交通省などもお薦めのGISであり、大学などのユーザーも多い²。

事前準備

　QGISを用いて推定を行うに当たり、次のものを準備した。

1. 背景地図（地理院タイルをmbtiles形式でダウンロードしたもの）
2. 国土地理院の基盤地図情報
3. 総務省統計局のメッシュ統計データと境界データ
   　１の背景地図は表示されるデータが正しいか人間が見て判断するためのものである。地理院タイルはラスタ形式で紙の地形図と同じような見た目で表示できる。これをダウンロードするのは手間なので２の基盤地図情報をセットアップして代替しても良いが、こちらはベクタ形式で見た目が素っ気なくわかりにくい。１と２については、自宅で準備したものを使った。
   　なお、３が今回最も大切なデータである。

QGISでの作業

地理情報の登録・表示

　QGISでの作業は地理情報をグループ毎にレイヤに追加することから始まる。
　まず背景地図をレイヤに追加する。
　次に250mメッシュをレイヤに追加する。
　最後に被災世帯のデータはCSVのため、そのままでは地図上に表現できないので、次の手順で追加する。

1. メニューの【レイヤ】から【レイヤを追加】→【CSVテキストレイヤを追加】を選ぶ。
   2.「 データソースマネージャ|CSVテキスト」というダイアログが開くので登録するファイル名を選択し、レイヤ名を入力する。ファイル選択は【ファイル名】のテキストボックスの右隣にある「・・・」のボタンをクリックするとファイル選択のダイアログが開く。
2. 【ファイル形式】は変更不要、【ジオメトリ定義】ではポイント座標を選択し「X値」に経度、「Y値」に緯度を示す項目を指定する。
3. このデータは単純に緯度・経度の情報があるだけなので「ジオメトリのCRS」はWGS84のまま変更しない。
4. 【追加】をクリックする。
   　これでCSV形式の地理情報を地図上のポイントとして追加できる。
   　ちなみに地理院地図タイルをダウンロードしたものはCRSがEPSG3857(WGS84/Pseude-Mercator)という形式で投影されているが、先ほど準備した２と３はEPSG2452(JGD2000/Japan Plane Rectangular CS X)という方法で投影されていて、本来は地図の表現形式が異なるため重ねることができない。
   　しかし、QGISにはオンザフライ投影機能という便利な機能があり、異なる投影方式の地図を自動で重ねることができる。
   　なおざっくり説明するとEPSG3857は極地を除く地球全体を表現できるが、EPSG2452は日本が定めた公共座標系の第10系を意味し、東北地方を第10系の基準点からの距離で表現できる。
   　

世帯数集計対象メッシュの選択と集計

　正しく被災世帯のデータを取り込めると、地図上に○で表示される。今回はこの情報についてはダミーデータを使用する。次の地図上の赤いポイントが氏被災世帯と考えてほしい。

　世帯数集計対象のメッシュは、メッシュ毎に内部に含まれる○の数を数え、１以上のメッシュを選択すれば良い（条件は世帯数 > 0とする）。
こんな感じ。

ポイントが含まれるメッシュを色づけすると次のようになる。

　QGISで集計対象のメッシュを選択したレイヤを作成しshapefileに保存すると、メッシュに割り振られたコードがdbfファイルに記録されているので、MS AccessなどRDBMSで取り込み、公開されているCSVとリンクすれば集計ができる。

　QGISで、この集計を行うには次のとおり操作する。
　QGISだけで最後まで集計できるが、複数の表を結合処理のはMS Accessのほうがやりやすいので、QGISでは対象メッシュの選択まで行う。
　図については本当のデータを使えないのでダミーデータを用いる（秋田市役所の周囲にポリゴンを設定し、その内部にランダムなポイントを設定した）。

1. メニューの【ベクタ】から【解析ツール】→【ポリゴン内の点の数】を選択する。
2. 「ポリゴン内の点の数」というダイアログが開く。【ポリゴン】に250mメッシュのレイヤを【ポイント】に被災世帯のレイヤを指定する。
3. 実行すると「カウント（count）」というレイヤが新しく作成される。これは250mメッシュのレイヤに内部の点を計数した項目が追加されている。追加された項目目は先のダイアログの【カウント属性名】で指定されたものとなる。
4. このままでは使いにくいので、「被災世帯数 > 0」のメッシュを抽出したレイヤを作る。まず、レイヤパネルで追加された「カウント(count)」レイヤを選択し、マウスの右ボタンをクリックし、【属性テーブルを開く】を実行すると属性テーブルのダイアログが開く。
5. ダイアログのメニューバーにΣのアイコン（QGIS式による選択）があるのでクリックすると、式を入力するダイアログが開く。左側に文字入力する領域（式入力ペイン）があるので「"NUMPOINTS" > 0」と入力する。これは先の集計結果の項目NUMPOINTSが0より大きいレコードを選択することを意味する。ついで右下にある【地物を選択】というコンボリストをクリックすると対象レコードが選択される。
6. 選択レコードのみの新しいレイヤを作成する。レイヤを選択したままマウス右ボタンをクリックし【エクスポート】を選び【新規ファイルに選択地物を保存】をクリックする。ダイアログが開くので、保存するファイル名を指定する。shapefile形式で出力されたファイルはいくつかあり、地物の属性はdBase形式のファイルに保存されている。MS Accessでこのファイルを読込、メッシュ統計の属性値ファイル（CSVファイル）と紐付ければ人口、世帯数など好きな項目を集計できる。この作業はMS Accessの基本的な操作のため省略する。
   　以上の作業で推定すると、約24,000世帯ほどだった。地図上で250mメッシュと被災世帯の地理的分布を見ると、メッシュが荒すぎる感じがしたため、250mメッシュを50mに分割したメッシュをmieruneが公開しているツールを用いて作成し推定してみたが、結果はあまり変わらなかった。

総括など

　GISを用いた被災世帯数の推定作業は、必要なデータの準備が大変だが、考え方はわかりやすく応用も広い。今まで地理的な広がりも含めて統計情報をまとめることは大変だったがGISが使えることで実用的な手法になったと思う。
　例えば学校や何かの施設の位置を元に距離や領域のポリゴンデータを準備することで、影響を受ける人口等の推計も行える。
　一方GISを使うこと自体が、あまり一般的では無いため対応できる人員が限られることは難点と言える。
　また、被災世帯情報をマッピングしてみると明らかに間違った地点を示しているデータが散見された。これは被災情報を集約するための某大企業製システムで地理情報を扱えるが、テキストの住所情報と地理情報が連動しておらず、手入力で地理情報を入力したことが原因と考えられる。
　実際に罹災証明書の発行数は約7500件だった。推定世帯数と大幅に異なる原因は字（街区）全体が浸水した地域よりも、一部が周囲と比べて特に水深が深くなったところが多かったことが考えられる。
　窓口対応で被害状況を聞き取りしていると、自宅敷地が周囲よりも低いため、周囲の家から水が流れ込んできて自宅だけ浸水したという話を何度もきいた。地形や盛り土、排水路などの状況で、そこの一角だけ浸水したというところが意外と多い印象がある。

感謝

　ポリゴン内部のポイントを数える方法は、共創プラットフォームで名古屋市役所の方から教えて頂いた。

1. 某大手企業製システムからの出力ファイルだが、表頭の名称内に改行文字が含まれ、そのままではデータ加工できない邪悪な作りだった。
   　この被災世帯データを元に支援業務を行うが、データ項目を見ると緯度・経度形式の地理情報が含まれていた。 ↩

2. 稼働環境も多い。MS Windowsの他、Unixクローン環境なら、ほぼ使うことができる。私はx86ではFreeBSDやLinux、Macで使用し、さらにSBCのRaspberry pi3Bでも動作を確認している。無償で利用でき、さらに使用方法がわかりやすいため、地理情報を共有する上で強力なツールと考えている。 ↩