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

Last change on this file since 5587 was 2229, checked in by steve, 18 years ago

Moved directories into production and development parent directories

File size: 1.8 KB
Line 
1"""
2Author: John Jakeman
3Created: 12/12/2005
4Program used to convert a commer seperated .xya file consiting of
5X and Y corrdinates in degrees and an elevation in meters into an .xya file
6in which all quantities are measured in meters
7"""
8
9from Numeric import array, zeros, Float, asarray
10
11# Open .xya file to be converted
12filename  = 'cairns.xya'
13fid = open(filename)
14
15# Skip first line
16line = fid.readline()
17
18# Read remaining lines
19lines = fid.readlines()
20fid.close()
21
22r_earth = 6.378135E+06
23# Origin of cartesian system is the point (lat_origin, long_origin)
24lat_origin = 145.0
25long_origin = -24.0
26
27from math import cos, pi
28
29# fields[0]: latitude
30# fields[1]: longitude
31# fields[2]: elevation
32
33# Open .xya file to be used to store new coordinates in meters
34filename = 'cairns_test.xya'
35fid = open(filename, "w")
36fid.write('elevation\n')
37
38
39for i, line in enumerate(lines):
40    fields = line.split(',')
41   
42    # Evalulates distance in X-direction from origin for points in
43    # Western Hemispehere
44    if float(fields[0]) < 0:
45        X = (360.0+float(fields[0])-lat_origin)/360.0*abs(cos(float(fields[1])*pi/180.0))*r_earth*2.0*pi
46   
47    elif float(fields[0]) < lat_origin: # Only needed when lat_origin > 0
48        X = (360.0-float(fields[0])-lat_origin)/360.0*abs(cos(float(fields[1])*pi/180.0))*r_earth*2.0*pi
49       
50    # Evaluates distance from origin in X-direction for points in
51    # Eastern Hemisphere
52    else:
53        X = (float(fields[0])-lat_origin)/360.0*abs(cos(float(fields[1])*pi/180.0))*r_earth*2.0*pi
54       
55    # Calulates distance in Y-direction from a point to origin
56    Y = abs((float(fields[1])-long_origin)/360.0*2*pi*r_earth)
57    if Y < 0 :
58        Y = 2*pi*r_earth + Y
59
60    # Write X and Y corrdinates and corresponding elevation to .xya file
61    fid.write(repr(X) + ', ')
62    fid.write(repr(Y) + ', ')
63    fid.write(fields[2])
64
65fid.close()
Note: See TracBrowser for help on using the repository browser.