1#!/usr/bin/env python 

2# http://pyrocko.org - GPLv3 

3# 

4# The Pyrocko Developers, 21st Century 

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

6 

7import sys 

8import logging 

9import numpy as num 

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

11from pyrocko import moment_tensor as pmt 

12from pyrocko.plot import automap 

13 

14from optparse import OptionParser 

15 

16km = 1000. 

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

18 

19program_name = 'automap' 

20 

21usage = ''' 

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

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

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

25 

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

27 

28 

29def latlon_arrays(locs): 

30 return num.array( 

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

32 

33 

34def main(args=None): 

35 if args is None: 

36 args = sys.argv[1:] 

37 

38 parser = OptionParser( 

39 usage=usage, 

40 description=description) 

41 

42 parser.add_option( 

43 '--width', 

44 dest='width', 

45 type='float', 

46 default=20.0, 

47 metavar='FLOAT', 

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

49 

50 parser.add_option( 

51 '--height', 

52 dest='height', 

53 type='float', 

54 default=15.0, 

55 metavar='FLOAT', 

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

57 

58 parser.add_option( 

59 '--topo-resolution-min', 

60 dest='topo_resolution_min', 

61 type='float', 

62 default=40.0, 

63 metavar='FLOAT', 

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

65 

66 parser.add_option( 

67 '--topo-resolution-max', 

68 dest='topo_resolution_max', 

69 type='float', 

70 default=200.0, 

71 metavar='FLOAT', 

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

73 

74 parser.add_option( 

75 '--no-grid', 

76 dest='show_grid', 

77 default=True, 

78 action='store_false', 

79 help='don\'t show grid lines') 

80 

81 parser.add_option( 

82 '--no-topo', 

83 dest='show_topo', 

84 default=True, 

85 action='store_false', 

86 help='don\'t show topography') 

87 

88 parser.add_option( 

89 '--no-cities', 

90 dest='show_cities', 

91 default=True, 

92 action='store_false', 

93 help='don\'t show cities') 

94 

95 parser.add_option( 

96 '--no-illuminate', 

97 dest='illuminate', 

98 default=True, 

99 action='store_false', 

100 help='deactivate artificial illumination of topography') 

101 

102 parser.add_option( 

103 '--illuminate-factor-land', 

104 dest='illuminate_factor_land', 

105 type='float', 

106 metavar='FLOAT', 

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

108 

109 parser.add_option( 

110 '--illuminate-factor-ocean', 

111 dest='illuminate_factor_ocean', 

112 type='float', 

113 metavar='FLOAT', 

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

115 

116 parser.add_option( 

117 '--theme', 

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

119 default='topo', 

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

121 

122 parser.add_option( 

123 '--download-etopo1', 

124 dest='download_etopo1', 

125 action='store_true', 

126 help='download full ETOPO1 topography dataset') 

127 

128 parser.add_option( 

129 '--download-srtmgl3', 

130 dest='download_srtmgl3', 

131 action='store_true', 

132 help='download full SRTMGL3 topography dataset') 

133 

134 parser.add_option( 

135 '--make-decimated-topo', 

136 dest='make_decimated', 

137 action='store_true', 

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

139 

140 parser.add_option( 

141 '--stations', 

142 dest='stations_fn', 

143 metavar='FILENAME', 

144 help='load station coordinates from FILENAME') 

145 

146 parser.add_option( 

147 '--events', 

148 dest='events_fn', 

149 metavar='FILENAME', 

150 help='load event coordinates from FILENAME') 

151 

152 parser.add_option( 

153 '--debug', 

154 dest='debug', 

155 action='store_true', 

156 default=False, 

157 help='print debugging information to stderr') 

158 

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

160 

161 if options.debug: 

162 util.setup_logging(program_name, 'debug') 

163 else: 

164 util.setup_logging(program_name, 'info') 

165 

166 if options.download_etopo1: 

167 import pyrocko.datasets.topo.etopo1 

168 pyrocko.datasets.topo.etopo1.download() 

169 

170 if options.download_srtmgl3: 

171 import pyrocko.datasets.topo.srtmgl3 

172 pyrocko.datasets.topo.srtmgl3.download() 

173 

174 if options.make_decimated: 

175 import pyrocko.datasets.topo 

176 pyrocko.datasets.topo.make_all_missing_decimated() 

177 

178 if (options.download_etopo1 or options.download_srtmgl3 or 

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

180 

181 sys.exit(0) 

182 

183 if options.theme == 'soft': 

184 color_kwargs = { 

185 'illuminate_factor_land': options.illuminate_factor_land or 0.2, 

186 'illuminate_factor_ocean': options.illuminate_factor_ocean or 0.15, 

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

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

189 'topo_cpt_wet': 'light_sea_uniform', 

190 'topo_cpt_dry': 'light_land_uniform'} 

191 elif options.theme == 'topo': 

192 color_kwargs = { 

193 'illuminate_factor_land': options.illuminate_factor_land or 0.5, 

194 'illuminate_factor_ocean': options.illuminate_factor_ocean or 0.25} 

195 

196 events = [] 

197 if options.events_fn: 

198 events = model.load_events(options.events_fn) 

199 

200 stations = [] 

201 

202 if options.stations_fn: 

203 stations = model.load_stations(options.stations_fn) 

204 

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

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

207 

208 parser.print_help() 

209 sys.exit(1) 

210 

211 if len(args) == 4: 

212 try: 

213 lat = float(args[0]) 

214 lon = float(args[1]) 

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

216 except Exception: 

217 parser.print_help() 

218 sys.exit(1) 

219 else: 

220 lats, lons = latlon_arrays(stations+events) 

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

222 radius = float( 

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

224 radius *= 1.1 

225 

226 m = automap.Map( 

227 width=options.width, 

228 height=options.height, 

229 lat=lat, 

230 lon=lon, 

231 radius=radius, 

232 topo_resolution_max=options.topo_resolution_max, 

233 topo_resolution_min=options.topo_resolution_min, 

234 show_topo=options.show_topo, 

235 show_grid=options.show_grid, 

236 illuminate=options.illuminate, 

237 **color_kwargs) 

238 

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

240 

241 if options.show_cities: 

242 m.draw_cities() 

243 

244 if stations: 

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

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

247 

248 m.gmt.psxy( 

249 in_columns=(lons, lats), 

250 S='t8p', 

251 G='black', 

252 *m.jxyr) 

253 

254 for s in stations: 

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

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

257 

258 if events: 

259 beachball_symbol = 'mt' 

260 beachball_size = 20.0 

261 for ev in events: 

262 if ev.moment_tensor is None: 

263 m.gmt.psxy( 

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

265 S='c12p', 

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

267 W='1p,black', 

268 *m.jxyr) 

269 

270 else: 

271 devi = ev.moment_tensor.deviatoric() 

272 mt = devi.m_up_south_east() 

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

274 * pmt.magnitude_to_moment(5.0) 

275 m6 = pmt.to6(mt) 

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

277 

278 if m.gmt.is_gmt5(): 

279 kwargs = dict( 

280 M=True, 

281 S='%s%g' % ( 

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

283 else: 

284 kwargs = dict( 

285 S='%s%g' % ( 

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

287 

288 m.gmt.psmeca( 

289 in_rows=[data], 

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

291 E='white', 

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

293 *m.jxyr, 

294 **kwargs) 

295 

296 m.save(args[-1]) 

297 

298 

299if __name__ == '__main__': 

300 main(sys.argv[1:])