1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6import math 

7import os.path as op 

8 

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 

14 

15__doc__ = util.parse_md(__file__) 

16 

17positive_region = tile.positive_region 

18 

19earthradius = 6371000.0 

20r2d = 180./math.pi 

21d2r = 1./r2d 

22km = 1000. 

23d2m = d2r*earthradius 

24m2d = 1./d2m 

25 

26topo_data_dir = config.config().topo_dir 

27 

28srtmgl3 = SRTMGL3( 

29 name='SRTMGL3', 

30 data_dir=op.join(topo_data_dir, 'SRTMGL3')) 

31 

32srtmgl3_d2 = dataset.DecimatedTiledGlobalDataset( 

33 name='SRTMGL3_D2', 

34 base=srtmgl3, 

35 ndeci=2, 

36 data_dir=op.join(topo_data_dir, 'SRTMGL3_D2')) 

37 

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')) 

43 

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')) 

51 

52etopo1 = ETOPO1( 

53 name='ETOPO1', 

54 data_dir=op.join(topo_data_dir, 'ETOPO1')) 

55 

56etopo1_d2 = dataset.DecimatedTiledGlobalDataset( 

57 name='ETOPO1_D2', 

58 base=etopo1, 

59 ndeci=2, 

60 data_dir=op.join(topo_data_dir, 'ETOPO1_D2')) 

61 

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')) 

67 

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')) 

73 

74srtmgl3_all = [ 

75 srtmgl3, 

76 srtmgl3_d2, 

77 srtmgl3_d4, 

78 srtmgl3_d8] 

79 

80etopo1_all = [ 

81 etopo1, 

82 etopo1_d2, 

83 etopo1_d4, 

84 etopo1_d8] 

85 

86dems = srtmgl3_all + etopo1_all 

87 

88 

89def make_all_missing_decimated(): 

90 for dem in dems: 

91 if isinstance(dem, dataset.DecimatedTiledGlobalDataset): 

92 dem.make_all_missing() 

93 

94 

95def comparison(region, dems=dems): 

96 import matplotlib.pyplot as plt 

97 

98 west, east, south, north = tile.positive_region(region) 

99 

100 fig = plt.gcf() 

101 

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) 

110 

111 plt.show() 

112 

113 

114class UnknownDEM(Exception): 

115 pass 

116 

117 

118def dem_names(): 

119 return [dem.name for dem in dems] 

120 

121 

122def dem(dem_name): 

123 for dem in dems: 

124 if dem.name == dem_name: 

125 return dem 

126 

127 raise UnknownDEM(dem_name) 

128 

129 

130def get(dem_name, region): 

131 return dem(dem_name).get(region) 

132 

133 

134def elevation(lat, lon): 

135 ''' 

136 Get elevation at given point. 

137 

138 Tries to use SRTMGL3, falls back to ETOPO01 if not available. 

139 ''' 

140 

141 for dem in ['SRTMGL3', 'ETOPO1']: 

142 r = get(dem, (lon, lat)) 

143 if r is not None and r != 0: 

144 return r 

145 

146 

147def select_dem_names(kind, dmin, dmax, region, mode='highest'): 

148 assert kind in ('land', 'ocean') 

149 assert mode in ('highest', 'lowest') 

150 

151 step = -1 if mode == 'lowest' else 1 

152 

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 

159 

160 for dem in etopo1_all[::step]: 

161 if dem.is_suitable(region, dmin, dmax): 

162 ok.append(dem.name) 

163 break 

164 

165 return ok 

166 

167 

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])