Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/dataset/topo/__init__.py: 73%

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6''' 

7Multiresolution topography and bathymetry datasets. 

8 

9This module provides access to gridded topography data from the following 

10sources: 

11 

12* NOAA `ETOPO1 <https://www.ngdc.noaa.gov/mgg/global/>`_ Global Relief Map. 

13* NASA Shuttle Radar Topography Mission Global 30 arc second (`SRTMGL3 

14 <https://lpdaac.usgs.gov/dataset_discovery/measures/measures_products_table/srtmgl3_v003>`_) 

15 topography data. 

16 

17''' 

18 

19import math 

20import os.path as op 

21 

22from pyrocko import config, util 

23from .srtmgl3 import SRTMGL3, AuthenticationRequired # noqa 

24from .etopo1 import ETOPO1 

25from . import dataset, tile 

26from pyrocko.plot.cpt import get_cpt_path as cpt # noqa 

27 

28positive_region = tile.positive_region 

29 

30earthradius = 6371000.0 

31r2d = 180./math.pi 

32d2r = 1./r2d 

33km = 1000. 

34d2m = d2r*earthradius 

35m2d = 1./d2m 

36 

37topo_data_dir = config.config().topo_dir 

38 

39_srtmgl3 = SRTMGL3( 

40 name='SRTMGL3', 

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

42 

43_srtmgl3_d2 = dataset.DecimatedTiledGlobalDataset( 

44 name='SRTMGL3_D2', 

45 base=_srtmgl3, 

46 ndeci=2, 

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

48 

49_srtmgl3_d4 = dataset.DecimatedTiledGlobalDataset( 

50 name='SRTMGL3_D4', 

51 base=_srtmgl3_d2, 

52 ndeci=2, 

53 data_dir=op.join(topo_data_dir, 'SRTMGL3_D4')) 

54 

55_srtmgl3_d8 = dataset.DecimatedTiledGlobalDataset( 

56 name='SRTMGL3_D8', 

57 base=_srtmgl3_d4, 

58 ndeci=2, 

59 ntx=1001, 

60 nty=1001, 

61 data_dir=op.join(topo_data_dir, 'SRTMGL3_D8')) 

62 

63_etopo1 = ETOPO1( 

64 name='ETOPO1', 

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

66 

67_etopo1_d2 = dataset.DecimatedTiledGlobalDataset( 

68 name='ETOPO1_D2', 

69 base=_etopo1, 

70 ndeci=2, 

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

72 

73_etopo1_d4 = dataset.DecimatedTiledGlobalDataset( 

74 name='ETOPO1_D4', 

75 base=_etopo1_d2, 

76 ndeci=2, 

77 data_dir=op.join(topo_data_dir, 'ETOPO1_D4')) 

78 

79_etopo1_d8 = dataset.DecimatedTiledGlobalDataset( 

80 name='ETOPO1_D8', 

81 base=_etopo1_d4, 

82 ndeci=2, 

83 data_dir=op.join(topo_data_dir, 'ETOPO1_D8')) 

84 

85_srtmgl3_all = [ 

86 _srtmgl3, 

87 _srtmgl3_d2, 

88 _srtmgl3_d4, 

89 _srtmgl3_d8] 

90 

91_etopo1_all = [ 

92 _etopo1, 

93 _etopo1_d2, 

94 _etopo1_d4, 

95 _etopo1_d8] 

96 

97_dems = _srtmgl3_all + _etopo1_all 

98 

99 

100def make_all_missing_decimated(): 

101 for dem in _dems: 

102 if isinstance(dem, dataset.DecimatedTiledGlobalDataset): 

103 dem.make_all_missing() 

104 

105 

106def comparison(region, dems=_dems): 

107 import matplotlib.pyplot as plt 

108 

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

110 

111 fig = plt.gcf() 

112 

113 for idem, dem_ in enumerate(dems): 

114 fig.add_subplot(len(dems), 1, idem+1) 

115 t = dem_.get(region) 

116 if t: 

117 plt.pcolormesh(t.x(), t.y(), t.data) 

118 plt.title(dem_.name) 

119 plt.xlim(west, east) 

120 plt.ylim(south, north) 

121 

122 plt.show() 

123 

124 

125class UnknownDEM(Exception): 

126 pass 

127 

128 

129def dem_names(): 

130 return [dem.name for dem in _dems] 

131 

132 

133def dem(dem_name): 

134 for dem in _dems: 

135 if dem.name == dem_name: 

136 return dem 

137 

138 raise UnknownDEM(dem_name) 

139 

140 

141def get(dem_name, region): 

142 return dem(dem_name).get(region) 

143 

144 

145def elevation(lat, lon): 

146 ''' 

147 Get elevation at given point. 

148 

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

150 ''' 

151 

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

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

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

155 return r 

156 

157 

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

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

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

161 

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

163 

164 ok = [] 

165 if kind == 'land': 

166 for dem in _srtmgl3_all[::step]: 

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

168 ok.append(dem.name) 

169 break 

170 

171 for dem in _etopo1_all[::step]: 

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

173 ok.append(dem.name) 

174 break 

175 

176 return ok 

177 

178 

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