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