1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import os
7import re
8import math
9import os.path as op
11from pyrocko import config, util
12from .srtmgl3 import SRTMGL3, AuthenticationRequired # noqa
13from .etopo1 import ETOPO1
14from . import dataset, tile
16__doc__ = util.parse_md(__file__)
18positive_region = tile.positive_region
20earthradius = 6371000.0
21r2d = 180./math.pi
22d2r = 1./r2d
23km = 1000.
24d2m = d2r*earthradius
25m2d = 1./d2m
27topo_data_dir = config.config().topo_dir
29srtmgl3 = SRTMGL3(
30 name='SRTMGL3',
31 data_dir=op.join(topo_data_dir, 'SRTMGL3'))
33srtmgl3_d2 = dataset.DecimatedTiledGlobalDataset(
34 name='SRTMGL3_D2',
35 base=srtmgl3,
36 ndeci=2,
37 data_dir=op.join(topo_data_dir, 'SRTMGL3_D2'))
39srtmgl3_d4 = dataset.DecimatedTiledGlobalDataset(
40 name='SRTMGL3_D4',
41 base=srtmgl3_d2,
42 ndeci=2,
43 data_dir=op.join(topo_data_dir, 'SRTMGL3_D4'))
45srtmgl3_d8 = dataset.DecimatedTiledGlobalDataset(
46 name='SRTMGL3_D8',
47 base=srtmgl3_d4,
48 ndeci=2,
49 ntx=1001,
50 nty=1001,
51 data_dir=op.join(topo_data_dir, 'SRTMGL3_D8'))
53etopo1 = ETOPO1(
54 name='ETOPO1',
55 data_dir=op.join(topo_data_dir, 'ETOPO1'))
57etopo1_d2 = dataset.DecimatedTiledGlobalDataset(
58 name='ETOPO1_D2',
59 base=etopo1,
60 ndeci=2,
61 data_dir=op.join(topo_data_dir, 'ETOPO1_D2'))
63etopo1_d4 = dataset.DecimatedTiledGlobalDataset(
64 name='ETOPO1_D4',
65 base=etopo1_d2,
66 ndeci=2,
67 data_dir=op.join(topo_data_dir, 'ETOPO1_D4'))
69etopo1_d8 = dataset.DecimatedTiledGlobalDataset(
70 name='ETOPO1_D8',
71 base=etopo1_d4,
72 ndeci=2,
73 data_dir=op.join(topo_data_dir, 'ETOPO1_D8'))
75srtmgl3_all = [
76 srtmgl3,
77 srtmgl3_d2,
78 srtmgl3_d4,
79 srtmgl3_d8]
81etopo1_all = [
82 etopo1,
83 etopo1_d2,
84 etopo1_d4,
85 etopo1_d8]
87dems = srtmgl3_all + etopo1_all
90def make_all_missing_decimated():
91 for dem in dems:
92 if isinstance(dem, dataset.DecimatedTiledGlobalDataset):
93 dem.make_all_missing()
96def cpt(name):
97 if os.path.exists(name):
98 return name
100 if not re.match(r'[A-Za-z0-9_]+', name):
101 raise Exception('invalid cpt name')
103 fn = util.data_file(os.path.join('colortables', '%s.cpt' % name))
104 if not os.path.exists(fn):
105 raise Exception('cpt file does not exist: %s' % fn)
107 return fn
110def comparison(region, dems=dems):
111 import matplotlib.pyplot as plt
113 west, east, south, north = tile.positive_region(region)
115 fig = plt.gcf()
117 for idem, dem_ in enumerate(dems):
118 fig.add_subplot(len(dems), 1, idem+1)
119 t = dem_.get(region)
120 if t:
121 plt.pcolormesh(t.x(), t.y(), t.data)
122 plt.title(dem_.name)
123 plt.xlim(west, east)
124 plt.ylim(south, north)
126 plt.show()
129class UnknownDEM(Exception):
130 pass
133def dem_names():
134 return [dem.name for dem in dems]
137def dem(dem_name):
138 for dem in dems:
139 if dem.name == dem_name:
140 return dem
142 raise UnknownDEM(dem_name)
145def get(dem_name, region):
146 return dem(dem_name).get(region)
149def elevation(lat, lon):
150 '''
151 Get elevation at given point.
153 Tries to use SRTMGL3, falls back to ETOPO01 if not available.
154 '''
156 for dem in ['SRTMGL3', 'ETOPO1']:
157 r = get(dem, (lon, lat))
158 if r is not None and r != 0:
159 return r
162def select_dem_names(kind, dmin, dmax, region):
163 assert kind in ('land', 'ocean')
164 ok = []
165 if kind == 'land':
166 for dem in srtmgl3_all:
167 if dem.is_suitable(region, dmin, dmax):
168 ok.append(dem.name)
169 break
171 for dem in etopo1_all:
172 if dem.is_suitable(region, dmin, dmax):
173 ok.append(dem.name)
174 break
176 return ok
179if __name__ == '__main__':
180 # comparison((-180., 180., -90, 90), dems=[etopo1_d8])
181 util.setup_logging('topo', 'info')
182 comparison((30, 31, 30, 31), dems=[srtmgl3, srtmgl3_d2])