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