Pyraid + GeoAlchemy2 + GeoFormAlchemy2 (3)
（承前）
Python コード内で GeoAlchemy2 を使う
サンプルとして下図のようなデータを作成した。
PostGIS 上で見ると、
gistest=> select * from plots;
id | name | geom
----+--------+----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
1 | points | 010700002004120000070000000101000000F2FFFFFF23932140DDFFFFFFF9364A400101000000E4FFFFFF47941E40DDFFFFFFF90C47400101000000E4FFFFFF47D81640DDFFFFFFF96647400101000000BAFFFFFF23C52B40DDFFFFFF797D47400101000000DDFFFFFF916B3B40DDFFFFFFB9B547400101000000DDFFFFFF91DC3340DDFFFFFF393449400101000000DDFFFFFF91553340DDFFFFFF79614540
(1 行)
一方、Pyramid 上でビュー函数を定義し、そこにデバッグ用 print 文を仕込んでみる。
@view_config(route_name='item', renderer='templates/plot.pt')
def item_view(request):
item_id = request.matchdict['id']
item = DBSession.query(Plot).filter(Plot.id == item_id).first()
print item.geom
すると、pserve を実行したコンソール上に同じ文字列が出力される。
0107000000070000000101000000f2ffffff23932140ddfffffff9364a400101000000e4ffffff47941e40ddfffffff90c47400101000000e4ffffff47d81640ddfffffff96647400101000000baffffff23c52b40ddffffff797d47400101000000ddffffff916b3b40ddffffffb9b547400101000000ddffffff91dc3340ddffffff393449400101000000ddffffff91553340ddffffff79614540
これは OpenGIS の WKB (Well-Known Binary) フォーマットのデータを hex 表現したもので、PostGIS に格納されている元々のデータは バイナリである（と思われる）。 item.geom オブジェクトは geoalchemy2.elements.WKBElement クラスの インスタンスで、このクラスは data と desc の属性を持つ。 data 属性はバイナリ表現のデータを保持し、desc は data の内容を hex 化して 返す、@property でデコレーションされたメソッドである。
WKB を WKT (Well-Known Text) 形式に変換するには主にふたつの方法がある。 PostGIS サーバに変換させる方法と、ローカルマシンで処理する方法である。
PostGIS サーバに変換させるには、geom カラムに ST_AsText() 函数を適用すれば良い。 ビュー函数のデバッグ文を以下のように変更する。
print DBSession.scalar(item.geom.ST_AsText())
すると出力は、
GEOMETRYCOLLECTION(POINT(8.7873840332031 52.429504394531),POINT(7.6448059082031 46.101379394531),POINT(5.7112121582031 46.804504394531),POINT(13.885040283203 46.980285644531),POINT(27.420196533203 47.419738769531),POINT(19.861602783203 50.408020019531),POINT(19.334259033203 42.761535644531))
コレクションの中身を分割させるには ST_Dump() 函数を使う。
a = DBSession.execute(item.geom.ST_Dump().geom.ST_AsText())
for i in a.fetchall():
print i
出力は、
(u'POINT(8.7873840332031 52.429504394531)',)
(u'POINT(7.6448059082031 46.101379394531)',)
(u'POINT(5.7112121582031 46.804504394531)',)
(u'POINT(13.885040283203 46.980285644531)',)
(u'POINT(27.420196533203 47.419738769531)',)
(u'POINT(19.861602783203 50.408020019531)',)
(u'POINT(19.334259033203 42.761535644531)',)
各行がタプルになっている点に注意が必要だ。
ローカルマシンに変換させる場合は geoalchemy2.shape.to_shape() を使う。 ビュー函数を修正してアプリケーション再起動。
from geoalchemy2.shape import to_shape
@view_config(route_name='item', renderer='templates/plot.pt')
def item_view(request):
item_id = request.matchdict['id']
item = DBSession.query(Plot).filter(Plot.id == item_id).first()
print to_shape(item.geom)
コンソール出力は以下のようになる。
GEOMETRYCOLLECTION (POINT (8.7873840332031 52.429504394531), POINT (7.6448059082031 46.101379394531), POINT (5.7112121582031 46.804504394531), POINT (13.885040283203 46.980285644531), POINT (27.420196533203 47.419738769531), POINT (19.861602783203 50.408020019531), POINT (19.334259033203 42.761535644531))
geoalchemy2.shape.to_shape() は shapely というパッケージを利用しており、 shapely は OS の libgeos を呼び出しているようだ。
to_shape() の返り値はオブジェクトで、この場合は <class ‘shapely.geometry.collection.GeometryCollection’> になる。このクラスはリストとして動作するのでスライス操作で要素を取り出すことが できる。 個々の要素は、この場合は <class ‘shapely.geometry.point.Point’> で、 POINT データの内容を更に操作することができる。たとえば X座標値だけを取り出すには
for i in shape.to_shape(item.geom):
print i.x
などとすれば良い。
Python コードとの相性は geoalchemy2.shape の方が良さそうだ。 ただ、PostGIS と libgeos がトポロジー計算などでまったく同じ精度の 計算を行うのかどうかは確かめないといけない。