PostGisで緯度経度から地域情報を取得する

緯度経度から何県？とか知りたい。
PostGisでできるかやってみる。
基本的にココを参考にさせていただきました。
http://yamakk.com/blog/2010/07/21/import-polygon-data-into-postgis/

（手順）
　①国土数値情報ダウンロードサービスから元データをダウンロード
　②shapefileに変換
　③PostGisへインポート

①国土数値情報ダウンロードサービスから元データをダウンロード
国土数値情報ダウンロードサービス
http://nlftp.mlit.go.jp/ksj/jpgis/jpgis_datalist.html
->（JPGIS準拠）データのダウンロード
->行政区域（面）
に行って
「ダウンロードするデータの選択」で47都道府県をすべて選択する。
次の画面で「世界測地系」かつ最新のzipファイルをすべてダウンロードする。
zipファイルを解凍すると
N03-110331_01.xml
とかN03-....xmlというファイルがたくさんできる。

②shapefileに変換
ココから変換ツールをダウンロードする。(Windows版のみ)
http://nlftp.mlit.go.jp/ksj/jpgis/jpgis_tool.html
インストールして起動してN03-....xmlファイルを一気にshapefileに変換する。
N03-....shp、N03-....shx、N03-....dbfなど３種類のファイルができる。

③PostGisへインポート
スクリプトを作る。
N03-....shpファイルがある同じディレクトリにこのファイルをおく。

# -*- coding: utf-8 -*-
import sys
import logging
import os
import glob

# DB_HOST='***.***.***.***'
DB_HOST='localhost'

def start_trans(fp):
    try:
        logging.debug('start_trans start')

        for shapefile in glob.glob('N03*.shp'):
            create_sql = '%s_createtable.sql' % shapefile
            create_cmd = 'shp2pgsql -p %s japan > %s' % (shapefile, create_sql)
            insert_sql = '%s_insert.sql' % shapefile
            insert_cmd = 'shp2pgsql -W CP932 -a %s japan > %s' % (shapefile, insert_sql)
            os.system(insert_cmd)

            load_create_cmd = 'psql -h %s -U hoge areadb < %s\n' % (DB_HOST, create_sql)
            load_insert_cmd = 'psql -h %s -U hoge areadb < %s\n' % (DB_HOST, insert_sql)

            print load_create_cmd
            fp.writelines(load_create_cmd)
            print load_insert_cmd
            fp.writelines(load_insert_cmd)

        logging.debug('start_trans end')
        return
    except:
        logging.error('start_trans %s', sys.exc_info())
        raise            os.system(create_cmd)

if __name__ == "__main__":
    try:
#        logging.getLogger().setLevel(logging.ERROR)
#        logging.getLogger().setLevel(logging.DEBUG)
        logging.getLogger().setLevel(logging.INFO)
        filename = 'polygon.txt'
        fp = open(filename, 'w')
        start_trans(fp)
        fp.close()

    except HTTPError, e:
        print e.code
        print e.read()

これで、polygon.txtファイルができる.
実行する前にデータベースを作っておく。

createdb -E UTF8 template_postgis
$psql
# postgres=create role hoge superuser createdb login password '*****';
# postgres=\q
$createlang -d template_postgis plpgsql
$psql -d postgres -c "UPDATE pg_database SET datistemplate='true' WHERE datname='template_postgis';"
$psql -d template_postgis -f /usr/local/pgsql/share/contrib/postgis-1.5/postgis.sql 
$psql -d template_postgis -c "GRANT ALL ON geometry_columns TO PUBLIC;" 
$psql -d template_postgis -c "GRANT ALL ON spatial_ref_sys TO PUBLIC;"
$createdb -U username -E utf-8 -W -e areadb

それからさきほどのファイルを実行する。

source polygon.txt

インデックス登録。

areadb=#CREATE INDEX japan_gis_idx ON japan USING GIST (the_geom GIST_GEOMETRY_OPS);

これで終わり。

実際に検索してみる。

areadb=# SELECT prn, cn2 FROM japan WHERE ST_Contains(the_geom, GeomFromText('POINT(135.182734 34.690516)'));
  prn   |  cn2   
--------+--------
 兵庫県 | 中央区
(1 row)

うむ、できてる。

参考：