source: anuga_work/development/cairns_2006/convert_lat_long.py @ 6823

"""
Author: John Jakeman
Created: 12/12/2005
Program used to convert a commer seperated .xya file consiting of
X and Y corrdinates in degrees and an elevation in meters into an .xya file
in which all quantities are measured in meters
"""

from Numeric import array, zeros, Float, asarray

# Open .xya file to be converted
filename  = 'cairns.xya'
fid = open(filename)

# Skip first line
line = fid.readline()

# Read remaining lines
lines = fid.readlines()
fid.close()

r_earth = 6.378135E+06
# Origin of cartesian system is the point (lat_origin, long_origin)
lat_origin = 145.0
long_origin = -24.0

from math import cos, pi

# fields[0]: latitude
# fields[1]: longitude
# fields[2]: elevation

# Open .xya file to be used to store new coordinates in meters
filename = 'cairns_test.xya'
fid = open(filename, "w")
fid.write('elevation\n')


for i, line in enumerate(lines):
    fields = line.split(',')
    
    # Evalulates distance in X-direction from origin for points in
    # Western Hemispehere 
    if float(fields[0]) < 0:
        X = (360.0+float(fields[0])-lat_origin)/360.0*abs(cos(float(fields[1])*pi/180.0))*r_earth*2.0*pi
    
    elif float(fields[0]) < lat_origin: # Only needed when lat_origin > 0
        X = (360.0-float(fields[0])-lat_origin)/360.0*abs(cos(float(fields[1])*pi/180.0))*r_earth*2.0*pi
        
    # Evaluates distance from origin in X-direction for points in
    # Eastern Hemisphere
    else:
        X = (float(fields[0])-lat_origin)/360.0*abs(cos(float(fields[1])*pi/180.0))*r_earth*2.0*pi
        
    # Calulates distance in Y-direction from a point to origin
    Y = abs((float(fields[1])-long_origin)/360.0*2*pi*r_earth)
    if Y < 0 :
        Y = 2*pi*r_earth + Y

    # Write X and Y corrdinates and corresponding elevation to .xya file
    fid.write(repr(X) + ', ')
    fid.write(repr(Y) + ', ')
    fid.write(fields[2])

fid.close()