フューチャー Advent Calendar 2020 17日目です。
昨日は @alfe_below さんによる「package.json dependencies メンテの仕方最短ルート」でした。

はじめに

地理情報を扱うアプリケーションを開発するにあたって、基礎知識として頭に入れておいたほうが良さそうなことと、PostGISの簡単な使い方をまとめてみました。

Google Maps JavaScript APIやLeaflet.jsなど、フロントエンドのライブラリの使い方については触れていません。また、測量に関する細かい話も書いていません。

基礎知識

ここでは、地理情報のデータとしての表現方法や、システムの設計・開発にあたって覚えておいたほうが良さそうなことを簡単に説明します。

そもそも地理情報とは？ということについては、以下のリンク先が参考になるかと思います。

地理情報とGISデータモデル
https://www.esrij.com/gis-guide/gis-datamodel/gis-datamodel/

CRSとEPSGコード

(139.767, 35.681) という点があったとして、それが地球のどこを指し示すのかは原点や単位系といったパラメータの取り方によって異なるものになります。

このパラメータ群をまとめてCRS（空間参照系）と呼び、その組み合わせについてEPSG（European Petroleum Survey Group）という団体がユニークなIDを採番したものがEPSGコードです。よく使われているものの例としてEPSG:4326があり、GPSで使われています。

CRSは計算によって変換（射影）が可能で、各言語でライブラリが開発されています。

参考リンク

EPSG:4326 仕様
https://epsg.io/4326
proj4js : JavascriptのCRS変換ライブラリ
https://github.com/proj4js/proj4js
pyproj : PythonのCRS変換ライブラリ
https://github.com/pyproj4/pyproj

Well-Known Text (WKT)

駅の位置、道路の始点と終点、土地の形状など地理情報といった情報はベクトルによって表現されます。たとえば駅の位置であれば、2次元のベクトルが1個として表現できます。

Well-Known Textは、これらベクトルの表記方法です。例えば駅の位置（点）であれば、POINT(139.767 35.681) というように書くことが出来ます。

テキストではなくバイナリデータで表す場合、Well-Known Binary (WKB)という表記方法があります。

参考リンク

Wikipediaでの解説
https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry

各種データフォーマット

地理情報オブジェクトを保存する際のフォーマットにはいくつか種類があります。

Shapefile

地理情報システムの老舗であるESRIが開発した形式です。仕様も公開されており、デファクトスタンダードと言って良いかと思います。

また政府発行の地理情報データも、だいたいShapefile形式での提供があります。

一方シェープ”ファイル”と言いながらも、拡張子違いで最低3つのファイルを同じディレクトリに配置して読み込み、という独特の扱いにくさがあります。

参考リンク

シェープファイルとは
https://desktop.arcgis.com/ja/arcmap/10.3/manage-data/shapefiles/what-is-a-shapefile.htm
政府発行の町・字境界データ。Shapefileでの提供があります。
https://www.e-stat.go.jp/gis/statmap-search?type=2

KML

Google Earth、Google Mapsで主に使われる形式です。XMLです。

参考リンク

KMLとは
https://desktop.arcgis.com/ja/arcmap/10.3/manage-data/kml/what-is-kml-.htm

GeoJSON

その名の通りJSON形式で地理情報データを表現します。その仕様はRFC7946で定められています。

Web Friendlyなフォーマットで、たとえばLeaflet.jsではGeoJSONオブジェクトを簡単にレンダリング出来るようになっています。
https://leafletjs.com/examples/geojson/

他にも多くのサービス、ツール類がGeoJSONに対応しています。

参考リンク

公式サイト
https://geojson.org/
GeoJSONを扱えるサービスやライブラリなどのまとめ
https://github.com/tmcw/awesome-geojson

位置情報のハッシュ化

GeohashやQuadkeyといったアルゴリズムを使用すると、位置情報をハッシュ化することが可能です。

位置情報を扱う上で近傍検索というのはよくあるユースケースだと思いますが、これを緯度・経度から1件ずつ距離を計算していては計算コストが非常に高くなってしまいます。位置情報をハッシュ化しておくと、これが文字列の前方一致検索で可能になります。

参考リンク

位置情報データの扱い方（ジオハッシュとか）
https://qiita.com/yabooun/items/da59e47d61ddff141f0c

標準地域メッシュ

国が地域を緯度・経度に基づいて分割したものです。日本ではその範囲の広さに応じて1次メッシュ、2次メッシュ、3次メッシュ…がJISで定められています。

緯度・経度から計算でメッシュコードが求められ、また近傍のメッシュは前方一致するので、ハッシュコードのような扱い方も可能です。

参考リンク

標準地域メッシュ
https://www.esrij.com/gis-guide/gis-other/mesh/
地域メッシュ統計について
http://www.stat.go.jp/data/mesh/m_tuite.html
緯度経度からメッシュコードを作成する方法
http://white-bear.info/archives/1400

PostGIS

PostGISは、PostgreSQLで地理情報を扱えるようにするための拡張機能です。PostGISを有効化すると、いくつかのデータ型と関数が使えるようになります。

公式サイトで紹介されている次のクエリでは、ゴッサム・シティに所在するスーパーヒーローの名前を取得することが出来ます。

SELECT superhero.name
FROM city, superhero
WHERE ST_Contains(city.geom, superhero.geom)
AND city.name = 'Gotham';

ST_Contains() は、第1引数のオブジェクトが第2引数のオブジェクトを内包する場合にTrueを返します。はみ出していても被っていればTrueとしたい場合はST_WithIn() が使用可能です。他にもPostGISは多数の関数を提供しています。

ST_Contains
https://postgis.net/docs/ST_Contains.html
ST_Within
https://postgis.net/docs/ST_Within.html
PostGIS Reference
https://postgis.net/docs/reference.html

GeometryとGeography

PostGISで扱えるデータ型に、GeometryとGeographyというふたつの型があります。これらの違いは、その座標をデカルト座標として扱うか否かです。

緯度と経度は角度ですので、Geographyとして扱うのが適切です。しかしGeometry型に緯度経度を入れることも可能で、用途次第では特に問題は生じず計算コストが低いというメリットがあります。

公式ドキュメントでは、「正直どっちがいいかわからないならGeographyがいいよ」とされています。

参考リンク

Geography
https://postgis.net/workshops/postgis-intro/geography.html
When to use Geography Data type over Geometry data type
https://postgis.net/docs/using_postgis_dbmanagement.html#PostGIS_GeographyVSGeometry

とりあえず使ってみる

何はともあれPostGISを動かして見るときはDockerを使うと便利です。
postgis/postgisというイメージがDockerHubに登録されているのでこれを使います。

$ docker run -ti --rm -p 15432:5432 -e POSTGRES_PASSWORD='postgres' postgis/postgis:11-2.5

このイメージの使用方法はpostgresとほぼ同じです。
コンテナが立ち上がれば、あとはお好みのクライアントから接続してください。

DBeaverというクライアントを使用すると、Geography型のデータをOpenStreetMap上でプレビューしてくれてとても便利です。

InsertとSelectの例

空港を保存するテーブルを例として考えます。ここではGeography型を使います。

create table airports (
    id serial primary key
    , name text
    , iata_code text
    , location geography(POINT, 4326)
);

このテーブルに羽田空港、成田空港を登録するときは、次のようなINSERT文になります。

insert into airports (name, iata_code, location)
values
    ('東京国際空港', 'HND', ST_GeogFromText('SRID=4326;POINT(139.7776446 35.5493975)'))
    , ('成田空港', 'NRT', ST_GeogFromText('SRID=4326;POINT(140.3906614 35.7719867)'));

Geography(Geometry)型のカラムにデータを投入する際は、ST_GeogFromText(ST_GeomFromText)のような関数を使用します。この関数によって、EWKT形式のテキストからデータをInsertすることが出来ます。

羽田空港（HND）から最も近い空港を取得するクエリは次のようになります。

select
    airports.iata_code
    , ST_AsText(airports.location) as location
    , ST_Distance(airports.location, ref_airport.location, false) as distance
from
    airports
    , (
        select location
        from airports
        where iata_code = 'HND'
        limit 1
    ) ref_airport
where
    airports.iata_code != 'HND'
order by
    distance asc
limit 1

結果は次のようになります。Geography型を使っている場合、distanceの単位はメートルとなります。地理院地図で羽田空港と成田空港の直線距離を計測したところ約60kmだったので、正しく計算が出来ていると言えます。

iata_code|location                      |distance      |
---------|------------------------------|--------------|
NRT      |POINT (140.3906614 35.7719867)|60661.39893749|

空間インデックス

PostGISではGeometry / Geography型のカラムにインデックスを作成することが出来ます。先の羽田から近い空港を探すようなクエリであれば、locationカラムに空間インデックスを貼ることでパフォーマンスの向上が見込めます。

create index gix_airports_location on airports using gist(location)

Pythonとの連携

PostgreSQLを扱うpsycopg2、ジオメトリオブジェクトを扱うShapely、Pandasのようにジオメトリオブジェクトを含むデータを扱えるGeoPandasという組み合わせを紹介します。

Shapelyを使って点を定義

from shapely.geometry import Point

p = Point(139.7776446, 35.5493975)

psycopg2でデータベースに保存

import psycopg2

conn = psycopg2.connect("postgresql://...")

with conn:
    sql = """
    insert into airports (name, iata_code, location)
    values (%s, %s, ST_GeogFromText(%s))
    """
    with conn.cursor() as cur:
        cur.execute(sql, ("東京国際空港", "HND", "SRID=4326;" + p.wkt))

conn.close()

データベースからGeoPandasのDataFrameにデータをロード

import geopandas as gpd

sql = """
    select id, iata_code, location from airports
"""
conn = psycopg2.connect(_dsn)
gdf = gpd.read_postgis(sql, conn, geom_col="location")
_logger.info(gdf)

#    id iata_code                    location
# 0   1       HND  POINT (139.77764 35.54940)
# 1   2       NRT  POINT (140.39066 35.77199)

GeoPandas.DataFrameからGeoJSONへ変換

gdf.to_json(indent=2)

# {
#   "type": "FeatureCollection",
#   "features": [
#     {
#       "id": "0",
#       "type": "Feature",
#       "properties": {
#         "iata_code": "HND",
#         "id": 1
#       },
#       "geometry": {
#         "type": "Point",
#         "coordinates": [
#           139.7776446,
#           35.5493975
#         ]
#       }
#     },
#     ...
#   ]
# }

参考リンク

psycopg2
https://pypi.org/project/psycopg2/
Shapely
https://pypi.org/project/Shapely/
GeoPandas
https://geopandas.org/

おわりに

地理情報データおよび測量はとっても奥が深い世界でした。アリストテレスの時代から知見が積み重ねられてきた分野なので当然といえば当然です。

この記事がこれから位置情報を扱うアプリケーションを開発する誰かの役に立てば幸いです。