はじめに

場所的な情報、量的な情報、両者の時間的な情報、それら全てを考慮する必要があるため、GISデータ、人流データ等と呼ばれるデータは、少々扱いが特殊だと感じました。

そのため、自身の学習のノートとしてまとめます。（2024-02-26時点）

また、倫理的配慮から、匿名化されたデータであることが多く、別なデータと一意に結びつける事が難しいため、有用に扱うには、それらを加味した上で有効な方法を検討する必要がありそうです。

私が調べた限り、人流データと呼ばれるデータにおいて、既存の調査・研究等では、何かしら目的に適した「前処理」をしている（せずに扱う手法を見たことがないです）

たとえば以下のようなもの

「地理」的な情報を欠落させ、場所はあくまでカテゴリ（名義尺度）として扱う
（仮想的な）個人を特定し「時間」を欠落させ、（移動の）軌跡を抽出する（a→d→b→cの様な）
「時間・空間」で丸め、存在する（匿名化され、仮想的な）人数や属性データを扱う

扱う事柄、目的を意識せずしては何も出来ない。見たいものを見ることも容易ではないようです。

また、統計で多群を比較する際は留意が必要なようです。

統計検定を理解せずに使っている人のために III

参考：

度々具体的なデータを扱いますが、出処として、1kmメッシュ別将来推計人口データから取得した、宮城の1kmメッシュ別将来推計人口（H30国政局推計を使用しています。

話はメッシュ形式のデータを前提に進みます。

# e.g. データを取ってくる

import shutil
import requests

r = requests.get('https://nlftp.mlit.go.jp/ksj/gml/data/m1kh30/m1kh30-18/1km_mesh_suikei_2018_shape_04.zip')
with open('1km_mesh_suikei_2018_shape_04.zip', 'wb') as f:
  f.write(r.content)

# データの解凍
shutil.unpack_archive('1km_mesh_suikei_2018_shape_04.zip')

# e.g. データをGeoJSONに変換する
$ docker run --rm -v `pwd`:/WORK ghcr.io/osgeo/gdal:ubuntu-full-latest ogr2ogr -f GeoJSON /WORK/foobar.geojson /WORK/1km_mesh_2018_04.shp

GISデータとは
データの読み込み・加工
- GDAL
- geopandas
Plotlyによる可視化
参考

GISデータとは

GISデータとは、場所と、それに対応する任意のデータが結びついたもの全般を指す言葉のようです。

人流データというと時間や移動そのものが絡んでくることが多いようです。

GISとは

地理情報システム（GIS：Geographic Information System）は、地理的位置を手がかりに、位置に関する情報を持ったデータ（空間データ）を総合的に管理・加工し、視覚的に表示し、高度な分析や迅速な判断を可能にする技術である。
引用: 国土地理院 - GISとは・・・

地理空間情報とは

地理空間情報とは、空間上の特定の地点又は区域の位置を示す情報（位置情報）とそれに関連付けられた様々な事象に関する情報、もしくは位置情報のみからなる情報をいう。
引用: 国土地理院 - GISとは・・・

人流データとは

人流データとは、人がどのように動いたかを数値化したものです。
https://www.ntt.com/business/sdpf/knowledge/archive_94.html

GISデータに何かしらの形で時間を意味するデータがくっついた形

おおよそ知られているフォーマット

GPS等で得られたヒト、モノ（ID）等の移動軌跡そのもの
設置されたセンサー等のポイントに対して、得られた計測値
1.から特定エリア内で集計されたデータ
1.から生成された一定時間内の起点終点データ

に加え、計測対象（人の属性）や計測地点（場所の属性）を付加する場合がある。

そのままで処理するのはかなり大変？

メッシュについて

GIS・人流データ関連において、「メッシュ」というと、地域標準メッシュと呼ばれるだいたい1km四方のメッシュを指すことが多いです。

一方で、独自のメッシュデータで100m以下のメッシュ情報を提供する業者も存在するようです。

メッシュデータとは、地図上の情報をデジタル化したり各種統計情報をとるために地図上の経緯度方眼として定められた地域メッシュのことです。
国土数値情報のメッシュデータは、総務省（旧総務庁）が定めた「統計に用いる標準地域メッシュおよび標準地域メッシュコード」に従って、それぞれの区域に関する統計データを編集したものです。
引用: 国土数値情報ダウンロードサイト - メッシュデータについて

地域標準メッシュコードの一覧は下記から入手できます。

総務省統計局 - 市区町村別メッシュ・コード一覧

※「だいたい」なのは、地球が球であるが、それを平面とみなして処理しようとしている弊害。

この辺を参考QGIS Documentation - やさしい GIS 入門 8. 座標参照系

主な地域メッシュの区分方

地域メッシュの区分方法
区画の種類	区分方法	緯度の間隔	経度の間隔	一辺の長さ	地図との関係
第1次地域区画	全国の地域を偶数緯度及びその間隔（120分）を3等分した緯度における緯線並びに1度ごとの経線とによって分割してできる区域	40分	1度	約80km	20万分の1地勢図の1図葉の区画
第2次地域区画	第1次地域区画を緯線方向及び経線方向に8等分してできる区域	5分	7分	約10km	2万5千分の1地勢図の1図葉の区画
（統合地域メッシュ）			30秒
基準地域メッシュ	第2次地域区画を緯線方向及び経線方向に10等分してできる区域	30秒	45秒	約1km
（第3次地域区画）
2分の1地域メッシュ	基準地域メッシュ（第3次地域区画）を緯線方向、経線方向に2等分してできる区域	15秒	22.5秒	約500m
（分割地域メッシュ）
4分の1地域メッシュ	2分の1地域メッシュを緯線方向、経線方向に2等分してできる区域	7.5秒	11.25秒	約250m
（分割地域メッシュ）

参考：

位置を表すデータの大別

wikipedia - GeoJSON ジオメトリの項に分かりやすくまとめられている。

点・線・ポリゴン（多角形）で表される。

大体、一意なIDが振られていて、対応するデータと合わせ利用される。

メッシュデータの場合、メッシュ一つ一つに対応するIDが振られており、メッシュコードをキーにし、地図上の位置的な情報と、目的の数値データを対応付けている。

位置を表すデータと、それに対応するデータ両方が含まれるデータもあれば、別々に与えられるデータもある。

引用: visualizing.jp - GeoJsonとは

参考：

GeoJSON

御託を並べるより見た方が早いので、サンプルとして宮城県の人口推移予測データGeoJSONに変換しものをみてみる。

以下の様にfeatures要素以下にリストとして、辞書が存在し、各辞書は、properties下の任意の情報と、geometry下に、座標データを持つ。

このデータは、位置を表すデータと、それに対応するデータ両方が含まれている例。

{
    "type": "FeatureCollection",
    "name": "1km_mesh_2018_04",
    "crs": {
        "type": "name",
        "properties": {
            "name": "urn:ogc:def:crs:EPSG::4612"
        }
    },
    "features": [
        {
            "type": "Feature",
            "properties": {
                "MESH_ID": 56405557,
                "SHICODE": 4341,
                "PTN_2015": 3.9816,
                "HITOKU2020": "*",
                "GASSAN2020": 56405556,
                "PTN_2020": 2.8717,
                "PT0_2020": 0.0,
...etc...
                "RTE_2050": 0.0
            },
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [
                        [
                            140.7125,
                            37.791666666666664
                        ],
                        [
                            140.7125,
                            37.8
                        ],
                        [
                            140.725,
                            37.8
                        ],
                        [
                            140.725,
                            37.791666666666664
                        ],
                        [
                            140.7125,
                            37.791666666666664
                        ]
                    ]
                ]
            }
        },
        {
            "type": "Feature",
            "properties": {
                "MESH_ID": 56405558,
                "SHICODE": 4341,
...etc...

参考:

座標系についての簡単な説明

※「だいたい」なのは、地球が球であるが、それを平面とみなして処理しようとしている弊害。
この辺を参考QGIS Documentation - やさしい GIS 入門 8. 座標参照系

ものすごくさらっと流しているけどとても大切なこと。

また、2点間・ラインストリングの距離を出したいときなどにも変換が必要になる

TODO write

日本（政府公開データで）採用している座標系
google mapで採用している座標系　WGS84
mapboxで採用している座標系

データの読み込み・加工

目的に沿った適切なデータに前処理するのにも、使いたいライブラリで扱いやすい形にするのにも役に立つ

GDAL

GDAL is a translator library for raster and vector geospatial data formats that is released under an MIT style Open Source License by the Open Source Geospatial Foundation.
引用: GDAL

ものすごく多機能である一方、数字に用がある場合使うのはVector Programsぐらい？

コンテナイメージが配布されているので、それを使うのが便利

ファイルの変換等

geopandas

GeoPandas is an open source project to make working with geospatial data in python easier
引用: https://geopandas.org/en/stable/

install

pip install geopandas

以下のようにするだけで、ものすごく多彩なデータをそのまま、pandasの様に扱える

import geopandas
gdf = geopandas.read_file("path_to_data")

参考:

GeoPandas - Documentation

Plotlyによる可視化

以の様な書き方で画像の様なplotlyによる地図・メッシュの描写が出来る。

データはgeopandasで読みこませると、GeoJSON、.shp(シェイプファイル)でも一様に扱う事ができる様子。

import geopandas as gps
import plotly.express as px

pop = gps.read_file('1km_mesh_2018_04.shp')

fig = px.choropleth_mapbox(pop[['MESH_ID','PT0_2050']], # 描写するデータフレーム
                           geojson=pop.geometry,# 描写の元となるポリゴン等のデータ
                           locations=pop.index, # 描写するデータに含まれた位置を一意に定めるID
                           color='PT0_2050',    # 塗りつぶしに使うデータフレームの列名
                           color_continuous_scale="Viridis",
                           #featureidkey='properties.MESH_ID',# GeoJSONのポリゴンに振られた一意なIDデフォルト値'id'、大体properties以下
                           range_color=(500, 1500), #色付けの範囲
                           mapbox_style="carto-positron",
                           #color_discrete_map={},
                           zoom=7, 
                           center = {"lon": 140.8693, "lat": 38.2682},
                           opacity=0.3,
                           labels={'PT0_2050':'2050年の人口予測'}
                          )
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.show()

下のように書くことも出来る。

引数、geojson,locations,featureidkeyに留意。

If the GeoJSON you are using either does not have an id field or you wish you use one of the keys in the properties field, you may use the featureidkey parameter to specify where to match the values of locations.

引用: Mapbox Choropleth Maps in Python - Indexing by GeoJSON Properties

fig = px.choropleth_mapbox(pop[['MESH_ID','PT0_2050']], # 描写するデータフレーム
                           geojson=pop, 
                           locations='MESH_ID', # 描写するデータに含まれた位置を一意に定めるIDの列名
                           color='PT0_2050',    # 塗りつぶしに使うデータフレームの列名
                           color_continuous_scale="Viridis",
                           featureidkey='properties.MESH_ID',# GeoJSONのポリゴンに振られた一意なIDデフォルト値'id'、大体properties以下
                           range_color=(500, 1500), #色付けの範囲
                           mapbox_style="carto-positron",
                           zoom=7, 
                           center = {"lon": 140.8693, "lat": 38.2682},
                           opacity=0.3,
                           labels={'PT0_2050':'2050年の人口予測'}
                          )
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.show()

例えば、人数ごとに区分を作るなどする場合を考える。

地図の描写ライブラリ側で、値に対応する色等を定義することもできる（ライブライもある）。

下のように、pandasのメゾットを使い、区分分けした列を作成し、それを与えることも考えられる。

pd.cut(pop['PT0_2050'],bins=[50,100,500,1000,2000,pop[candidate].max()],labels=False)

参考:

参考

GIS・人流データに関して

統計

Plotly & dash

メモ

以下完全にメモ

Plotlyでアニメーションなコロプレス地図を作る

新規感染者数（県/日）のデータから１週間合計のデータを作って、地図上にプロットする。

合計したり、日数削ってるのは重いから。

データからわかる－新型コロナウイルス感染症情報
から新規陽性者数の推移（日別）のデータ
japonyol.net - 47都道府県のポリゴンデータ geojson
から県境のポリゴンデータを取ってきて遊んでみる。

データの用意

!wget https://japonyol.net/editor/article/geo/prefectures.geojson
!wget https://covid19.mhlw.go.jp/public/opendata/newly_confirmed_cases_daily.csv
import json
import pandas as pd
import geopandas as gps
import plotly.express as px

# geodataframeは読み混んでるけど実際は前処理にしか使っていない
pref = gps.read_file('prefectures.geojson')
df   = pd.read_csv('newly_confirmed_cases_daily.csv')

with open('prefectures.geojson') as f:
    geojson = json.load(f)

若干面倒な所として、列としてフレームを与えなければならない。

GeoJSON, geodataframeで扱う場合、同じポリゴンで、日付が違う要素を複数作る必要があるのでちょっと面倒。

そのため、DataFrameと県境を表すGeoJSONは別に与える。

日付	県名1	県名2
2022/8/1	0	1
2022/8/2	1	2

上を下の形に変換している

日付	県名	値
2022/8/1	県名1	0
2022/8/1	県名2	1
2022/8/2	県名1	1
2022/8/2	県名2	2

データの整形

# データの整形
l =  []
df['Date'] = pd.to_datetime(df['Date'])
df = df.set_index('Date')

# 週ごとに集計する
tempdf = df.resample('w').agg('sum')

# 列の数だけ繰り返す
for c in tempdf.iloc[:,1:]:
    tdf = pd.DataFrame(
        {
          'pref': c,
          'value': tempdf[c]
        },
          index=tdf.index
    )
    l.append(tdf)

adf = pd.concat(l).sort_index()
adf

# Hokkaido -> 1 のような県名（アルファベット）から整数に変換

pref['alpha']=df.T.drop('ALL').reset_index()['index']
d = pref['alpha'].to_dict()
# k,v = v,kになっていることに留意
pref_num_map = {v: k for k, v in d.items()}

# 期間でスライスする(resampleしたときにやれというツッコミはさておき)
pdf = adf.loc['2022-8-1':'2022-9-1'].reset_index()
# アニメーションフレームのキーにするのにdateのままだと都合が悪いので
pdf['Date']=pdf['Date'].astype(str)
# print(pdf)
# Hokkaido -> 1 のような県名（アルファベット）から整数に変換
pdf['pref']=pdf['pref'].apply(lambda x: pref_num_map[x])
pdf

描写

# 描写処理

px.choropleth_mapbox(
    pdf,
    locations='pref',
    featureidkey="properties.pref",
    geojson=geojson,
    color='value',
    range_color=[100,300000],
    mapbox_style="carto-positron",
    center = {"lon": 139.6917, "lat": 35.6894},
    animation_frame='Date',
    zoom=5,
    height=900
).show()

表現する期間や幅、色味を調節しないと何を表したいかよくわからない
（アニメーションも簡単に作れると指名したい）

とても重い（処理時間がかかる・データ量が多い）

dashに無理やり人口動態統計を描写した時のコード

動けばよろしい、とりあえず画面さえできれば良いというというノリ

読み込みに時間がかかり、long_callbackにしないとタイムアウトした気がする。

long_callbackデコレータを使うには、マネージャーを用意する必要があり、CeleryかDiskCacheを使う。

参考： PythonのウェブフレームワークDashのver2.0に関して

from dash import Dash, dcc, html, Input, Output
import plotly.express as px
import geopandas as gps
import pandas as pd
from dash.long_callback import DiskcacheLongCallbackManager
from dash.dependencies import Input, Output

## Diskcache
import diskcache
cache = diskcache.Cache("./cache")
long_callback_manager = DiskcacheLongCallbackManager(cache)

app = Dash(__name__)


token = 'mapboxtoken'
pop = gps.read_file('1km_mesh_2018_04.shp')
# 総人口だけ抜き出す
p = [n for n in pop.columns if 'PT0' in n]
print(pop)
print(p)



app.layout = html.Div([
    html.H4('人口動態統計'),
    html.P("Select year:"),
    dcc.RadioItems(
        id='mapbox-county-choropleth-x-candidate', 
        options=p,
        value="PT0_2050",
        inline=True
    ),
    dcc.Graph(id="mapbox-county-choropleth-x-graph"),
    dcc.Graph(id='mapbox-county-histgram-x-graph')
])


#@app.callback(
@app.long_callback(
    Output("mapbox-county-choropleth-x-graph", "figure"), 
    Input("mapbox-county-choropleth-x-candidate", "value"),
    manager=long_callback_manager,
    )
def display_choropleth(candidate):

    pop['tmp']=pd.cut(pop[candidate],bins=[50,100,500,1000,2000,pop[candidate].max()],labels=False)
    fig = px.choropleth_mapbox(
        pop[p+['tmp','MESH_ID']], 
        geojson=pop.geometry, 
        color='tmp',
        locations=pop.index, 
        mapbox_style="carto-positron",
        zoom=8, 
        center = {"lon": 140.8703, "lat": 38.2682},
        opacity=0.3,
        hover_data=candidate,
        width=1000,
        height=600
        )
    #fig.update_xaxes(type="log")
    fig.update_layout(
        margin={"r":0,"t":0,"l":0,"b":0},
        mapbox_accesstoken=token
    )

    return fig

@app.long_callback(
    Output("mapbox-county-histgram-x-graph", "figure"), 
    Input("mapbox-county-choropleth-x-candidate", "value"),
    manager=long_callback_manager,
    )
def display_histgram(candidate):
    pop['tmp2']=pd.cut(pop[candidate],bins=[50,100,500,1000,2000,pop[candidate].max()]).astype(str)
    fig = px.histogram(pop['tmp2'])
    fig.update_yaxes(range=(0,2700))
    return fig

if __name__ == "__main__":
    app.run_server(debug=True)