How to plot seismic activity with Anaconda and Jupyter Notebooks
We will be using Anaconda as our distribution platform for this activity. To get started with Anaconda, download the distribution from Anaconda, install according to your OS instructions, make a new environment by clicking the create button in the bottom left, and then search, select, and install the following packages [matplotlib, basemap, pillow, pandas]
Next, go to the homepage of the Anaconda Navigator, make sure that the top left drop down menu has your newly created environment selected, and install Jupyter Notebook.
If you already have Anaconda setup, please make sure you’ve updated Matplotlib to version 3.0.2, or downgraded to 3.0.0, 3.0.1 introduces a bug with Basemap when setting attributes with .plot(), which we will be doing later. Reference the github issue here for more information. If you’re unsure how to go about updating conda and it’s packages, see this Stack Overflow post, although you may have to resort to using pip if the distro is not updated at the time of reading. (it has not been as of 11/17/2018)
# open your environment with cmd prompt/shell
pip install matplotlib --upgrade# conda version will still display what it's distro is set to
# feel free to check your version
import matplotlib
print(matplotlib.__version__)
Great, now select Environments on the left hand menu and you will see an icon on your newly created environment, you will see several options, the one you want is “Open with Jupyter Notebook”
Go ahead and run it, and a new window will open in your default web browser listing your root directory. Go ahead and pick a place you want to store your Jupyter Notebook by navigating through your folders and create a new notebook by clicking the “new” button on the menu above your folders. If you followed the instructions correctly you should see something like this.
Awesome, you’ve made your first notebook. Now we can start having fun! Go ahead and open the notebook. You’ll be treated to an interactive shell, and I encourage you to read up on using the Jupyter Notebook if you’re not familiar with doing so.
First, we will create a basic orthographic map with color. For a list of all supported projections see here.
# import dependencies
# note that Basemap will be supported until 2020 only
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np# make sure to set your plot size before you start rendering to screen
plt.figure(figsize=(8, 8))# by default, only crude and low resolutions are installed, if you wanted
# higher fidelity images, refer to the documentation.default_map = Basemap(projection='ortho', lat_0=45, lon_0=-105,
resolution='l', area_thresh=1000.0)default_map.drawcoastlines()
default_map.drawcountries()
# matplotlib provides color creation
default_map.fillcontinents(color='navajowhite')
default_map.drawmapboundary()plt.show()
Next, let’s add lines of latitude and longitude by adding the following code before plt.show()
default_map.drawmeridians(np.arange(0, 360, 30))
default_map.drawparallels(np.arange(-90, 90, 30))
This is starting to look quite handsome! Sadly, as pretty as it is, it may not be the best projection to see seismic activity on a global scale. I’d encourage you to play around with projections at this point and see which one you like! As for myself, I’ve chosen to go with a Hammer projection and I’ve adjusted the lat/long values to my liking.
# change to hammer projection, fix lat/long
# import dependencies
# note that Basemap will be supported until 2020 only
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np# make sure to set your plot size before you start rendering to screen
plt.figure(figsize=(18, 18))# by default, only crude and low resolutions are installed, if you wanted
# higher fidelity images, refer to the documentation.default_map = Basemap(projection='hammer', lat_0=0, lon_0=-100,
resolution='l', area_thresh=1000.0)default_map.drawcoastlines()
default_map.drawcountries()
# matplotlib provides color creation
default_map.fillcontinents(color='navajowhite')
default_map.drawmapboundary()default_map.drawmeridians(np.arange(0, 360, 30))
default_map.drawparallels(np.arange(-90, 90, 30))plt.show()
I know what you’re thinking at this point! “Look ma, I can see my house from here.” Well, why don’t we try zooming in to our own geographical regions and taking a closer look. The easiest way of doing so would be changing your projection type to a “Mercator Projection”, which supports zooming. I looked up my long, lat coordinates online, and adjusted outwards to get a clearer picture. Go ahead and install basemap-data-hires from Anaconda as well if you want higher resolution for the zoomed in maps.
For this next step, I went ahead and took the liberty of throwing the merc projector into a function that accepts two integers, to create a regional map maker. Now there might be some issues if you exceed the lat/long limits, but you can easily adjust your inputs to get around this. This is quick and dirty code with no exception handling.
# what if we want to look at a specific area of the world?from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as npdef create_merc(lat, long):
plt.figure(figsize=(10, 10))# make sure the value of resolution is a lowercase L,
# for 'low', not a numeral 1
my_map = Basemap(projection='merc', lat_0=(lat-20), lon_0=(long+15),
resolution = 'h', area_thresh = .1,
# lower left hand corner longitude, lat
# upper right hand corner long, lat
llcrnrlon=(long-10), llcrnrlat=(lat-10),
urcrnrlon=(long+10), urcrnrlat=(lat+10))my_map.drawcoastlines()
my_map.drawcountries()
my_map.fillcontinents(color='navajowhite', lake_color='paleturquoise')
my_map.drawmapboundary()my_map.drawmeridians(np.arange(0, 360, 30))
my_map.drawparallels(np.arange(-90, 90, 30))plt.show()create_merc(39, -76)
Wow. Can we just appreciate all the hard work that went into these libraries to allow us to make something like this in just a few hours of reading through documentation? (and admittedly many more hours for me fighting to overcome some bugs)
Alright, moving on! Let’s try plotting specific coordinates. Like where I live currently for instance.
# plot the initial coordinate input
x, y = merc_map(long, lat)
merc_map.plot(x, y, 'go', markersize=12)plt.show()
Okay, we’ve proven we can plot figures. But we are probably going to want multiple, so let’s remove the assignment from the function inputs, and use a container to store what we want (remember, we want to read in data and store it to be plotted later) Let’s also go ahead and add some labels for our plots. Same idea here as the coordinates! Let’s test it manually first. (and a little bit of goofing off with yours truly)
# what if we want to look at a specific area of the world?from mpl_toolkits.basemap import Basemap
import matplotlib
import matplotlib.pyplot as plt
import numpy as npdef create_merc(lat, long):
plt.figure(figsize=(10, 10), dpi=100)
# increase font size
matplotlib.rcParams.update({'font.size': 16})# make sure the value of resolution is a lowercase L,
# for 'low', not a numeral 1
merc_map = Basemap(projection='merc', lat_0=(lat-20), lon_0= (long+15),
resolution = 'h', area_thresh = .1,
# lower left hand corner longitude, lat
# upper right hand corner long, lat
llcrnrlon=(long-10), llcrnrlat=(lat-10),
urcrnrlon=(long+10), urcrnrlat=(lat+10))
longs = [-77, -74]
lats = [39, 40]
x, y = merc_map(longs, lats)
merc_map.plot(x, y, 'go', markersize=8)
locations = ['Columbia, MD', 'New York City']
for label, xc, yc in zip(locations, x, y):
plt.text(xc-100000, yc+55000, label, color='teal')
merc_map.drawcoastlines(color='grey', linestyle=':')
merc_map.drawrivers(linewidth=0.3, color="cornflowerblue")
merc_map.fillcontinents(color='navajowhite', lake_color='paleturquoise')
merc_map.drawcountries()
merc_map.drawmapboundary() merc_map.drawmeridians(np.arange(0, 360, 30))
merc_map.drawparallels(np.arange(-90, 90, 30))
plt.show()create_merc(39, -77)
Now that we’ve explored mapping, plotting, and annotating, let’s start pulling in some data. For the purpose of this tutorial, i’ve downloaded the seismic data into a CSV format and hosted it on github for you to access in it’s raw format.
Go ahead and load the raw data into pandas and start exploring the data. What we will actually want to plot seismic activity with is it’s geographic location, so check for any missing values in the columns latitude and longitude.
import pandas as pd
seismic_activity_30_days = pd.read_csv('https://raw.githubusercontent.com/edwardauron/working_data/master/all_week.csv')seismic_activity_30_days.head()print(seismic_activity_30_days['latitude'].isnull().values.any())
print(seismic_activity_30_days['longitude'].isnull().values.any())False
False
Wonderful, there is no missing lat or long data, so we can skip cleaning, and go straight into wrangling data from the columns into our previously defined lists in our function. We reuse our existing code from our previous function to take a dataframe, and add some ways to handle it.
def create_hammer(df):
plt.figure(figsize=(10, 10), dpi=100)
# increase font size
matplotlib.rcParams.update({'font.size': 16})# make sure the value of resolution is a lowercase L,
# for 'low', not a numeral 1
merc_map = Basemap(projection='hammer', lat_0=0, lon_0=-130,
resolution = 'l', area_thresh = 1000)
working_df = df
longs = working_df['longitude'].tolist()
lats = working_df['latitude'].tolist()
x, y = merc_map(longs, lats)
merc_map.plot(x, y, marker='o', color='lavenderblush', markersize=4)
merc_map.drawcoastlines(color='grey', linestyle=':')
merc_map.drawrivers(linewidth=0.3, color="cornflowerblue")
merc_map.fillcontinents(color='navajowhite', lake_color='paleturquoise')
merc_map.drawcountries()
merc_map.drawstates(linewidth=0.3, color='saddlebrown')
merc_map.drawmapboundary()merc_map.drawmeridians(np.arange(0, 360, 30))
merc_map.drawparallels(np.arange(-90, 90, 30))
plt.show()create_hammer(seismic_activity_30_days)
If you’re following along you’ll notice that this creates a hot mess….
We need to plot all these points uniquely rather than all at once. Yup, it’s time for a for loop. Also it turns out “lavenderblush” is not the easiest color to see!
def create_hammer(df):
plt.figure(figsize=(10, 10), dpi=100)
# increase font size
matplotlib.rcParams.update({'font.size': 16})# make sure the value of resolution is a lowercase L,
# for 'low', not a numeral 1
hammer_map = Basemap(projection='hammer', lat_0=0, lon_0=-130,
resolution = 'l', area_thresh = 1000)
working_df = df
longs = working_df['longitude'].tolist()
lats = working_df['latitude'].tolist()
min_marker_size = 2.5
for longitude, latitude in zip(longs, lats):
x, y = hammer_map(longitude, latitude)
marker_size = min_marker_size
hammer_map.plot(x, y, marker='o', color='maroon', markersize=marker_size)
hammer_map.drawcoastlines(color='grey', linestyle=':')
hammer_map.drawrivers(linewidth=0.3, color="cornflowerblue")
hammer_map.fillcontinents(color='navajowhite', lake_color='paleturquoise')
hammer_map.drawcountries()
hammer_map.drawstates(linewidth=0.3, color='saddlebrown')
hammer_map.drawmapboundary()hammer_map.drawmeridians(np.arange(0, 360, 30))
hammer_map.drawparallels(np.arange(-90, 90, 30))
plt.show()
Now that’s more like it!
There is one final thing that we need to do. Can you guess? We need a way to differentiate these plot points! Can you find any columns in the dataframe that might help us? What sort of expression should we map that column to?Feel free to explore and come up with an idea before looking at my solution.
Well, here is my finished product. It still needs some refining, like a legend. But I went ahead and added a formatted title, as well as some toying around with the plotted circles sizes. It could use some hover events, and I plan to add an option to zoom in and project a user designated area with mercator projections, but I feel satisfied with what I have here. For now.
# plotting seismic activityfrom mpl_toolkits.basemap import Basemap
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import pandas as pddef create_hammer(df):
plt.figure(figsize=(18, 12), dpi=100)
# increase font size
matplotlib.rcParams.update({'font.size': 16})# make sure the value of resolution is a lowercase L,
# for 'low', not a numeral 1
hammer_map = Basemap(projection='hammer', lat_0=0, lon_0=-130,
resolution = 'l', area_thresh = 1000)
working_df = df
# generate data for plotting from our df
longs = working_df['longitude'].tolist()
lats = working_df['latitude'].tolist()
mags = working_df['mag'].tolist()
# handle dates for title
start = working_df.time.iloc[0]
end = working_df.time.iloc[1825]
min_marker_size = 2.5
for longitude, latitude, magnitude in zip(longs, lats, mags):
x, y = hammer_map(longitude, latitude)
marker_size = min_marker_size * magnitude\
if magnitude < 3.0:
hammer_map.plot(x, y, marker='o', color='maroon', markersize=marker_size)
elif magnitude < 5.0:
hammer_map.plot(x, y, marker='o', color='slateblue', markersize=marker_size)
else:
hammer_map.plot(x, y, marker='o', color='goldenrod', markersize=marker_size)
hammer_map.drawcoastlines(color='grey', linestyle=':')
hammer_map.drawrivers(linewidth=0.3, color="cornflowerblue")
hammer_map.fillcontinents(color='navajowhite', lake_color='paleturquoise')
hammer_map.drawcountries()
hammer_map.drawstates(linewidth=0.3, color='saddlebrown')
hammer_map.drawmapboundary()hammer_map.drawmeridians(np.arange(0, 360, 30))
hammer_map.drawparallels(np.arange(-90, 90, 30))
title = "Seismic activity with Magnitudes more than 1.0\n"
# string manipulation, list slicing
title += "%s through %s" % (start[:10], end[:10])
plt.title(title)
plt.show()create_hammer(seismic_activity_30_days)
