1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6import sys 

7import logging 

8import numpy as num 

9from pyrocko import util, model, orthodrome as od, gmtpy 

10from pyrocko import moment_tensor as pmt 

11from pyrocko.plot import automap 

12 

13from optparse import OptionParser 

14 

15km = 1000. 

16logger = logging.getLogger('pyrocko.apps.automap') 

17 

18program_name = 'automap' 

19 

20usage = ''' 

21usage: %s [options] [--] <lat> <lon> <radius_km> <output.(pdf|png)> 

22 %s [--download-etopo1] [--download-srtmgl3] 

23'''.strip() % (program_name, program_name) 

24 

25description = '''Create a simple map with topography.''' 

26 

27 

28def latlon_arrays(locs): 

29 return num.array( 

30 [(x.effective_lat, x.effective_lon) for x in locs]).T 

31 

32 

33def main(args=None): 

34 if args is None: 

35 args = sys.argv[1:] 

36 

37 parser = OptionParser( 

38 usage=usage, 

39 description=description) 

40 

41 parser.add_option( 

42 '--width', 

43 dest='width', 

44 type='float', 

45 default=20.0, 

46 metavar='FLOAT', 

47 help='set width of output image [cm] (%default)') 

48 

49 parser.add_option( 

50 '--height', 

51 dest='height', 

52 type='float', 

53 default=15.0, 

54 metavar='FLOAT', 

55 help='set height of output image [cm] (%default)') 

56 

57 parser.add_option( 

58 '--topo-resolution-min', 

59 dest='topo_resolution_min', 

60 type='float', 

61 default=40.0, 

62 metavar='FLOAT', 

63 help='minimum resolution of topography [dpi] (%default)') 

64 

65 parser.add_option( 

66 '--topo-resolution-max', 

67 dest='topo_resolution_max', 

68 type='float', 

69 default=200.0, 

70 metavar='FLOAT', 

71 help='maximum resolution of topography [dpi] (%default)') 

72 

73 parser.add_option( 

74 '--no-grid', 

75 dest='show_grid', 

76 default=True, 

77 action='store_false', 

78 help="don't show grid lines") 

79 

80 parser.add_option( 

81 '--no-topo', 

82 dest='show_topo', 

83 default=True, 

84 action='store_false', 

85 help="don't show topography") 

86 

87 parser.add_option( 

88 '--no-cities', 

89 dest='show_cities', 

90 default=True, 

91 action='store_false', 

92 help="don't show cities") 

93 

94 parser.add_option( 

95 '--no-illuminate', 

96 dest='illuminate', 

97 default=True, 

98 action='store_false', 

99 help='deactivate artificial illumination of topography') 

100 

101 parser.add_option( 

102 '--illuminate-factor-land', 

103 dest='illuminate_factor_land', 

104 type='float', 

105 metavar='FLOAT', 

106 help='set factor for artificial illumination of land (0.5)') 

107 

108 parser.add_option( 

109 '--illuminate-factor-ocean', 

110 dest='illuminate_factor_ocean', 

111 type='float', 

112 metavar='FLOAT', 

113 help='set factor for artificial illumination of ocean (0.25)') 

114 

115 parser.add_option( 

116 '--theme', 

117 choices=['topo', 'soft'], 

118 default='topo', 

119 help='select color theme, available: topo, soft (topo)"') 

120 

121 parser.add_option( 

122 '--download-etopo1', 

123 dest='download_etopo1', 

124 action='store_true', 

125 help='download full ETOPO1 topography dataset') 

126 

127 parser.add_option( 

128 '--download-srtmgl3', 

129 dest='download_srtmgl3', 

130 action='store_true', 

131 help='download full SRTMGL3 topography dataset') 

132 

133 parser.add_option( 

134 '--make-decimated-topo', 

135 dest='make_decimated', 

136 action='store_true', 

137 help='pre-make all decimated topography datasets') 

138 

139 parser.add_option( 

140 '--stations', 

141 dest='stations_fn', 

142 metavar='FILENAME', 

143 help='load station coordinates from FILENAME') 

144 

145 parser.add_option( 

146 '--events', 

147 dest='events_fn', 

148 metavar='FILENAME', 

149 help='load event coordinates from FILENAME') 

150 

151 parser.add_option( 

152 '--debug', 

153 dest='debug', 

154 action='store_true', 

155 default=False, 

156 help='print debugging information to stderr') 

157 

158 (options, args) = parser.parse_args(args) 

159 

160 if options.debug: 

161 util.setup_logging(program_name, 'debug') 

162 else: 

163 util.setup_logging(program_name, 'info') 

164 

165 if options.download_etopo1: 

166 import pyrocko.dataset.topo.etopo1 

167 pyrocko.dataset.topo.etopo1.download() 

168 

169 if options.download_srtmgl3: 

170 import pyrocko.dataset.topo.srtmgl3 

171 pyrocko.dataset.topo.srtmgl3.download() 

172 

173 if options.make_decimated: 

174 import pyrocko.dataset.topo 

175 pyrocko.dataset.topo.make_all_missing_decimated() 

176 

177 if (options.download_etopo1 or options.download_srtmgl3 or 

178 options.make_decimated) and len(args) == 0: 

179 

180 sys.exit(0) 

181 

182 if options.theme == 'soft': 

183 color_kwargs = { 

184 'illuminate_factor_land': options.illuminate_factor_land or 0.2, 

185 'illuminate_factor_ocean': options.illuminate_factor_ocean or 0.15, 

186 'color_wet': (216, 242, 254), 

187 'color_dry': (238, 236, 230), 

188 'topo_cpt_wet': 'light_sea_uniform', 

189 'topo_cpt_dry': 'light_land_uniform'} 

190 elif options.theme == 'topo': 

191 color_kwargs = { 

192 'illuminate_factor_land': options.illuminate_factor_land or 0.5, 

193 'illuminate_factor_ocean': options.illuminate_factor_ocean or 0.25} 

194 

195 events = [] 

196 if options.events_fn: 

197 events = model.load_events(options.events_fn) 

198 

199 stations = [] 

200 

201 if options.stations_fn: 

202 stations = model.load_stations(options.stations_fn) 

203 

204 if not (len(args) == 4 or ( 

205 len(args) == 1 and (events or stations))): 

206 

207 parser.print_help() 

208 sys.exit(1) 

209 

210 if len(args) == 4: 

211 try: 

212 lat = float(args[0]) 

213 lon = float(args[1]) 

214 radius = float(args[2])*km 

215 except Exception: 

216 parser.print_help() 

217 sys.exit(1) 

218 else: 

219 lats, lons = latlon_arrays(stations+events) 

220 lat, lon = map(float, od.geographic_midpoint(lats, lons)) 

221 radius = float( 

222 num.max(od.distance_accurate50m_numpy(lat, lon, lats, lons))) 

223 radius *= 1.1 

224 

225 m = automap.Map( 

226 width=options.width, 

227 height=options.height, 

228 lat=lat, 

229 lon=lon, 

230 radius=radius, 

231 topo_resolution_max=options.topo_resolution_max, 

232 topo_resolution_min=options.topo_resolution_min, 

233 show_topo=options.show_topo, 

234 show_grid=options.show_grid, 

235 illuminate=options.illuminate, 

236 **color_kwargs) 

237 

238 logger.debug('map configuration:\n%s' % str(m)) 

239 

240 if options.show_cities: 

241 m.draw_cities() 

242 

243 if stations: 

244 lats = [s.lat for s in stations] 

245 lons = [s.lon for s in stations] 

246 

247 m.gmt.psxy( 

248 in_columns=(lons, lats), 

249 S='t8p', 

250 G='black', 

251 *m.jxyr) 

252 

253 for s in stations: 

254 m.add_label(s.lat, s.lon, '%s' % '.'.join( 

255 x for x in s.nsl() if x)) 

256 

257 if events: 

258 beachball_symbol = 'mt' 

259 beachball_size = 20.0 

260 for ev in events: 

261 if ev.moment_tensor is None: 

262 m.gmt.psxy( 

263 in_rows=[[ev.lon, ev.lat]], 

264 S='c12p', 

265 G=gmtpy.color('scarletred2'), 

266 W='1p,black', 

267 *m.jxyr) 

268 

269 else: 

270 devi = ev.moment_tensor.deviatoric() 

271 mt = devi.m_up_south_east() 

272 mt = mt / ev.moment_tensor.scalar_moment() \ 

273 * pmt.magnitude_to_moment(5.0) 

274 m6 = pmt.to6(mt) 

275 data = (ev.lon, ev.lat, 10) + tuple(m6) + (1, 0, 0) 

276 

277 if m.gmt.is_gmt5(): 

278 kwargs = dict( 

279 M=True, 

280 S='%s%g' % ( 

281 beachball_symbol[0], beachball_size / gmtpy.cm)) 

282 else: 

283 kwargs = dict( 

284 S='%s%g' % ( 

285 beachball_symbol[0], beachball_size*2 / gmtpy.cm)) 

286 

287 m.gmt.psmeca( 

288 in_rows=[data], 

289 G=gmtpy.color('chocolate1'), 

290 E='white', 

291 W='1p,%s' % gmtpy.color('chocolate3'), 

292 *m.jxyr, 

293 **kwargs) 

294 

295 m.save(args[-1]) 

296 

297 

298if __name__ == '__main__': 

299 main()