pythonで緯度経度座標↔UTM座標変換 #Python

緯度経度座標からUTM (ユニバーサル横メルカトル) 座標への変換式はちょっと面倒．
任意の緯度が与えられた場合の赤道から子午線弧長を求めるためにちょっと面倒な計算が必要となる．
pythonならpyprojのモジュールを用いることで簡単に座標変換ができる．逆も可．
そのままでは使えないので少し手を加えてやる必要がある．

EQA → UTM

eqa2utm.py

from pyproj import Proj

#An arbitral point in EQA coordinate
lon=139.767120  #[deg.]
lat=35.681079  #[deg.]

##Compute UTM zone
#"int": Extract only integer value
#31: Offset for UTM zone definition
#6: Angle in a UTM zone for the longitude direction
e2u_zone=int(divmod(lon, 6)[0])+31

#Define EQA2UTM converter
e2u_conv=Proj(proj='utm', zone=e2u_zone, ellps='WGS84')
#Apply the converter
utmx, utmy=e2u_conv(lon, lat)
#Add offset if the point in the southern hemisphere
if lat<0:
    utmy=utmy+10000000

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

$python3 eqa2utm.py
#OUTPUT
 UTM zone is  54  
 UTM Easting is 388435.0160657919 [m]
 UTM Northing is  3949276.5698600397 [m]

確認

出典: https://www.latlong.net/lat-long-utm.html

UTM → EQA

UTM座標からの変換にはUTM zoneと南北半球の指定が必要．

utm2eqa.py

from pyproj import Proj

#Define UTM zone
e2u_zone=54
#Define the northern('N')/southern('S') hemisphere
hemi='N'
utmx=388435.0160657919 #[m]
utmy=4060202.406410741 #[m]

#Add offset if the point in the southern hemisphere
if hemi=='S':
    utmy=utmy-10000000

#Define coordinate converter
e2u_conv=Proj(proj='utm', zone=e2u_zone, ellps='WGS84')
#Convert UTM2EQA
lon, lat=e2u_conv(utmx, utmy, inverse=True)

print("Longitude is ",lon," [deg.] \n",\
      "Latitude is ", lat, "[deg.]")

#OUTPUT
Longitude is  139.75135058449402 [deg.]
 Latitude is  36.68091466691 [deg.]

Mesh丸ごと変換

EQA2UTM

おおむね同じです．

EQA2UTM_mesh.py

from pyproj import Proj
import numpy as np
import matplotlib.pyplot as plt

#np.linspace(start, end, number_of_data)
lon_1d=np.linspace(130.3, 131.05, 2701)
lat_1d=np.linspace(31.18333, 31.9333, 2701)
lon, lat=np.meshgrid(lon_1d, lat_1d)

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

#Option:
#Display the UTM Easting grid
plt.pcolormesh(utmx/1000, utmy/1000, utmx, cmap='jet')
plt.colorbar() 
plt.xlabel('UTM Easting [km]')
plt.ylabel('UTM Northing [km]')
plt.show

色はUTM Easting [m] のgrid