Import polygon data into PostGIS

国土数値情報統一フォーマットのダウンロードサービス 行政界・海岸線（面）のデータを
国土数値情報データ変換ツールをつかってShapeに変換したものをPostGISにinsertする。

参考 国土交通省にあるGISデータをPostGISへインポート(改訂版)

流れとしては,
国土数値情報データ–>(Ksjtool)–>Shapefile–(PostGIS shp2pgsql)–>SQL
こんな感じ。

47都道府県のPolygonをバッチでPostGIS DB入れられるようスクリプトを書いた.
import_polygon.py

#coding:utf-8
import os
import glob

for dir in glob.glob('N03*'):
    basename = os.path.basename(dir)
    _shp = '%s*.shp' % (os.path.join(dir,basename))
    shapefile = shapefile = glob.glob(_shp)[0]
    
    create_sql = '%s_createtable.sql' % os.path.join(dir,basename)
    create_cmd = 'shp2pgsql -p %s japan > %s' % (shapefile, create_sql)
    print create_cmd
    os.system(create_cmd)

    insert_sql = '%s_insert.sql' % os.path.join(dir,basename)
    insert_cmd = 'shp2pgsql -W cp932 -a %s japan > %s' % (shapefile, insert_sql)
    print insert_cmd
    os.system(insert_cmd)
    
    load_create_cmd = 'psql -U postgres gistest < %s' % create_sql
    load_insert_cmd = 'psql -U postgres gistest < %s' % insert_sql

    print load_create_cmd
    os.system(load_create_cmd)
    print load_insert_cmd
    os.system(load_insert_cmd)

実行

$ python import_polygon.py

SQLで入っているか確認

$ psql -U postgres gistest
psql (8.4.4)
Type "help" for help.

gistest=# SELECT prn, cn2, SUBSTRING(AsText(the_geom), 0, 100) from japan limit 5;
 兵庫県 | 豊岡市 | MULTIPOLYGON(((134.793994 35.668108,134.793978 35.668144,134.793915 35.668139,134.793854 35.668166,
 兵庫県 | 豊岡市 | MULTIPOLYGON(((134.762469 35.674204,134.762507 35.674247,134.762601 35.67427,134.762675 35.674276,1
 兵庫県 | 豊岡市 | MULTIPOLYGON(((134.762035 35.673419,134.762049 35.673454,134.762098 35.673476,134.762188 35.673457,
 兵庫県 | 豊岡市 | MULTIPOLYGON(((134.761883 35.672992,134.761889 35.673025,134.761921 35.67304,134.761951 35.673035,1
 兵庫県 | 香美町 | MULTIPOLYGON(((134.535152 35.66762,134.535201 35.66763,134.535221 35.667637,134.535186 35.667659,13

吉祥寺駅付近の(x,y)=(35.7031,139.5797)で検索 //PostGISの座標表記はすべてyx順

$ psql -U postgres gistest
psql (8.4.4)
Type "help" for help.

gistest=#SELECT prn, cn2 FROM japan WHERE ST_Contains(the_geom, GeomFromText('POINT(139.5731 35.7031)'));
  prn   |   cn2    
--------+----------
 東京都 | 武蔵野市
(1 row)

Spatial Indexを作成する

$ psql -U postgres gistest
psql (8.4.4)
Type "help" for help.

gistest=#CREATE INDEX inx_japan ON japan USING GIST (the_geom GIST_GEOMETRY_OPS);
CREATE INDEX

gistest=# VACUUM ANALYZE japan;
VACUUM

MapServerで確認
japan.map

MAP
    SIZE 800 800           #画像size
    EXTENT 128 33 150 44   #出力範囲の座標
    STATUS ON              #地図を表示するか
    UNITS DD               #地図の単位(DD は緯度経度)
    IMAGECOLOR 255 255 255 #背景色R G B
    IMAGETYPE PNG          #地図画像を保存する形式
    LAYER
         NAME "japan"
         CONNECTIONTYPE POSTGIS
         CONNECTION "user=postgres password='' dbname=gistest host=localhost"
         DATA "the_geom FROM japan" #select文
         TYPE LINE
         STATUS ON
         CLASS
             COLOR 0 0 0
         END
    END
END

$ sudo yum install gd-devel
$ wget http://download.osgeo.org/mapserver/mapserver-5.6.5.tar.gz
$ ./configure --with-postgis=/usr/local/pgsql/bin/pg_config --with-proj=/usr/local
$ make
$ ./shp2img -m japan.map -o japan.png

japan.png