Assumptions and Expectations:
During this demo, you will be exposed to:
What is IPython?
Those interfaces can include:
You'll need git installed - this may require administrator privileges. Recent versions of Mac OSX come with git pre-installed. Type
which git
in a Terminal window to see if you have this software already.
Next you need Python, with some additional libraries that makes it interesting for scientists.
I recommend Anaconda from Continuum Analytics. http://continuum.io/downloads
There are other sources for "scientific" Python distributions, listed here: http://www.scipy.org/install.html
You will also need to install the following packages from a command line (OSX Terminal or Windows Command Prompt):
Once Jupyter is installed, you need to start it up (from a Terminal or Command Prompt) in the directory where the IPythonDemo.ipynb file lives:
jupyter notebook
This is a markdown cell.
This is the second line.
Hello
This is a test
You cause the browser to render the markdown text in this cell by hitting "Shift-Return".
#This is a code cell
#You execute all of the code in these cells by hitting "Shift-Return".
foo = 5
print(foo)
print('goodbye')
5 goodbye
The following "magic" command (anything prefaced with "%" is referred to as magic) allows plots to be shown in-line with the code and markdown cells.
%matplotlib inline
#The import statement allows you to use code from other modules
import numpy as np
import matplotlib.pyplot as plt
x = np.random.rand(1000,1) #this results in a Numpy array - more on this later
y = np.random.rand(1000,1)
fig = plt.figure(figsize=(10,6))
plt.plot(x,y,'bx')
th = plt.title('Title')
xl = plt.xlabel('X')
yl = plt.ylabel('Y')
import mpld3
fig = plt.figure(figsize=(8,5))
plt.plot(x,y,'b.')
th = plt.title('Random Points')
xl = plt.xlabel('X')
yl = plt.ylabel('Y')
mpld3.display()
NB: The SciPy stack includes a package called pandas, which provides R-like data analysis and tools for reading in spreadsheet-like data (CSV, Excel, etc.)
import pandas as pd
df = pd.read_excel('EXPO_CAT_2007_12.xlsx',sheetname='Sheet1')
f = plt.figure(figsize=(10,6))
ph = plt.semilogy(df['PAGER_prefMag'].get_values(),df['PAGER_prefShakingDeaths'].get_values(),'r.')
th = plt.xlabel('Magnitude')
th = plt.ylabel('Shaking Deaths')
th = plt.title('Magnitude vs Shaking Deaths')
...and this is why we don't model fatalities this way!
pandas reads the spreadsheet into a data structure called a DataFrame. You can do many things with a DataFrame, the least of which is viewing the contents in a notebook, as below.
df
eqID | eqName | PAGER_dateStr | PAGER_dateNum | PAGER_lat | PAGER_lon | PAGER_depth | PAGER_prefMag | PAGER_prefMagType | ISO_country | ... | R055 | R060 | R065 | R070 | R075 | R080 | R085 | R090 | R095 | R100 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 196002292340 | Agadir, Morocco | 1960-02-29 23:40:12 | 715935.986250 | 30.500 | -9.500 | NaN | 6.8 | UK | MOROCCO | ... | 386373 | 157250 | 58444 | 72540 | 9437 | 2616 | 2043 | 1794 | 0 | 0 |
1 | 196005221911 | Concepcion, Chile | 1960-05-22 19:11:17 | 716018.799502 | -38.294 | -73.054 | 35.0 | 9.6 | Mw | CHILE | ... | 169057 | 147534 | 296107 | 642706 | 312787 | 3124 | 733 | 885 | 108 | 373 |
2 | 196209011920 | Buyin-Zara, Iran | 1962-09-01 19:20:40 | 716850.806019 | 35.563 | 49.817 | 19.0 | 6.9 | Ms | IRAN; ISLAMIC REPUBLIC OF | ... | 187662 | 71227 | 26724 | 43298 | 45794 | 12233 | 6586 | 672 | 0 | 0 |
3 | 196307260417 | Skopje, Yugoslavia | 1963-07-26 04:17:00 | 717178.178472 | 42.100 | 21.500 | 5.0 | 6.0 | UK | MACEDONIA; THE FORMER YUGOSLAV REPUBLIC OF | ... | 419853 | 101583 | 24669 | 9559 | 9775 | 560 | 0 | 0 | 0 | 0 |
4 | 196403280336 | Prince William Sound, Alaska | 1964-03-28 03:36:12 | 717424.150139 | 61.019 | -147.626 | 6.3 | 9.2 | Mw | UNITED STATES MINOR OUTLYING ISLANDS | ... | 4330 | 6425 | 20069 | 32395 | 8741 | 7207 | 0 | 0 | 0 | 0 |
5 | 196504291528 | Puget Sound, Washington | 1965-04-29 15:28:45 | 717821.644965 | 47.317 | -122.333 | 65.7 | 6.5 | mb | UNITED STATES | ... | 82770 | 88293 | 96935 | 44713 | 11616 | 53 | 0 | 0 | 0 | 0 |
6 | 196606280426 | Parkfield, California | 1966-06-28 04:26:15 | 718246.184896 | 35.875 | -120.487 | 1.2 | 6.1 | Mw | UNITED STATES | ... | 4261 | 2097 | 274 | 96 | 226 | 68 | 128 | 0 | 0 | 0 |
7 | 196608191222 | Varto, Turkey | 1966-08-19 12:22:10 | 718298.515394 | 39.165 | 41.575 | 17.0 | 6.8 | Mw | TURKEY | ... | 128780 | 117175 | 66620 | 33889 | 14404 | 6760 | 0 | 0 | 0 | 0 |
8 | 196702091524 | Huila, Colombia | 1967-02-09 15:24:47 | 718472.642211 | 2.890 | -74.801 | 41.4 | 7.2 | Mw | COLOMBIA | ... | 132270 | 93826 | 53462 | 24331 | 12691 | 5082 | 1744 | 84 | 0 | 0 |
9 | 196707221656 | Mudurnu Valley, Turkey | 1967-07-22 16:56:55 | 718635.706192 | 40.629 | 30.740 | 4.4 | 7.4 | Mw | TURKEY | ... | 332488 | 149770 | 94191 | 66201 | 37965 | 18016 | 19324 | 16157 | 8302 | 26912 |
10 | 196707300000 | Caracas, Venezuela | 1967-07-30 00:00:02 | 718643.000023 | 10.555 | -67.310 | 23.1 | 6.6 | Mw | VENEZUELA | ... | 53493 | 4545 | 6298 | 3922 | 302 | 28245 | 18778 | 0 | 0 | 0 |
11 | 196712102251 | Koyna, India | 1967-12-10 22:51:21 | 718776.952326 | 17.390 | 73.774 | 4.5 | 6.3 | Mw | INDIA | ... | 864854 | 95226 | 26255 | 19916 | 11493 | 9 | 0 | 0 | 0 | 0 |
12 | 196804090229 | Borrego Mountain, California | 1968-04-09 02:29:01 | 718897.103484 | 33.157 | -116.194 | 15.0 | 6.6 | Mw | UNITED STATES | ... | 75572 | 24267 | 6080 | 3851 | 1827 | 990 | 1213 | 268 | 112 | 132 |
13 | 196805231724 | Inangahua, New Zealand | 1968-05-23 17:24:20 | 718941.725231 | -41.743 | 172.123 | 46.8 | 7.2 | Mw | NEW ZEALAND | ... | 31282 | 4417 | 2108 | 1079 | 940 | 4926 | 567 | 416 | 64 | 381 |
14 | 196808311047 | Dasht-e Bayaz, Iran | 1968-08-31 10:47:40 | 719041.449769 | 34.037 | 58.962 | 11.8 | 7.2 | Mw | IRAN; ISLAMIC REPUBLIC OF | ... | 20499 | 84456 | 8283 | 3184 | 4054 | 3540 | 5106 | 208 | 0 | 0 |
15 | 196810140258 | Meckering, Australia | 1968-10-14 02:58:50 | 719085.124190 | -31.522 | 116.978 | 5.0 | 6.8 | Ms | AUSTRALIA | ... | 9869 | 3888 | 425 | 832 | 159 | 89 | 185 | 17 | 0 | 0 |
16 | 196907252249 | Yangjiang, China | 1969-07-25 22:49:38 | 719369.951134 | 21.611 | 111.803 | 3.7 | 5.9 | Ms | CHINA | ... | 468075 | 200547 | 120040 | 106790 | 48788 | 27734 | 36218 | 6716 | 0 | 0 |
17 | 196909292003 | Ceres, South Africa | 1969-09-29 20:03:30 | 719435.835764 | -33.191 | 19.335 | 15.0 | 6.3 | Ms | SOUTH AFRICA | ... | 44213 | 29145 | 9968 | 7396 | 8657 | 6712 | 2523 | 0 | 0 | 0 |
18 | 197001041700 | Tonghai, China | 1970-01-04 17:00:40 | 719532.708796 | 24.148 | 102.462 | 14.1 | 7.2 | Mw | CHINA | ... | 594048 | 389025 | 408284 | 355877 | 340600 | 135622 | 69758 | 81777 | 36177 | 36360 |
19 | 197003101715 | Calingiri, Australia | 1970-03-10 17:15:12 | 719597.718889 | -31.093 | 116.513 | 15.0 | 5.5 | mb | AUSTRALIA | ... | 404 | 14 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
20 | 197003282102 | Gediz, Turkey | 1970-03-28 21:02:25 | 719615.876678 | 39.174 | 29.545 | 24.4 | 6.9 | Mw | TURKEY | ... | 354498 | 196479 | 100403 | 70369 | 41779 | 43276 | 23391 | 16205 | 10199 | 17174 |
21 | 197005312023 | Peru | 1970-05-31 20:23:32 | 719679.849676 | -9.248 | -78.842 | 73.2 | 7.5 | mB | PERU | ... | 640684 | 428048 | 457396 | 361988 | 84788 | 27566 | 0 | 0 | 0 | 0 |
22 | 197012100434 | Peru-Ecuador Border | 1970-12-10 04:34:39 | 719872.190729 | -4.079 | -80.660 | 19.5 | 7.1 | Mw | PERU | ... | 209312 | 217192 | 487314 | 178338 | 163376 | 1418 | 0 | 0 | 0 | 0 |
23 | 197102091400 | San Fernando, California | 1971-02-09 14:00:41 | 719933.583808 | 34.400 | -118.391 | 5.0 | 6.6 | Mw | UNITED STATES | ... | 28862 | 13045 | 2422 | 881 | 425 | 270 | 426 | 219 | 0 | 0 |
24 | 197107090303 | Valparaiso, Chile | 1971-07-09 03:03:19 | 720083.127303 | -32.558 | -71.085 | 59.0 | 7.8 | Mw | CHILE | ... | 239605 | 139548 | 155303 | 122191 | 117201 | 173507 | 130 | 0 | 0 | 0 |
25 | 197204100206 | Ghir, Iran | 1972-04-10 02:06:51 | 720359.088090 | 28.411 | 52.789 | 5.6 | 6.9 | Ms | IRAN; ISLAMIC REPUBLIC OF | ... | 209924 | 17064 | 12987 | 13963 | 6194 | 11589 | 4919 | 1532 | 617 | 50 |
26 | 197212230629 | Managua, Nicaragua | 1972-12-23 06:29:44 | 720616.270648 | 12.351 | -86.127 | 6.8 | 6.2 | Ms | NICARAGUA | ... | 127586 | 12359 | 1554 | 14 | 0 | 0 | 0 | 0 | 0 | 0 |
27 | 197301031431 | NaN | 1973-01-03 14:31:07 | 720627.604942 | 39.163 | 71.915 | 42.5 | 5.5 | Mb | TAJIKISTAN | ... | 3818 | 43 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
28 | 197301061552 | NaN | 1973-01-06 15:52:42 | 720630.661597 | -14.653 | 166.484 | 27.1 | 6.1 | Mb | VANUATU | ... | 2538 | 222 | 659 | 319 | 915 | 0 | 0 | 0 | 0 | 0 |
29 | 197301130605 | NaN | 1973-01-13 06:05:40 | 720637.253935 | -16.817 | 28.137 | 10.0 | 5.0 | Mb | ZAMBIA | ... | 120 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
5607 | 200711251602 | NaN | 2007-11-25 16:02:15 | 733371.668229 | -8.292 | 118.370 | 20.0 | 6.5 | Mw | INDONESIA | ... | 282625 | 155074 | 172412 | 32902 | 15439 | 9371 | 2717 | 4107 | 0 | 0 |
5608 | 200711251741 | NaN | 2007-11-25 17:41:38 | 733371.737245 | -2.238 | 100.413 | 35.0 | 5.9 | Mw | INDONESIA | ... | 47355 | 3263 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5609 | 200711251953 | NaN | 2007-11-25 19:53:05 | 733371.828530 | -8.224 | 118.467 | 18.0 | 6.5 | Mw | INDONESIA | ... | 268828 | 286116 | 85785 | 3787 | 1520 | 6166 | 2369 | 858 | 0 | 0 |
5610 | 200711252312 | NaN | 2007-11-25 23:12:17 | 733371.966863 | 28.555 | 77.057 | 10.0 | 4.7 | Mb | INDIA | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5611 | 200711261351 | NaN | 2007-11-26 13:51:39 | 733372.577535 | 37.392 | 141.586 | 42.0 | 5.9 | Mw | JAPAN | ... | 4220 | 1496 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5612 | 200711270426 | NaN | 2007-11-27 04:26:58 | 733373.185394 | 16.058 | 119.838 | 35.0 | 5.9 | Mw | PHILIPPINES | ... | 448019 | 211334 | 166941 | 358146 | 0 | 0 | 0 | 0 | 0 | 0 |
5613 | 200711271149 | NaN | 2007-11-27 11:49:58 | 733373.493032 | -10.950 | 162.149 | 16.0 | 6.6 | Mw | SOLOMON ISLANDS | ... | 7808 | 3307 | 2939 | 2106 | 3839 | 547 | 417 | 0 | 0 | 0 |
5614 | 200712010144 | NaN | 2007-12-01 01:44:31 | 733377.072581 | 1.981 | 97.883 | 44.0 | 5.9 | Mw | INDONESIA | ... | 94266 | 131551 | 53313 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5615 | 200712032152 | NaN | 2007-12-03 21:52:45 | 733379.911632 | -6.374 | 148.463 | 29.0 | 5.9 | Mw | PAPUA NEW GUINEA | ... | 7950 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5616 | 200712032233 | NaN | 2007-12-03 22:33:13 | 733379.939734 | -6.452 | 148.623 | 37.0 | 5.8 | Mw | PAPUA NEW GUINEA | ... | 5365 | 598 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5617 | 200712062143 | NaN | 2007-12-06 21:43:47 | 733382.905405 | 12.263 | 125.386 | 37.0 | 5.9 | Mw | PHILIPPINES | ... | 157388 | 192098 | 107301 | 77402 | 0 | 0 | 0 | 0 | 0 | 0 |
5618 | 200712070141 | NaN | 2007-12-07 01:41:01 | 733383.070150 | -13.516 | -76.656 | 22.0 | 5.6 | Mw | PERU | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5619 | 200712081955 | NaN | 2007-12-08 19:55:19 | 733384.830081 | -7.558 | 37.646 | 6.0 | 5.6 | Mw | TANZANIA; UNITED REPUBLIC OF | ... | 49479 | 15839 | 9918 | 3041 | 596 | 0 | 0 | 0 | 0 | 0 |
5620 | 200712090203 | NaN | 2007-12-09 02:03:29 | 733385.085752 | -15.048 | -44.231 | 10.0 | 4.9 | Mb | BRAZIL | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5621 | 200712090522 | NaN | 2007-12-09 05:22:11 | 733385.223738 | 12.277 | 125.420 | 39.0 | 5.6 | Mw | PHILIPPINES | ... | 168342 | 53885 | 73126 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5622 | 200712101210 | NaN | 2007-12-10 12:10:07 | 733386.507025 | -3.568 | 144.439 | 10.0 | 5.7 | Mw | PAPUA NEW GUINEA | ... | 1495 | 274 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5623 | 200712130520 | NaN | 2007-12-13 05:20:21 | 733389.222465 | -23.157 | -70.479 | 15.0 | 6.0 | Mw | CHILE | ... | 3501 | 889 | 507 | 308 | 2131 | 753 | 21 | 0 | 0 | 0 |
5624 | 200712130723 | NaN | 2007-12-13 07:23:39 | 733389.308090 | -23.202 | -70.549 | 16.0 | 6.2 | Mw | CHILE | ... | 24791 | 2230 | 705 | 511 | 2554 | 118 | 112 | 0 | 0 | 0 |
5625 | 200712131452 | NaN | 2007-12-13 14:52:16 | 733389.619630 | 51.601 | 16.073 | 5.0 | 4.5 | Mb | POLAND | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5626 | 200712151822 | NaN | 2007-12-15 18:22:27 | 733391.765590 | -32.689 | -71.695 | 25.0 | 5.9 | Mw | CHILE | ... | 60977 | 24569 | 24868 | 26524 | 18413 | 0 | 0 | 0 | 0 | 0 |
5627 | 200712160809 | NaN | 2007-12-16 08:09:17 | 733392.339780 | -22.954 | -70.182 | 45.0 | 6.7 | Mw | CHILE | ... | 10307 | 30943 | 6858 | 4827 | 3479 | 0 | 0 | 0 | 0 | 0 |
5628 | 200712170605 | NaN | 2007-12-17 06:05:41 | 733393.253947 | 27.357 | 107.751 | 10.0 | 4.6 | Mb | CHINA | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5629 | 200712191134 | NaN | 2007-12-19 11:34:56 | 733395.482593 | -6.952 | 155.778 | 45.0 | 5.9 | Ms | SOLOMON ISLANDS | ... | 5163 | 6594 | 7553 | 339 | 0 | 0 | 0 | 0 | 0 | 0 |
5630 | 200712200306 | NaN | 2007-12-20 03:06:56 | 733396.129815 | -32.713 | -71.788 | 30.0 | 5.5 | Mw | CHILE | ... | 19533 | 20027 | 24957 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5631 | 200712200755 | NaN | 2007-12-20 07:55:15 | 733396.330035 | -39.011 | 178.291 | 20.0 | 6.6 | Mw | NEW ZEALAND | ... | 3338 | 535 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5632 | 200712200948 | NaN | 2007-12-20 09:48:29 | 733396.408669 | 39.417 | 33.212 | 10.0 | 5.7 | Mw | TURKEY | ... | 43608 | 19583 | 3460 | 2834 | 2510 | 0 | 0 | 0 | 0 | 0 |
5633 | 200712220711 | NaN | 2007-12-22 07:11:08 | 733398.299398 | -2.407 | 139.067 | 20.0 | 6.2 | Mw | INDONESIA | ... | 38829 | 32845 | 10353 | 1913 | 1413 | 584 | 600 | 53 | 0 | 0 |
5634 | 200712221226 | NaN | 2007-12-22 12:26:17 | 733398.518252 | 2.087 | 96.806 | 23.0 | 6.1 | Mw | INDONESIA | ... | 11849 | 27534 | 9038 | 12514 | 1135 | 0 | 0 | 0 | 0 | 0 |
5635 | 200712231256 | NaN | 2007-12-23 12:56:13 | 733399.539039 | -4.043 | 39.253 | 10.0 | 4.9 | Mw | KENYA | ... | 172 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5636 | 200712262347 | NaN | 2007-12-26 23:47:10 | 733402.991088 | 39.446 | 33.162 | 8.0 | 5.6 | Mw | TURKEY | ... | 37158 | 10131 | 6902 | 1935 | 1289 | 0 | 0 | 0 | 0 | 0 |
5637 rows × 51 columns
x = [1,2,3,4,5] #a homogenous python list
y = np.array(x) #this is a numpy array, created with a list as input
#To multiply every element in a list by 5, you have to loop over the list
for i in range(0,len(x)):
x[i] = x[i] * 7
print(x)
#To multiply every element in a numpy array by 5, just do this:
y = y * 7
print(y)
[7, 14, 21, 28, 35] [ 7 14 21 28 35]
There are a number of packages and methods for optimizing Python code. Here is a very basic progression.
#Looping over two lists of a million values and doing calculations on each element
import datetime
a = range(0,1000000)
b = range(0,1000000)
t1 = datetime.datetime.now()
for i in range(0,len(a)):
a[i]**2 + b[i]**2 + 2*a[i]*b[i]
t2 = datetime.datetime.now()
dt1 = ((t2-t1).microseconds/1e6)
print('Looping over list: %.2e seconds' % dt1)
Looping over list: 5.38e-01 seconds
#Vectorizing with numpy
a = np.arange(1e6)
b = np.arange(1e6)
t3 = datetime.datetime.now()
a**2 + b**2 + 2*a*b
t4 = datetime.datetime.now()
dt2 = ((t4-t3).microseconds/1e6)
print('Numpy array: %.2e seconds' % dt2)
print('Numpy is %.1fx faster than looping over a list' % (dt1/dt2))
Numpy array: 2.50e-02 seconds Numpy is 21.5x faster than looping over a list
#Running on multiple cores with numexpr
import numexpr as ne
t5 = datetime.datetime.now()
ne.evaluate("a**2 + b**2 + 2*a*b") #Turn simple numpy expression into a string
t6 = datetime.datetime.now()
dt3 = ((t6-t5).microseconds/1e6)
print('numexpr eval: %.2e seconds' % dt3)
print('Numexpr is %.1fx faster than numpy' % (dt2/dt3))
numexpr eval: 4.70e-03 seconds Numexpr is 5.3x faster than numpy
from mapio.shake import ShakeGrid #ShakeGrid is a class that holds one variable of a ShakeMap grid XML file
shakemap = ShakeGrid.load('grid.xml',adjust='res')
eventdict = shakemap.getEventDict()
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0,0,1.0,1.0])
p = plt.imshow(shakemap.getLayer('mmi').getData(),cmap='jet',vmin=1,vmax=10)
plt.xlabel('Columns')
plt.ylabel('Rows')
locstr = eventdict['event_description']
mag = eventdict['magnitude']
datestr = eventdict['event_timestamp'].strftime('%b %d, %Y %H:%M:%S')
th = plt.title('PGA for %s - %s M%.1f' % (locstr,datestr,mag))
ch=plt.colorbar(shrink=0.7)
#To save the above figure to a file, you can call the matplotlib savefig() function
#NB: This function MUST be in the same cell as the code generating the figure
#The formats supported by savefig() can vary, but generally speaking PNG, JPG, and PDF are well supported.
plt.savefig('shakemap.png')
#extra code cell to clean the PNG file made above - delete or don't execute this cell if you want to keep it
#Any system command can be run from IPython by preceding it with "!"
!rm shakemap.png
There are a couple of fundamental packages for dealing with vector (points, lines, polygons) data in Python. They are:
First, let's read in a shapefile of the boundaries of the US (from US Census: https://www.census.gov/geo/maps-data/data/cbf/cbf_nation.html)
import fiona
from shapely.geometry import shape, Point
recs = fiona.open('cb_2015_us_nation_20m.shp','r')
#turn a record into a shapely shape (MultiPolygon, in this case)
uspoly = shape(recs[0]['geometry'])
#This multipolygon has the continental US, plus lots of islands.
#Let's just grab the continental US (should be the one with most vertices)
npoints = []
for i in range(0,len(uspoly)):
poly = uspoly[i]
npoints.append(len(poly.exterior.coords))
npoints = np.array(npoints)
imax = npoints.argmax()
continent = uspoly[imax]
recs.close()
You can look at a shape in Jupyter simply by typing the name of the variable.
continent
Let's try looking at a list of recent (1 day) earthquakes and figuring out which ones are inside the continental US...
import urllib.request as request
import json
feedurl = 'http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_day.geojson'
fh = request.urlopen(feedurl)
data = fh.read().decode('utf-8')
jdict = json.loads(data)
insidelat = []
insidelon = []
outsidelat = []
outsidelon = []
for event in jdict['features']:
lon,lat,depth = event['geometry']['coordinates']
p = Point(lon,lat)
if continent.contains(p):
print('Event %s is inside continental US' % event['id'])
insidelat.append(lat)
insidelon.append(lon)
else:
outsidelat.append(lat)
outsidelon.append(lon)
fh.close()
Event mb80154014 is inside continental US Event us20006482 is inside continental US Event ci37385847 is inside continental US
We can also plot the inside/outside points on top of our continent...
cx,cy = zip(*continent.exterior.coords)
f = plt.figure(figsize=(12,8))
plt.plot(cx,cy,'b')
plt.plot(insidelon,insidelat,'r*')
plt.plot(outsidelon,outsidelat,'gd')
plt.axis([-140,-60,25,50])
lh=plt.legend(['Continental US','Earthquakes inside','Earthquakes outside'],numpoints=1)
ObsPy is a Python package which provides tools to access, process, and visualize seismological data.
Notes:
Streams are list-like objects which contain multiple Trace objects, i.e. gap-less continuous time series and related header/meta information.
import warnings
warnings.simplefilter("ignore", DeprecationWarning)
from obspy.clients.neic import Client
from obspy import UTCDateTime
from obspy import read
client = Client()
t = UTCDateTime() - 5 * 3600 # 5 hours before now
stream = client.get_waveforms("IU", "ANMO", "00", "BH?", t, t + 120) #get 2 minutes worth of data
fig = plt.figure(figsize=(12,4))
ph = stream.plot()
<matplotlib.figure.Figure at 0x11961e390>
#First we need to write a quick function to get the last 7 or 30 days's worth of focal mechanisms.
#To do this, we'll pull down the last 7 days of events from the web site, and extract the strike,
#dip and rake values for those events which have had a moment tensor generated for them.
import json
def getRecentFM():
feed7 = 'http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_week.geojson'
feed30 = 'http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_month.geojson'
fh = request.urlopen(feed7)
data = fh.read().decode('utf-8')
fh.close()
jdict = json.loads(data)
focals = []
for event in jdict['features']:
etypes = event['properties']['types'].strip(',').split(',')
if 'moment-tensor' in etypes or 'focal-mechanism' in etypes:
lon,lat,depth = event['geometry']['coordinates']
eurl = event['properties']['detail']
fh = request.urlopen(eurl)
edata = fh.read().decode('utf-8')
fh.close()
edict = json.loads(edata)
mt = edict['properties']['products']['moment-tensor'][0]
if 'nodal-plane-1-strike' not in mt['properties']:
continue
strike = float(mt['properties']['nodal-plane-1-strike'])
dip = float(mt['properties']['nodal-plane-1-dip'])
if 'nodal-plane-1-rake' in mt['properties']:
rake = float(mt['properties']['nodal-plane-1-rake'])
else:
rake = float(mt['properties']['nodal-plane-1-slip'])
focals.append((lat,lon,strike,dip,rake))
return focals
eventlist = getRecentFM()
len(eventlist)
12
Set up a map with Basemap, a mapping extension to the matplotlib package - handy for replacing GMT :)
from mpl_toolkits.basemap import Basemap
from obspy.imaging.beachball import beach
fig = plt.figure(figsize=(12,12))
ax = fig.add_axes([0,0,1.0,1.0])
m = Basemap(projection='cyl', lon_0=142.36929, lat_0=38.3215,
resolution='c',ax=ax)
m.drawcoastlines()
m.fillcontinents(lake_color='#2E2EFE')
m.drawparallels(np.arange(-90., 120., 30.))
m.drawmeridians(np.arange(0., 420., 60.))
m.drawmapboundary(fill_color='#2E2EFE')
for event in eventlist:
lat,lon = event[0:2]
if lon < 0:
lon += 360 #cylindrical projections need longitudes 0-360
fm = event[2:]
x, y = m(lon, lat)
b = beach(fm, xy=(x, y), width=10, linewidth=1, alpha=0.85,facecolor='r')
b.set_zorder(10)
bc=ax.add_collection(b)
Below are the first few lines of a Turkish strong-motion data set. The smtools package provides a function for reading that data.
!head 20111023104120_4902.txt
import warnings
warnings.filterwarnings("ignore")
from smtools.turkey import readturkey
tracelist = readturkey('20111023104120_4902.txt')
ewtrace = tracelist[0][1] #Get the EW component of the station in this data file from Turkey
net = ewtrace.stats['network']
station = ewtrace.stats['station']
location = ewtrace.stats['location']
channel = ewtrace.stats['channel']
channel_id = '%s.%s.%s.%s' % (net,station,location,channel)
print(channel_id)
ewtrace.plot()
TR.4902..EW
ewtrace.detrend('linear')
ewtrace.detrend('demean')
#ewtrace.taper(max_percentage=0.05, type='cosine')
ewtrace.plot()
ewtrace.filter('highpass',freq=0.02,zerophase=True,corners=4)
pga = np.max(np.abs(ewtrace.data))
print('Peak acceleration of %s is %.2f m/s^2' % (channel_id,pga))
vtrace = ewtrace.copy()
vtrace.integrate()
pgv = np.max(np.abs(vtrace.data))
print('Peak velocity of %s is %.2f m/s' % (channel_id,pgv))
vtrace.plot()
Peak acceleration of TR.4902..EW is 0.56 m/s^2 Peak velocity of TR.4902..EW is 0.11 m/s
In the Terminal window where you started IPython, hit Ctrl-C and answer 'y' when presented with this question:
Shutdown this notebook server (y/[n])?
with plt.xkcd():
x = np.arange(0,12)
y = np.random.rand(12,1)
plt.plot(x,y,'b')
plt.xlabel('Months')
plt.ylabel('Profits')