pga绘全球地震分布

 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 Greatern"
	title_string += "%s through %s" % (starttime, endtime)
	plt.title(title_string)
	
	plt.show()