緯度経度のリストを作って一括でUTM座標変換→ファイル出力

EQAリストから一括変換してUTMリスト作成まで

以前pythonを用いてEQA (緯度経度) からUTM (ユニバーサルトランスバースメルカトル) 座標へ変換するpythonコードを書きました．
↓Pythonで緯度経度座標↔UTM座標変換
https://qiita.com/spicy_HotChicken/items/4cdf303493d73e24dc14

ひとつひとつ座標変換するときにはpythonスクリプトに座標値を入力してしました．
しかし，だんだん変換する数も増えてきて面倒になってきたので，緯度経度のリストを作って一括変換してテキストファイルで出力できれば楽なのではないかと思いました．
緯度経度リストは左から緯度，経度の順でスペース区切りにしたときの座標変換スクリプトになっています．

$ cat input.dat
139.2341321 35.534243
132.3214315 32.542664
138.4312145 36.534261

変換スクリプトは以下の通りになります．
冒頭で入力・出力ファイル名を指定しています．

eqa2utm_list.py

import numpy as np
import os 
from pyproj import Proj

inpf='input.dat' #Define input file name.
outf='output.dat' #Define output file name.
#Check existence of output file.
#If output file exists, remove it.
if os.path.isfile(outf):
    os.remove(outf)

tmpdata=open(inpf, "r")
data=tmpdata.read().split('\n')
for a in range(len(data)-1):
    ardata=np.array(data[a].split(), float)
    #datastr=data[a]
    lon=ardata[0]
    lat=ardata[1]
    #print lon, lat

    #eqa2utm.py
    e2u_zone=int(divmod(lon, 6)[0])+31
    e2u_conv=Proj(proj='utm', zone=e2u_zone, ellps='WGS84')
    utmx, utmy=e2u_conv(lon, lat)
    if lat<0:
        utmy=utmy+10000000

    with open(outf, 'a') as f:
        outp=str(utmx)+str(' ')+str(utmy)+str(' ')+str(e2u_zone)
        f.write("%s\n" % outp)

print '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n',\
 ' EQA2UTM transformation Completed. \n',\
 '>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n'

#Show output file on a screen
f=open(outf, 'r')
data=f.read()
print(data)

私の環境だとpythonでテキストデータを読み込むとstring型のlistになるので，*.read.split('\n')で改行ごとに分割してlistにします．
今回は1行ずつ読んでいく方法しかわからなかったので，ループで1行ずつ読んで変換→ファイル出力を繰り返しています．
ardataは緯度経度座標をfloat型のnparrayにしています．
あとは緯度と経度をそれぞれ定義して，EQA2UTMの座標変換をしています．

適用例

試しに動かしてみます．

$python eqa2utm_list.py
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 EQA2UTM transformation Completed. 
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

339907.669515 3933725.28372 54
248457.470634 3603752.79472 53
270030.490036 4046278.55358 54

inputファイルは左から緯度，経度の順に，outputは左からUTM Easting, UTM Northing, UTM Zoneの順にリストになっています．

(おまけ) 対話的にEQA2UTM変換をする

先ほど紹介したEQA2UTMの変換はスクリプトに座標したい緯度経度座標を直接入力していました．これは面倒．
今回は緯度経度リストをつくって一括変換しました．座標変換したい緯度経度リストさえ作ってしまえばいいのです．
しかし1点だけ変換したいけどリストを作るほどでなければ，対話的に座標を入力してターミナル画面上にUTM座標値を出力するようにしてしまおうという算段です．そんな人いないか．．．

eqa2utm_int.py

from pyproj import Proj

#変えたのは下の2行だけです↓
lon=float(input('>> Longitude [deg.] >>')) 
lat=float(input('>> Latitude [deg.] >>'))  

e2u_zone=int(divmod(lon, 6)[0])+31
e2u_conv=Proj(proj='utm', zone=e2u_zone, ellps='WGS84')
utmx, utmy=e2u_conv(lon, lat)

if lat<0:
    utmy=utmy+10000000

print "UTM zone is ", e2u_zone, " \n", \
      "UTM Easting is", utmx, "[m]\n",\
      "UTM Northing is ", utmy, "[m]"

#EOF

$ python eqa2utm_int.py
>> Longitude [deg.] >>132.55
>> Latitude [deg.] >>38.222
UTM zone is  53  
UTM Easting is 285531.137964 [m]
UTM Northing is  4233284.85815 [m]

pythonスクリプトを動かすと緯度と経度を聞かれるので，入力してEnter押下でUTM座標値をターミナル画面上に表示します．楽ちんです．

-