Source code for pyrocko.apps.automap

# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------

'''
Automap - create simple maps with GMT.
'''

import sys
import logging
import numpy as num
from pyrocko import util, model, orthodrome as od, gmtpy
from pyrocko import moment_tensor as pmt
from pyrocko.plot import automap

from optparse import OptionParser

km = 1000.
logger = logging.getLogger('pyrocko.apps.automap')

program_name = 'automap'

usage = '''
usage: %s [options] [--] <lat> <lon> <radius_km> <output.(pdf|png)>
       %s [--download-etopo1] [--download-srtmgl3]
'''.strip() % (program_name, program_name)

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


def latlon_arrays(locs):
    return num.array(
        [(x.effective_lat, x.effective_lon) for x in locs]).T


[docs]def main(args=None): ''' CLI entry point for Pyrocko's ``automap`` app. ''' if args is None: args = sys.argv[1:] parser = OptionParser( usage=usage, description=description) parser.add_option( '--width', dest='width', type='float', default=20.0, metavar='FLOAT', help='set width of output image [cm] (%default)') parser.add_option( '--height', dest='height', type='float', default=15.0, metavar='FLOAT', help='set height of output image [cm] (%default)') parser.add_option( '--topo-resolution-min', dest='topo_resolution_min', type='float', default=40.0, metavar='FLOAT', help='minimum resolution of topography [dpi] (%default)') parser.add_option( '--topo-resolution-max', dest='topo_resolution_max', type='float', default=200.0, metavar='FLOAT', help='maximum resolution of topography [dpi] (%default)') parser.add_option( '--no-grid', dest='show_grid', default=True, action='store_false', help="don't show grid lines") parser.add_option( '--no-topo', dest='show_topo', default=True, action='store_false', help="don't show topography") parser.add_option( '--no-cities', dest='show_cities', default=True, action='store_false', help="don't show cities") parser.add_option( '--no-illuminate', dest='illuminate', default=True, action='store_false', help='deactivate artificial illumination of topography') parser.add_option( '--illuminate-factor-land', dest='illuminate_factor_land', type='float', metavar='FLOAT', help='set factor for artificial illumination of land (0.5)') parser.add_option( '--illuminate-factor-ocean', dest='illuminate_factor_ocean', type='float', metavar='FLOAT', help='set factor for artificial illumination of ocean (0.25)') parser.add_option( '--theme', choices=['topo', 'soft'], default='topo', help='select color theme, available: topo, soft (topo)"') parser.add_option( '--download-etopo1', dest='download_etopo1', action='store_true', help='download full ETOPO1 topography dataset') parser.add_option( '--download-srtmgl3', dest='download_srtmgl3', action='store_true', help='download full SRTMGL3 topography dataset') parser.add_option( '--make-decimated-topo', dest='make_decimated', action='store_true', help='pre-make all decimated topography datasets') parser.add_option( '--stations', dest='stations_fn', metavar='FILENAME', help='load station coordinates from FILENAME') parser.add_option( '--events', dest='events_fn', metavar='FILENAME', help='load event coordinates from FILENAME') parser.add_option( '--debug', dest='debug', action='store_true', default=False, help='print debugging information to stderr') (options, args) = parser.parse_args(args) if options.debug: util.setup_logging(program_name, 'debug') else: util.setup_logging(program_name, 'info') if options.download_etopo1: import pyrocko.dataset.topo.etopo1 pyrocko.dataset.topo.etopo1.download() if options.download_srtmgl3: import pyrocko.dataset.topo.srtmgl3 pyrocko.dataset.topo.srtmgl3.download() if options.make_decimated: import pyrocko.dataset.topo pyrocko.dataset.topo.make_all_missing_decimated() if (options.download_etopo1 or options.download_srtmgl3 or options.make_decimated) and len(args) == 0: sys.exit(0) if options.theme == 'soft': color_kwargs = { 'illuminate_factor_land': options.illuminate_factor_land or 0.2, 'illuminate_factor_ocean': options.illuminate_factor_ocean or 0.15, 'color_wet': (216, 242, 254), 'color_dry': (238, 236, 230), 'topo_cpt_wet': 'light_sea_uniform', 'topo_cpt_dry': 'light_land_uniform'} elif options.theme == 'topo': color_kwargs = { 'illuminate_factor_land': options.illuminate_factor_land or 0.5, 'illuminate_factor_ocean': options.illuminate_factor_ocean or 0.25} events = [] if options.events_fn: events = model.load_events(options.events_fn) stations = [] if options.stations_fn: stations = model.load_stations(options.stations_fn) if not (len(args) == 4 or ( len(args) == 1 and (events or stations))): parser.print_help() sys.exit(1) if len(args) == 4: try: lat = float(args[0]) lon = float(args[1]) radius = float(args[2])*km except Exception: parser.print_help() sys.exit(1) else: lats, lons = latlon_arrays(stations+events) lat, lon = map(float, od.geographic_midpoint(lats, lons)) radius = float( num.max(od.distance_accurate50m_numpy(lat, lon, lats, lons))) radius *= 1.1 m = automap.Map( width=options.width, height=options.height, lat=lat, lon=lon, radius=radius, topo_resolution_max=options.topo_resolution_max, topo_resolution_min=options.topo_resolution_min, show_topo=options.show_topo, show_grid=options.show_grid, illuminate=options.illuminate, **color_kwargs) logger.debug('map configuration:\n%s' % str(m)) if options.show_cities: m.draw_cities() if stations: lats = [s.lat for s in stations] lons = [s.lon for s in stations] m.gmt.psxy( in_columns=(lons, lats), S='t8p', G='black', *m.jxyr) for s in stations: m.add_label(s.lat, s.lon, '%s' % '.'.join( x for x in s.nsl() if x)) if events: beachball_symbol = 'mt' beachball_size = 20.0 for ev in events: if ev.moment_tensor is None: m.gmt.psxy( in_rows=[[ev.lon, ev.lat]], S='c12p', G=gmtpy.color('scarletred2'), W='1p,black', *m.jxyr) else: devi = ev.moment_tensor.deviatoric() mt = devi.m_up_south_east() mt = mt / ev.moment_tensor.scalar_moment() \ * pmt.magnitude_to_moment(5.0) m6 = pmt.to6(mt) data = (ev.lon, ev.lat, 10) + tuple(m6) + (1, 0, 0) if m.gmt.is_gmt5(): kwargs = dict( M=True, S='%s%g' % ( beachball_symbol[0], beachball_size / gmtpy.cm)) else: kwargs = dict( S='%s%g' % ( beachball_symbol[0], beachball_size*2 / gmtpy.cm)) m.gmt.psmeca( in_rows=[data], G=gmtpy.color('chocolate1'), E='white', W='1p,%s' % gmtpy.color('chocolate3'), *m.jxyr, **kwargs) m.save(args[-1])
if __name__ == '__main__': main()