PostGISのよく使う機能をまとめた
PostGISの中でも頻繁に使用する機能について調べたので、参照用にSQLとその使用方法についてまとめておきます。
GISは空間参照系や楕円体のあたりなど技術的にもマニアックだが、とりあえず理解している範囲内で基本的なことをメモ。
PostGISについて
PosgreSQLには地理情報に関する機能を集めたPostGIS呼ばれる拡張がある。メリデメは以下の通り。
メリット
- 複雑な集計操作を豊富な関数で処理でき、GIS基盤のミドルウェアとして利用可能。
- データベースに格納できるので既存のデータとの結合やデータの共有が容易。
- 毎回必要に応じて集計すればよいので、粒度や集計単位の異なるshpファイルをたくさん持たなくてもよくなる。
- QGIS等とも相性がよく、手軽にSQLで集計して主題図を記述可能。
デメリット
- 毎回集計するのでメモリに乗るデータなら普通にGISで処理した方が早いかも。
- データベースへのロードやエクスポート作業が必要。
Geometry型
PostGISではGeometry型とよばれるデータ型のカラムを作成し、ジオメトリー操作用の関数(ST_XXXXXX)でこれらを操作し、地理的な検索や集計を行ったり、地図として表示させるポリゴンを取得します。
PostGISには、３つのジオメトリー型が存在する。点（位置、ランドマーク）を表すPOINT、線（道路や川）を表すLINESTRING、面（領域、町丁目）を表すPOLYGONの３つ。厳密にはそれらの複合型、集合型もある。
Geography型
Geometry型だけでなく、Geography型と呼ばれるデータ型もあるが、これは日本国内の地図を触っているだけなら必要なさそうなので無視することにした。
厳密な距離を算出したいような場合でも必要かもしれない。残念ながらGeography型は対応している関数が少ないのでGeometryで持っている方が無難。
必要な時だけGeometry型で持っているデータをGeographyに変換して操作するのがベター。
はまったポイント
PostGISを使う上で初心者である私がはまった注意点は以下の3つです。これから基盤を構築される型は気をつけましょう。
- SRIDを統一する。・・・SRIDとは、2次元座標（数値）を実際の位置と対応させる空間参照系のユニークなIDのことを指します。PostGISでは各ジオメトリー型カラムにSRIDを設定することができる。shp2pgsqlを使ってインポートする際に、SRIDを変換する機能があるのでそれを使ってあらかじめデータベース用の共通SRIDデータのみ含むようにしておくと無難。わからない人はとりあえずSRID＝4326に統一するようにするとよいと思われる。
- 文字コード・・・落としてきたshpファイルが大体Shift-Jisだったりするが、shp2pgsqlでDBへとインポートする際に文字コードを指定できるので、インポートする時点でDBと同じ文字コードに統一しておく。
- Oracleの情報でないか確認する・・・PostGISを触っていると色々と情報ない中、Googleで調べる羽目になるが、"Postgre"をキーワードに入れていても、"Oracle"の情報が入っていたりする。元々SQLが似ているので混同しないように気をつける。実際何度か間違えた。
インストール
インストールは少し古いが、こちらを参照すること。
Stack Builderで適当にやってたら入ると思います。万が一、はまった時は思いっきり悔しがりましょう。
データを作成する
ポリゴンデータのインポート・・・PostGISにはshpファイルをインポートするshp2pgsqlというコマンドがあります。
shp2pgsqlは、shpファイルを引数に取り、SQLを出力します。このSQLを実行することで、データをPostgreSQLに取り込むことが出来ます。
pgAdminでもインポートできるけど、文字コード周りはトラブルの種なので面倒がらずに、こちらを用いたほうが無難。
ここでは国勢調査のデータより、渋谷区の町丁目ポリゴン情報をインポートします。
- shpファイル用テーブルの作成
C:\Users\Test>shp2pgsql -s 4326 -W cp932 -c \Users\Test\Desktop\A00200521201
0DDSWC13113\h22ka13113.shp > createtable.sql
Shapefile type: Polygon
Postgis type: MULTIPOLYGON[2]
createtable.sqlの中身は以下の通り。これで"h22ka13113"にshpファイルと同様のテーブルが作成されます。
SET CLIENT_ENCODING TO UTF8; SET STANDARD_CONFORMING_STRINGS TO ON; BEGIN; CREATE TABLE "h22ka13113" (gid serial, "area" numeric, "perimeter" numeric, "h22ka07_" int4, "h22ka07_id" int4, "ken" varchar(2), "city" varchar(3), ... INSERT INTO "h22ka13113" ("area","perimeter","h22ka07_","h22ka07_id","ken","city",....B37EBD24140C10B47D5B2766140908F9FB0EBD24140'); COMMIT;
- 列の作成・・・テーブルschema.table_stに地理情報データ含むジオメトリ型の列を作成します。
ちなみに、今回のように一つのテーブルに複数の異なるタイプの地理情報型を持たせるのはあまり美しくないのでやめておきましょう。
（今回はテスト用のデータのために、複数のテーブルを作りたくなかったので一つにまとめた。）
SELECT AddGeometryColumn('schema','table_st','column_name_geom1',4326, 'POINT', 2); SELECT AddGeometryColumn('schema','table_st','column_name_geom2',4326, 'LINESTRING', 2); SELECT AddGeometryColumn('schema','table_st','column_name_geom2',4326, 'POLYGON', 2);
ポリゴンデータの作成・・・ポリゴンデータを作成します。
- 点の作成:渋谷駅
INSERT INTO table_st(id, name, column_name_geom1, column_name_geom2,column_name_geom3) VALUES (123, '渋谷駅', ST_GeomFromText('POINT(139.70066 35.65878)', 4326), NULL, NULL);
- 線の作成:LINESSTRING1
INSERT INTO table_st(id, name, column_name_geom1, column_name_geom2,column_name_geom3) VALUES (234, 'LINE-1', NULL,ST_GeomFromText('LINESTRING(139.66657 35.67380, 139.67651 35.67690, 139.68333 35.67981)',4326),NULL);
- 面の作成:POLYGON-1
INSERT INTO table_st(id, name, column_name_geom1, column_name_geom2,column_name_geom3) VALUES (345, 'PLOYGON-1', NULL, NULL, ST_GeomFromText('POLYGON((139.69962 35.67493,139.69615 35.65587,139.70990 35.65822,139.69962 35.67493))',4326));
- 索引の作成・・・ジオメトリーカラムに索引を作成します。索引を使うにはSQLの書き方で工夫することが必要だが、念のため作成しておく。
CREATE INDEX index_geom1 ON public.table_st USING gist (column_name_geom1);
- SRIDの設定・・・列の作成時に全て設定しているためここでは特に必要ないが、列のSRIDの設定に関するSQLは下記の通り。
SELECT UpdateGeometrySRID('schema','table_st','column_name_geom1',4326);
＊類似の関数でST_SetSRIDというのがあるが、これはテーブル定義で持っているSRIDの値を変えるだけで座標系を直接変更するものではないらしいので注意。
- ジオメトリー型のSRIDを変更した情報を取得する。
SELECT UpdateGeometrySRID('schema','table_st','column_name_geom1',4326);
上記のSQLで、作成した結果は以下の通りです。
ポリゴン情報の取得
- 重心の計測・・・ポリゴンから重心を取得します。
SELECT ST_Centroid(column_name_geom2) FROM table_st WHERE id = 234;
- 距離の計測・・・点と点の距離を測る。
SELECT ST_distance(yoyogi.column_name_geom1,shibuya.column_name_geom1) dist, ST_distance(yoyogi.column_name_geom1,shibuya.column_name_geom1, false) meter FROM table_st shibuya, table_st yoyogi WHERE shibuya.name = '渋谷駅' AND yoyogi.name = '代々木八幡駅';
| dist | meter |
| 0.016518376262419 | 1658.68010562005 |
- 面と面の距離を測る。
SELECT ST_distance(p2.column_name_geom3,p1.column_name_geom3) dist, ST_distance(p2.column_name_geom3,p1.column_name_geom3, false) meter FROM table_st p1, table_st p2 WHERE p1.name = 'POLYGON-1' AND p2.name = 'POLYGON-2';
| dist | meter |
| 0.00402964942162027 | 366.006648511634 |
＊値の結果より、面と面の距離は最短距離を計測しているようである。
- 線と面の距離を測る。
SELECT ST_distance(l1.column_name_geom2,p1.column_name_geom3) dist, ST_distance(l1.column_name_geom2,p1.column_name_geom3, false) meter FROM table_st p1, table_st l1 WHERE p1.name = 'POLYGON-1' and l1.name = 'LINE-1';
| dist | meter |
| 0.0170052491895865 | 1568.26832718958 |
＊値の結果より、線と面の距離は最短距離を計測しているようである。
- 単一の線の距離を測る。
SELECT ST_length(column_name_geom2) length, ST_length(column_name_geom2, false) meter FROM table_st WHERE id = 234;
| length | meter |
| 0.0178270688501295 | 1657.57701436593 |
- X、Yの座標値を取得する。
SELECT ST_X(column_name_geom1) X, ST_Y(column_name_geom1) Y FROM table_st WHERE name = '渋谷駅'
| x | y |
| 139.70066 | 35.65878 |
検索と集計
- 共有・・・値の集計を行う際に集計対象となるジオメトリーを指定するためにジオメトリーが重なっているかどうかに基づいて検索をかけます。
POLYGON-1と重複する渋谷区の町丁目ポリゴンを取り出すのが下記のSQLです。ST_Intersectsを使用します。
SELECT moji FROM h22ka13113,table_st WHERE ST_Intersects(column_name_geom3, geom) AND table_st.name = 'POLYGON-1';
上記のSQLの結果であるPOLYGON-1と渋谷区の共有部分のある地域は以下のとおりです。
- 包含・・・値の集計を行う際に集計対象となるジオメトリーを指定するためにジオメトリー間の包含関係に基づき検索を行います。
ここでは渋谷駅が含まれるポリゴンを検索します。ST_Withinを使用します。
geomにcolumn_name_geom1が含まれています。引数の順番に注意してください。
SELECT moji FROM h22ka13113,table_st WHERE ST_Within(column_name_geom1,geom) AND table_st.name = '渋谷駅';
- バッファ検索・・・値の集計を行う際に集計対象となるジオメトリーを指定するためにバッファ検索を行います。
今度は渋谷駅から半径500m以内のバッファに含まれる町丁目を検索します。ST_DWithinを使用します。
column_name_geom1から500m以内の範囲にgeomが含まれています。引数の順番に注意してください。
＊何故か私の環境だと最後のuse_spheroidのbooleanに値を入れないとメートル単位の検索にならず、これはtrueでもfalseでも同様の結果です。
SELECT moji FROM h22ka13113,table_st WHERE ST_DWithin(column_name_geom1,geom, 500, true) AND table_st.name = '渋谷駅';
以下の図は上記の二つのSQLの結果です。
濃い紫が「道玄坂1丁目」で渋谷駅を含む包含関係のSQLの結果です。
ピンクが渋谷駅から500m以内のバッファ検索のSQLの結果です。
- ポリゴンの結合・・・ジオメトリーを結合します。市区町村のポリゴンを都道府県に変更したりする場合などに用います。
GROUP BYの単位で結合されます。下記のSQLで渋谷区全体のポリゴンが返されます。
SELECT st_union(geom) FROM h22ka13113;
- ポリゴンの簡素化・・・町丁目レベルで保有しているような大量のポリゴンを結合すると、一つのポリゴンのサイズが大きくなりすぎて複雑すぎる形になるので簡素化する必要があります。
ポリゴンの簡素化にはST_Simplify、ST_SimplifyPreserveTopologyを用います。前者だと実際に扱うような汚いポリゴンだと、不正なポリゴンが返ってくることがあったので
後者をお勧めします。なお、アルゴリズムにはDouglas–Peucker algorithmが使用されており、簡素化を行う閾値となるパラメータεを指定する必要があります。
パラメータεは大きければ大きいほど簡素化されます。下の図は渋谷区のポリゴンをST_Unionで合成した後に、ST_Simplifyでパラメータを変えて簡素化した際の例です。
SELECT st_union(geom) FROM h22ka13113;#赤：簡素化なし SELECT st_simplify(st_union(geom), 0.0008) FROM h22ka13113;#緑：ε=0.0008 SELECT st_simplify(st_union(geom), 0.002) FROM h22ka13113;#青：ε=0.002
管理
- PostGISのバージョンを問い合わせる
SELECT postgis_full_version();
- 地理情報カラムの情報を問い合わせる
SELECT * FROM geometry_columns;
| postgres | public | table_st | column_name_geom1 | 2 | 4326 | POINT |
| postgres | public | table_st | column_name_geom2 | 2 | 4326 | LINESTRING |
| postgres | public | table_st | column_name_geom3 | 2 | 4326 | POLYGON |
| postgres | public | h22ka13113 | geom | 2 | 4326 | MULTIPOLYGON |
- 空間参照系の情報を問い合わせる・・・ここにIDの存在するSRIDしかPostGISは認識できないので注意、必要なものがなければ追加する。
SELECT * FROM spatial_ref_sys;
| 4324 | EPSG | 4324 | GEOGCS["WGS 72BE",DATUM[... | +proj=longlat +ellps=WGS... |
| 4326 | EPSG | 4326 | GEOGCS["WGS 84",DATUM["]... | +proj=longlat +datum=WGS... |
| 4463 | EPSG | 4463 | GEOGCS["RGSPM06",DATUM["... | +proj=longlat +ellps=GRS... |
| 4470 | EPSG | 4470 | GEOGCS["RGM04",DATUM["Re... | +proj=longlat +ellps=GRS... |
DBからRへとポリゴンデータをダウンロードする
最後に、Rで主題図を書いたりする時に必要になる、DB上に保存されているジオメトリーデータをRに読み込む方法について解説します。
やり方は少し厄介ですが、以下の方法で実現します。
- RからRPostgreSQLを使ってデフォルトドライバーでPostgreSQLにアクセスします。
- PostGISのST_AsText関数でポリゴンを地理情報を文字列で表すWKTに変換してSQLを実行しRにcharacterとして読み込みます。
- これをRのrgeos::readWKTでspクラスに変換します。
このとき、ODBCなんかを使うと文字列長制限で大きなジオメトリーデータはダウンロードできないため、RPostgreSQLを使用しています。
とりあえずもっと賢い方法はあるんでしょうけど、今のところこれで特に不便もしていないのでこれを使用しています。
下記はPostgreSQLに入ってる渋谷区のデータをダウンロードして、人口に関する主題図を作成するRコードです。
library(DBI)
library(RPostgreSQL)
library(rgeos)
#PostgreSQLのドライバーを使用して、DB接続を行う。
conn <- dbConnect(dbDriver("PostgreSQL"),user="postgres",host="localhost", password="password",dbname="dbname")
#文字コードを変更する。環境によっては不要
postgresqlpqExec(conn, "SET CLIENT_ENCODING TO 'SJIS'")
#手もとの環境では、2億文字までは文字列で取得可能であった。20億文字でメモリ不足で落ちた。
#使用する際はあまり文字が大きくなり過ぎないようにST_Simplifyなどを使うように注意する。
nchar(fetch(dbSendQuery(conn, "SELECT repeat('yeah!', 200000000);"), n=-1))
#渋谷区のポリゴン情報をテキストとして抜き出す。
df<-fetch(dbSendQuery(conn, "SELECT ST_AsText(geom) geom, jinko from h22ka13113"), n=-1)
plot(readWKT(text=df$geom[1], id=rownames(df)[1]))
#rgeos::readWKTを使ってDBから取得しcharacter型のWKTをRのspクラスに変換する。
spps<-SpatialPolygons(sapply(1:nrow(df), function(i){readWKT(df[i, "geom"], id=rownames(df)[i])@polygons}))
spdf<-SpatialPolygonsDataFrame(Sr=spps, data=df)
#1000人ずつで7色に塗り分ける。
plot(spdf,col=heat.colors(7)[as.integer(cut(spdf$jinko, breaks=0:7*1000))])
title("渋谷区人口")
legend(x=spps@bbox[1, 1]-0.00005, y=spps@bbox[2, 1]+0.02, title="人口",
legend=levels(cut(spdf$jinko, breaks=0:7*1000)), fill=heat.colors(7), bg="white", cex=0.8)
enjoy!
お役立ちリンク
| タイトル | URL | 摘要 |
| PostGIS 2.0.0マニュアル日本語訳 | http://www.finds.jp/docs/pgisman/2.0.0/postgis.html | マニュアルです |
| PostGIS入門 | http://workshops.boundlessgeo.com/postgis-intro-jp/index.html | PostGIS使い方に関する包括的なチュートリアル、一度は目を通すべき |
| GISのための測地成果、測地系、楕円体、投影座標系、EPSGコードのまとめ | http://d.hatena.ne.jp/tmizu23/20091215/1260868350 | 測地系、EPSGコードに関するよい解説 |
| PostGISとは？ | http://lets.postgresql.jp/documents/tutorial/PostGIS/1/ | PostGISに関する基礎的なチュートリアル、少し情報は古いが、インストール周りの解説が丁寧で読みやすい |
