謎の座標系をEPSGに総当りして特定する

謎の直交座標系のファイルを目前に頭を抱えた事はないだろうか。

...普通はないと思う。
これはそんなDEMファイルを処理しないといけなくなった人の物語。
Qiita初投稿です。本業は海底火山ですがそれはまた別の機会に。

XとYがこんな感じの数字で与えられている海底地形のDEMがあった。
(6811000, -2740000)
CRSが分からないのでQGISに投影できない。

まず勘で見当をつける。

数字的に直交座標系であることは間違いなさそう。単位をメートルとすると、1000kmのオーダーになるのでいかにもそんな感じである。
しかしUTM座標系であれば多くの場合Xが6桁、Yが7桁の正の数字になるところが、この座標はXが7桁、Yがなんと負の7桁である。
UTMではないか、UTMだとしても間違ったゾーン設定で作られた可能性が高い。
場所はニュージーランド沖なのでNZ国土地理院にある直交座標系をすべて試したが、全部違った。

やっぱり自動化

そういう訳で特定は諦め、機械に頼ることにした。
方針は以下の3点。

1. Projモジュールの座標変換機能(transform)を使ってEPSGコード(取り敢えず1000~40000まで)に総当りを行う
2. 事前情報としてこの座標がNZ北方、南緯-32~-30度の間にある事が分かっていたので、WGS84に変換を試みた後、これを用いて結果をフィルタリング(def isvalid)する
3. XとYが入れ替わっている可能性も考え、1つのEPSGコードごとに入力側と出力側合わせて4通りを検討する(def check)

　GIS_WarpWarp.py

from pyproj import Proj, transform
x1, y1 = 6811000, -2740000
WGS84 = Proj(init='epsg:4326')

## 出力された座標が妥当かを返す。事前情報はここに。
def isvalid(x,y):
    if y<-32.0 or -30.0<y: return False
    return True

def check(epsg):
    if(epsg % 10000 == 0): print('=> ',epsg)
    ## 第1関門：未定義EPSGだった場合はここでRuntimeErrorが返る。
    try:
        if(epsg<10**4): proj1 = Proj(init=('epsg:%04d' % epsg))
        elif(epsg<10**5): proj1 = Proj(init=('epsg:%05d' % epsg))
        elif(epsg<10**6): proj1 = Proj(init=('epsg:%06d' % epsg))
    except RuntimeError: return None
    ## 第2関門：座標変換を試みる。変換不可能ならRuntimeErrorが返る。
    try: x2,y2 = transform(proj1, WGS84, x1,y1)
    except RuntimeError: pass
    else:
        if isvalid(x2,y2): print(epsg, x1,y1, x2,y2)
        if isvalid(y2,x2): print(epsg, x1,y1, y2,x2)
    ## 第3関門：XとYを入れ替えて座標変換。
    try: x2,y2 = transform(proj1, WGS84, y1,x1)
    except RuntimeError: pass
    else:
        if isvalid(x2,y2): print(epsg, y1,x1, x2,y2)
        if isvalid(y2,x2): print(epsg, y1,x1, y2,x2)

## という感じでリスト内包表記を使って1000~40000に総当り。
void = [check(i) for i in range(1000,40000)]

実行すると、epsg, x1,y1, x2,y2という形でブアーッと出力される。

3795 6811000 -2740000 -12.248595429748068 -31.56469346681908
3950 6811000 -2740000 24.482815965032906 -31.5999781135824
3994 6811000 -2740000 -179.0469481850484 -31.112571615808832
4455 6811000 -2740000 1.938224862502011 -30.09909426737733
4855 -2740000 6811000 46.07814966253316 -31.334547234623837

さがしていたのはこれ(EPSG:3994 Mercator41)でした。なんじゃそりゃー。
どうもニュージーランド大気海洋局のCRSらしい。NZ国土地理院のサイトを探しても見つからないわけだ…。