[PYTHON] Make a list of latitude and longitude and convert UTM coordinates at once → File output

From EQA list to batch conversion to UTM list creation

Previously I wrote python code to convert from EQA (latitude and longitude) to UTM (Universal Transverse Mercator) coordinates using python. ↓ Latitude / longitude coordinates ↔ UTM coordinate conversion in Python https://qiita.com/spicy_HotChicken/items/4cdf303493d73e24dc14

When converting the coordinates one by one, I entered the coordinate values in the python script. However, the number of conversions has gradually increased and it has become troublesome, so I thought it would be easier if I could create a list of latitude and longitude, convert them all at once, and output them as a text file. The latitude / longitude list is a coordinate conversion script when the latitude and longitude are separated in order from the left.

$ cat input.dat
139.2341321 35.534243
132.3214315 32.542664
138.4312145 36.534261

The conversion script is as follows. The input / output file name is specified at the beginning.

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)

In my environment, when text data is read by python, it becomes a string type list, so it is divided by * .read.split ('\ n') for each line break to make a list. This time, I only knew how to read line by line, so I read line by line in a loop and repeat conversion → file output. ʻArdata` makes the latitude / longitude coordinates a float type nparray. After that, the latitude and longitude are defined respectively, and the coordinates of EQA2UTM are converted.

Application example

I will try it.

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

339907.669515 3933725.28372 54
248457.470634 3603752.79472 53
270030.490036 4046278.55358 54

The input files are listed in the order of latitude and longitude from the left, and the output is listed in the order of UTM Easting, UTM Northing, UTM Zone from the left.

(Bonus) Interactively perform EQA2UTM conversion

In the EQA2UTM conversion introduced earlier, the latitude and longitude coordinates you want to coordinate were directly input to the script. This is troublesome. This time, I made a latitude / longitude list and converted it all at once. All you have to do is create a latitude / longitude list for which you want to convert coordinates. However, if you want to convert only one point but it is not enough to make a list, it is a calculation to enter coordinates interactively and output UTM coordinate values on the terminal screen. Is there such a person? .. ..

eqa2utm_int.py

from pyproj import Proj

#I changed only the bottom two lines ↓
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]

When you run the python script, you will be asked for latitude and longitude, so enter it and press Enter to display the UTM coordinate values on the terminal screen. It's easy.