Coordonnées de latitude et de longitude ↔ Conversion de coordonnées UTM avec python

La formule de conversion des coordonnées latitude / longitude en coordonnées UTM (Universal Horizontal Mercator) est un peu gênante. Un calcul un peu compliqué est nécessaire pour trouver la longueur de l'arc méridien à partir de la ligne équatoriale étant donné une latitude arbitraire. Avec python, vous pouvez facilement convertir des coordonnées en utilisant le module pyproj. L'inverse est également possible. Il ne peut pas être utilisé tel quel, il doit donc être légèrement modifié.

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]

Vérification Source: https://www.latlong.net/lat-long-utm.html

UTM → EQA Pour convertir à partir des coordonnées UTM, il est nécessaire de spécifier la zone UTM et l'hémisphère nord-sud.

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.]

Conversion entière de maillage

• EQA2UTM

C'est presque pareil.

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 

La couleur est la grille UTM Est [m]