The conversion formula from latitude / longitude coordinates to UTM (Universal Transverse Mercator) coordinates is a bit cumbersome. A little complicated calculation is required to obtain the meridian arc length from the equator given an arbitrary latitude. With python, you can easily convert coordinates by using the pyproj module. The reverse is also possible. It cannot be used as it is, so it is necessary to make some modifications.
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]")
#OUTPUT
UTM zone is 54
UTM Easting is 388435.0160657919 [m]
UTM Northing is 3949276.5698600397 [m]
Verification Source: https://www.latlong.net/lat-long-utm.html
UTM → EQA To convert from UTM coordinates, it is necessary to specify the UTM zone and the northern and southern hemispheres.
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.]
It's almost the same.
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
The color is UTM Easting [m] grid
Recommended Posts