地理情報・PostGISことはじめ

[フューチャー Advent Calendar 2020](https://qiita.com/advent-calendar/2020/future) 17日目です。
昨日は @alfe_below さんによる「[package.json dependencies メンテの仕方 最短ルート](https://qiita.com/alfe_below/items/1141ec9acbb81b504855)」でした。

## はじめに

地理情報を扱うアプリケーションを開発するにあたって、基礎知識として頭に入れておいたほうが良さそうなことと、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](https://en.wikipedia.org/wiki/Geohash)や[Quadkey](https://docs.microsoft.com/en-us/bingmaps/articles/bing-maps-tile-system)といったアルゴリズムを使用すると、位置情報をハッシュ化することが可能です。

位置情報を扱う上で近傍検索というのはよくあるユースケースだと思いますが、これを緯度・経度から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](https://postgis.net/)は、PostgreSQLで地理情報を扱えるようにするための拡張機能です。PostGISを有効化すると、いくつかのデータ型と関数が使えるようになります。

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

```sql
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](https://hub.docker.com/r/postgis/postgis)というイメージがDockerHubに登録されているのでこれを使います。

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

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

[DBeaver](https://dbeaver.io/)というクライアントを使用すると、Geography型のデータをOpenStreetMap上でプレビューしてくれてとても便利です。

![dbeaver.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/148726/14342417-7a3e-a7c4-7854-27ba7427173e.png)


### InsertとSelectの例

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

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

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

```sql
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）から最も近い空港を取得するクエリは次のようになります。

```sql
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`カラムに空間インデックスを貼ることでパフォーマンスの向上が見込めます。

```sql
create index gix_airports_location on airports using gist(location)
```

### Pythonとの連携

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

Shapelyを使って点を定義

```python
from shapely.geometry import Point

p = Point(139.7776446, 35.5493975)
```

psycopg2でデータベースに保存

```python
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にデータをロード

```python
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へ変換

```python
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/

## おわりに

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

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