1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, division
7import os
8import re
9import math
10import os.path as op
12from pyrocko import config, util
13from .srtmgl3 import SRTMGL3, AuthenticationRequired # noqa
14from .etopo1 import ETOPO1
15from . import dataset, tile
17__doc__ = util.parse_md(__file__)
19positive_region = tile.positive_region
21earthradius = 6371000.0
22r2d = 180./math.pi
23d2r = 1./r2d
24km = 1000.
25d2m = d2r*earthradius
26m2d = 1./d2m
28topo_data_dir = config.config().topo_dir
30srtmgl3 = SRTMGL3(
31 name='SRTMGL3',
32 data_dir=op.join(topo_data_dir, 'SRTMGL3'))
34srtmgl3_d2 = dataset.DecimatedTiledGlobalDataset(
35 name='SRTMGL3_D2',
36 base=srtmgl3,
37 ndeci=2,
38 data_dir=op.join(topo_data_dir, 'SRTMGL3_D2'))
40srtmgl3_d4 = dataset.DecimatedTiledGlobalDataset(
41 name='SRTMGL3_D4',
42 base=srtmgl3_d2,
43 ndeci=2,
44 data_dir=op.join(topo_data_dir, 'SRTMGL3_D4'))
46srtmgl3_d8 = dataset.DecimatedTiledGlobalDataset(
47 name='SRTMGL3_D8',
48 base=srtmgl3_d4,
49 ndeci=2,
50 ntx=1001,
51 nty=1001,
52 data_dir=op.join(topo_data_dir, 'SRTMGL3_D8'))
54etopo1 = ETOPO1(
55 name='ETOPO1',
56 data_dir=op.join(topo_data_dir, 'ETOPO1'))
58etopo1_d2 = dataset.DecimatedTiledGlobalDataset(
59 name='ETOPO1_D2',
60 base=etopo1,
61 ndeci=2,
62 data_dir=op.join(topo_data_dir, 'ETOPO1_D2'))
64etopo1_d4 = dataset.DecimatedTiledGlobalDataset(
65 name='ETOPO1_D4',
66 base=etopo1_d2,
67 ndeci=2,
68 data_dir=op.join(topo_data_dir, 'ETOPO1_D4'))
70etopo1_d8 = dataset.DecimatedTiledGlobalDataset(
71 name='ETOPO1_D8',
72 base=etopo1_d4,
73 ndeci=2,
74 data_dir=op.join(topo_data_dir, 'ETOPO1_D8'))
76srtmgl3_all = [
77 srtmgl3,
78 srtmgl3_d2,
79 srtmgl3_d4,
80 srtmgl3_d8]
82etopo1_all = [
83 etopo1,
84 etopo1_d2,
85 etopo1_d4,
86 etopo1_d8]
88dems = srtmgl3_all + etopo1_all
91def make_all_missing_decimated():
92 for dem in dems:
93 if isinstance(dem, dataset.DecimatedTiledGlobalDataset):
94 dem.make_all_missing()
97def cpt(name):
98 if os.path.exists(name):
99 return name
101 if not re.match(r'[A-Za-z0-9_]+', name):
102 raise Exception('invalid cpt name')
104 fn = util.data_file(os.path.join('colortables', '%s.cpt' % name))
105 if not os.path.exists(fn):
106 raise Exception('cpt file does not exist: %s' % fn)
108 return fn
111def comparison(region, dems=dems):
112 import matplotlib.pyplot as plt
114 west, east, south, north = tile.positive_region(region)
116 fig = plt.gcf()
118 for idem, dem_ in enumerate(dems):
119 fig.add_subplot(len(dems), 1, idem+1)
120 t = dem_.get(region)
121 if t:
122 plt.pcolormesh(t.x(), t.y(), t.data)
123 plt.title(dem_.name)
124 plt.xlim(west, east)
125 plt.ylim(south, north)
127 plt.show()
130class UnknownDEM(Exception):
131 pass
134def dem_names():
135 return [dem.name for dem in dems]
138def dem(dem_name):
139 for dem in dems:
140 if dem.name == dem_name:
141 return dem
143 raise UnknownDEM(dem_name)
146def get(dem_name, region):
147 return dem(dem_name).get(region)
150def elevation(lat, lon):
151 '''
152 Get elevation at given point.
154 Tries to use SRTMGL3, falls back to ETOPO01 if not available.
155 '''
157 for dem in ['SRTMGL3', 'ETOPO1']:
158 r = get(dem, (lon, lat))
159 if r is not None and r != 0:
160 return r
163def select_dem_names(kind, dmin, dmax, region):
164 assert kind in ('land', 'ocean')
165 ok = []
166 if kind == 'land':
167 for dem in srtmgl3_all:
168 if dem.is_suitable(region, dmin, dmax):
169 ok.append(dem.name)
170 break
172 for dem in etopo1_all:
173 if dem.is_suitable(region, dmin, dmax):
174 ok.append(dem.name)
175 break
177 return ok
180if __name__ == '__main__':
181 # comparison((-180., 180., -90, 90), dems=[etopo1_d8])
182 util.setup_logging('topo', 'info')
183 comparison((30, 31, 30, 31), dems=[srtmgl3, srtmgl3_d2])