1# ******************************************************************************
2#  $Id: lookup.py d5e0d780b2b0b96cf6933d358ea7927be776dfbb 2018-05-11 19:04:41 +1000 Ben Elliston $
3#
4#  Project:  GDAL ECW Driver
5#  Purpose:  Script to lookup ECW (GDT) coordinate systems and translate
6#            into OGC WKT for storage in $GDAL_HOME/data/ecw_cs.wkt.
7#  Author:   Frank Warmerdam, warmerdam@pobox.com
8#
9# ******************************************************************************
10#  Copyright (c) 2003, Frank Warmerdam
11#
12#  Permission is hereby granted, free of charge, to any person obtaining a
13#  copy of this software and associated documentation files (the "Software"),
14#  to deal in the Software without restriction, including without limitation
15#  the rights to use, copy, modify, merge, publish, distribute, sublicense,
16#  and/or sell copies of the Software, and to permit persons to whom the
17#  Software is furnished to do so, subject to the following conditions:
18#
19#  The above copyright notice and this permission notice shall be included
20#  in all copies or substantial portions of the Software.
21#
22#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23#  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25#  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27#  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28#  DEALINGS IN THE SOFTWARE.
29# ******************************************************************************
30
31import osr
32
33##############################################################################
34# rtod(): radians to degrees.
35
36
37def r2d(rad):
38    return float(rad) / 0.0174532925199433
39
40##############################################################################
41# load_dict()
42
43
44def load_dict(filename):
45    lines = open(filename).readlines()
46
47    this_dict = {}
48    for line in lines:
49        if line[:8] != 'proj_name':
50            tokens = [token.strip() for token in line.split(',')]
51
52            this_dict[tokens[0]] = tokens
53
54    return this_dict
55
56##############################################################################
57# Mainline
58
59
60# dir = 'M:/software/ER Viewer 2.0c/GDT_Data/'
61directory = '/u/data/ecw/gdt/'
62
63dict_list = ['tranmerc', 'lambert1', 'lamcon2', 'utm', 'albersea', 'mercator',
64             'obmerc_b', 'grinten', 'cassini', 'lambazea', 'datum_sp']
65
66dict_dict = {}
67
68for item in dict_list:
69    dict_dict[item] = load_dict(directory + item + '.dat')
70    # print 'loaded: %s' % item
71
72pfile = open(directory + 'project.dat')
73pfile.readline()
74
75for line in pfile.readlines():
76    try:
77        tokens = line.strip().split(',')
78        if len(tokens) < 3:
79            continue
80
81        tokens = [token.strip() for token in tokens]
82
83        ident = tokens[0]
84        typ = tokens[1]
85
86        lsize = float(tokens[2])
87        lsize_str = tokens[2]
88
89        dline = dict_dict[typ][ident]
90
91        srs = osr.SpatialReference()
92
93        # Handle translation of the projection parameters.
94
95        if type != 'utm':
96            fn = float(dline[1])
97            fe = float(dline[2])
98
99        if type == 'tranmerc':
100            srs.SetTM(r2d(dline[5]), r2d(dline[4]), float(dline[3]), fe, fn)
101
102        elif type == 'mercator':
103            srs.SetMercator(r2d(dline[5]), r2d(dline[4]), float(dline[3]), fe, fn)
104
105        elif type == 'grinten':
106            srs.SetVDG(r2d(dline[3]), fe, fn)
107
108        elif type == 'cassini':
109            srs.SetCS(r2d(dline[4]), r2d(dline[3]), fe, fn)
110
111        elif type == 'lambazea':
112            srs.SetLAEA(r2d(dline[5]), r2d(dline[4]),
113                        fe, fn)
114
115        elif type == 'lambert1':
116            srs.SetLCC1SP(r2d(dline[5]), r2d(dline[4]),
117                          float(dline[3]), fe, fn)
118
119        elif type == 'lamcon2':
120            srs.SetLCC(r2d(dline[7]), r2d(dline[8]),
121                       r2d(dline[9]), r2d(dline[6]), fe, fn)
122
123# elif type == 'lambert2':
124#   false_en = '+y_0=%.2f +x_0=%.2f' \
125#   % (float(dline[12])*lsize, float(dline[13])*lsize)
126#     result = '+proj=lcc %s +lat_0=%s +lon_0=%s +lat_1=%s +lat_2=%s' \
127#           % (false_en, r2d(dline[3]), r2d(dline[4]),
128#           r2d(dline[7]), r2d(dline[8]))
129
130        elif type == 'albersea':
131            srs.SetACEA(r2d(dline[3]), r2d(dline[4]),
132                        r2d(dline[5]), r2d(dline[6]), fe, fn)
133
134# elif type == 'obmerc_b':
135#     result = '+proj=omerc %s +lat_0=%s +lonc=%s +alpha=%s +k=%s' \
136#              % (false_en, r2d(dline[5]), r2d(dline[6]), r2d(dline[4]), dline[3])
137
138        elif type == 'utm':
139            srs.SetUTM(int(dline[1]), dline[2] != 'S')
140
141        # Handle Units from projects.dat file.
142        if srs.IsProjected():
143            srs.SetAttrValue('PROJCS', ident)
144            if lsize_str == '0.30480061':
145                srs.SetLinearUnits('US Foot', float(lsize_str))
146            elif lsize_str != '1.0':
147                srs.SetLinearUnits('unnamed', float(lsize_str))
148
149        wkt = srs.ExportToWkt()
150        if wkt:
151            print('%s,%s' % (ident, srs.ExportToWkt()))
152        else:
153            print('%s,LOCAL_CS["%s - (unsupported)"]' % (ident, ident))
154
155    except KeyError:
156        print('%s,LOCAL_CS["%s - (unsupported)"]' % (ident, ident))
157
158    except:
159        print('cant translate: ', line)
160        raise
161
162# Translate datums to their underlying spheroid information.
163
164pfile = open(directory + 'datum.dat')
165pfile.readline()
166
167for line in pfile.readlines():
168    tokens = [token.strip() for token in line.strip().split(',')]
169
170    ident = tokens[0]
171
172    sp_name = tokens[2]
173    dline = dict_dict['datum_sp'][ident]
174    srs = osr.SpatialReference()
175
176    if ident == 'WGS84':
177        srs.SetWellKnownGeogCS('WGS84')
178    elif ident == 'NAD27':
179        srs.SetWellKnownGeogCS('NAD27')
180    elif ident == 'NAD83':
181        srs.SetWellKnownGeogCS('NAD83')
182    else:
183        srs.SetGeogCS(tokens[1], ident, sp_name, float(dline[2]), float(dline[4]))
184
185    print('%s,%s' % (ident, srs.ExportToWkt()))
186