1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import os.path as op
9from pyrocko import config, util
10from .srtmgl3 import SRTMGL3, AuthenticationRequired # noqa
11from .etopo1 import ETOPO1
12from . import dataset, tile
13from pyrocko.plot.cpt import get_cpt_path as cpt # noqa
15__doc__ = util.parse_md(__file__)
17positive_region = tile.positive_region
19earthradius = 6371000.0
20r2d = 180./math.pi
21d2r = 1./r2d
22km = 1000.
23d2m = d2r*earthradius
24m2d = 1./d2m
26topo_data_dir = config.config().topo_dir
28srtmgl3 = SRTMGL3(
29 name='SRTMGL3',
30 data_dir=op.join(topo_data_dir, 'SRTMGL3'))
32srtmgl3_d2 = dataset.DecimatedTiledGlobalDataset(
33 name='SRTMGL3_D2',
34 base=srtmgl3,
35 ndeci=2,
36 data_dir=op.join(topo_data_dir, 'SRTMGL3_D2'))
38srtmgl3_d4 = dataset.DecimatedTiledGlobalDataset(
39 name='SRTMGL3_D4',
40 base=srtmgl3_d2,
41 ndeci=2,
42 data_dir=op.join(topo_data_dir, 'SRTMGL3_D4'))
44srtmgl3_d8 = dataset.DecimatedTiledGlobalDataset(
45 name='SRTMGL3_D8',
46 base=srtmgl3_d4,
47 ndeci=2,
48 ntx=1001,
49 nty=1001,
50 data_dir=op.join(topo_data_dir, 'SRTMGL3_D8'))
52etopo1 = ETOPO1(
53 name='ETOPO1',
54 data_dir=op.join(topo_data_dir, 'ETOPO1'))
56etopo1_d2 = dataset.DecimatedTiledGlobalDataset(
57 name='ETOPO1_D2',
58 base=etopo1,
59 ndeci=2,
60 data_dir=op.join(topo_data_dir, 'ETOPO1_D2'))
62etopo1_d4 = dataset.DecimatedTiledGlobalDataset(
63 name='ETOPO1_D4',
64 base=etopo1_d2,
65 ndeci=2,
66 data_dir=op.join(topo_data_dir, 'ETOPO1_D4'))
68etopo1_d8 = dataset.DecimatedTiledGlobalDataset(
69 name='ETOPO1_D8',
70 base=etopo1_d4,
71 ndeci=2,
72 data_dir=op.join(topo_data_dir, 'ETOPO1_D8'))
74srtmgl3_all = [
75 srtmgl3,
76 srtmgl3_d2,
77 srtmgl3_d4,
78 srtmgl3_d8]
80etopo1_all = [
81 etopo1,
82 etopo1_d2,
83 etopo1_d4,
84 etopo1_d8]
86dems = srtmgl3_all + etopo1_all
89def make_all_missing_decimated():
90 for dem in dems:
91 if isinstance(dem, dataset.DecimatedTiledGlobalDataset):
92 dem.make_all_missing()
95def comparison(region, dems=dems):
96 import matplotlib.pyplot as plt
98 west, east, south, north = tile.positive_region(region)
100 fig = plt.gcf()
102 for idem, dem_ in enumerate(dems):
103 fig.add_subplot(len(dems), 1, idem+1)
104 t = dem_.get(region)
105 if t:
106 plt.pcolormesh(t.x(), t.y(), t.data)
107 plt.title(dem_.name)
108 plt.xlim(west, east)
109 plt.ylim(south, north)
111 plt.show()
114class UnknownDEM(Exception):
115 pass
118def dem_names():
119 return [dem.name for dem in dems]
122def dem(dem_name):
123 for dem in dems:
124 if dem.name == dem_name:
125 return dem
127 raise UnknownDEM(dem_name)
130def get(dem_name, region):
131 return dem(dem_name).get(region)
134def elevation(lat, lon):
135 '''
136 Get elevation at given point.
138 Tries to use SRTMGL3, falls back to ETOPO01 if not available.
139 '''
141 for dem in ['SRTMGL3', 'ETOPO1']:
142 r = get(dem, (lon, lat))
143 if r is not None and r != 0:
144 return r
147def select_dem_names(kind, dmin, dmax, region, mode='highest'):
148 assert kind in ('land', 'ocean')
149 assert mode in ('highest', 'lowest')
151 step = -1 if mode == 'lowest' else 1
153 ok = []
154 if kind == 'land':
155 for dem in srtmgl3_all[::step]:
156 if dem.is_suitable(region, dmin, dmax):
157 ok.append(dem.name)
158 break
160 for dem in etopo1_all[::step]:
161 if dem.is_suitable(region, dmin, dmax):
162 ok.append(dem.name)
163 break
165 return ok
168if __name__ == '__main__':
169 # comparison((-180., 180., -90, 90), dems=[etopo1_d8])
170 util.setup_logging('topo', 'info')
171 comparison((30, 31, 30, 31), dems=[srtmgl3, srtmgl3_d2])