PostGISにEPSGコードを追加する

##状況・背景

• PostGISで利用できる空間参照系は、spatial_ref_sysテーブルに登録されているもののみ
• マニアックな空間参照系（今回の例では、EPSG:55005）は定義されていない
• sparial_ref_sys登録には、投影情報のproj記法とwkt記法の両方が必要
• projまたはwktのどちらかだけが分かっている（今回の例では、proj記法だけが公開されている）

##環境

CentOS 7.1
PostgreSQL 9.4
PostGIS 2.94

##projからwktを求める

GDALのライブラリを利用する。下記参考サイトにproj/wkt双方向の変換が可能なpythonスクリプトがあったので拝借。
参考サイト：[Converting WKT projection info to Proj4 (and back)]
(http://spatialnotes.blogspot.jp/2010/11/converting-wkt-projection-info-to-proj4.html)

pythonがインストールされていることを確認し、

$ python -V
Python 2.7.5

yumでgdalをインストール。

$ sudo yum install -y gdal gdal-devel gdal-python

あとは、参考サイトに掲載されたスクリプト（↓↓↓）を適当なところに置き、

proj2wkt.py

#!/usr/bin/env python

import os
import sys
import string
import osgeo.osr

if (len(sys.argv) <> 2):
        print 'Usage: proj2wkt.py [Proj4 Projection Text]'
else:
        srs = osgeo.osr.SpatialReference()
        srs.ImportFromProj4(sys.argv[1])
        print srs.ExportToWkt()

パラメータにproj記法の文字列を指定して実行

$ ./proj2wkt.py "+proj=eqc +a=6377397.155 +b=6356079 +lat_ts=35.43094469 +lon_0=0 +lat_0=0 +k=1.0 +x_0=0 +y_0=0 +units=m no_defs"
PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6377397.155,299.1533345638937]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Equirectangular"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["standard_parallel_1",35.43094469],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

##PostGISに登録

上で求めたwktを、projと合わせてspatial_ref_sysにインポート。

INSERT INTO spatial_ref_sys(srid, auth_name, auth_srid, srtext, proj4text) VALUES (55005, 'EPSG', 55005, 'PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6377397.155,299.1533345638937]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Equirectangular"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",0],PARAMETER["standard_parallel_1",35.43094469],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]', '+proj=eqc +a=6377397.155 +b=6356079 +lat_ts=35.43094469 +lon_0=0 +lat_0=0 +k=1.0 +x_0=0 +y_0=0 +units=m no_defs');

##補足

今回のケースではprojからwktを求めて登録したが、逆のパターン（wktが分かっているけどprojが不明なケース）にも対応可能。その際は以下のスクリプトを利用する。

wkt2proj.py

#!/usr/bin/env python

import os
import sys
import string
import osgeo.osr

if (len(sys.argv) <> 2):
        print 'Usage: wkt2proj.py [WKT Projection Text]'
else:
        srs = osgeo.osr.SpatialReference()
        srs.ImportFromWkt(sys.argv[1])
        print srs.ExportToProj4()

以上