PGA绘全球地震分布
Contents
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 |
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 19 21:11:41 2015
@author: aicesun
"""
from obspy import readEvents
from obspy.fdsn import Client
import datetime
# Open the earthquake data file.
catalog = readEvents("http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month.quakeml")
# Create empty lists for the data we are interested in.
lats, lons = [], []
magnitudes = []
timestrings = []
# Read through the entire file, skip the first line,
# and pull out just the lats and lons.
for item in range(len(catalog)):
magnitudes.append(catalog[item].magnitudes[0].mag)
lats.append(catalog[item].origins[0].latitude)
lons.append(catalog[item].origins[0].longitude)
timestrings.append(catalog[item].origins[0].time)
# --- Build Map ---
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
def get_marker_color(magnitude):
# Returns green for small earthquakes, yellow for moderate
# earthquakes, and red for significant earthquakes.
if magnitude < 3.0:
return ('go')
elif magnitude < 5.0:
return ('yo')
else:
return ('ro')
# Make this plot larger.
plt.figure(figsize=(16,12))
eq_map = Basemap(projection='robin', resolution = 'l', area_thresh = 1000.0,lat_0=0, lon_0=-130)
eq_map.drawcoastlines()
eq_map.drawcountries()
eq_map.fillcontinents(color = 'gray')
eq_map.drawmapboundary()
eq_map.drawmeridians(np.arange(0, 360, 30))
eq_map.drawparallels(np.arange(-90, 90, 30))
min_marker_size = 2.5
for lon, lat, mag in zip(lons, lats, magnitudes):
x,y = eq_map(lon, lat)
msize = mag * min_marker_size
marker_string = get_marker_color(mag)
eq_map.plot(x, y, marker_string, markersize=msize)
starttime = timestrings[-1].date
endtime = timestrings[0].date
title_string = "Earthquakes of Magnitude 2.5 or Greater\n"
title_string += "%s through %s" % (starttime, endtime)
plt.title(title_string)
plt.show() |
Author aice
LastMod 2017-09-20