W4_inclass_loma_prieta_guerrero_anthony

W4_inclass_loma_prieta_guerrero_anthony October 2, 2023 1 Loma Prieta Earthquake, Earthquake Occurrence Statistics: Omori’s Law Last week we: - pandas DataFrames, indexing, and data cleaning. - Load marine geophysical data (bathymetry and marine magnetic anomalies) from two oceanic ridges. - Select data and drop rows with NaNs. - Plot bathymetry data and evaluate spreading rate. - Declare a functions to detrend data, calculate distance and estimate spreading rate - Plot marine magnetic anomaly data and compare spreading rates. Our goals for today: - Load a Bay Area seismic catalog of January 01,1989 to December 31, 1995. - Compute the distance and time interval between Loma Prieta quake and subsequant earthquakes to indentify aftershocks. - Filter the aftershocks from the catalog and look at their distribution. - Examine the Gutenberg-Richter statistic - Make and earthquake forecast 1.1 Setup Run this cell as it is to setup your environment. [1]: import numpy as np import pandas as pd import datetime import matplotlib import matplotlib.pyplot as plt import matplotlib.dates as mdates import cartopy.crs as ccrs import cartopy.feature as cfeature On October 17 at 5:15pm PDT (October 18 1989 at 04:15am UTC) the M6.9 Loma Prieta earth- quake occurred in the Santa Cruz mountains approximately 80 km southwest of the Berkeley Cam- pus. We will use this earthquake sequence to investigate the performance of catalog declustering algorithm. https://en.wikipedia.org/wiki/1989_Loma_Prieta_earthquake 1.2 Load the Earthquake Catalog Load the .csv data file of all the earthquakes from January 01,1989 to December 31, 1995 in the ANSS (Advanced National Seismic System) catalog from between latitudes 36.0-38.5° and longitude -123.0 to -121.0° ( http://ncedc.org/anss/catalog-search.html ). 1

[2]: # read data #pd.to_datetime? LP_catalog = pd . read_csv( 'data/bay_area_anss_1989_1995.csv' ) #the following function just reformats to a DateTime object format LP_catalog[ 'DateTime' ] = pd . to_datetime(LP_catalog[ 'DateTime' ]) LP_catalog . head() [2]: Unnamed: 0 DateTime Latitude Longitude Depth Magnitude \ 0 0 1989-01-01 03:51:37.170 36.5928 -121.1027 0.00 1.26 1 1 1989-01-01 11:39:34.870 36.9583 -121.5762 5.70 1.07 2 2 1989-01-01 21:05:18.030 37.5293 -121.6905 7.03 0.99 3 3 1989-01-01 21:54:27.470 37.5590 -121.6760 7.38 0.85 4 4 1989-01-01 21:57:12.240 36.6493 -121.2678 4.66 1.17 MagType NbStations Gap Distance RMS Source EventID 0 Md 4 150.0 6.0 0.13 NC 1160762.0 1 Md 34 36.0 1.0 0.08 NC 129467.0 2 Md 11 80.0 4.0 0.02 NC 129474.0 3 Md 12 66.0 1.0 0.04 NC 129475.0 4 Md 21 65.0 1.0 0.04 NC 129476.0 DateTime objects are very useful. There are functions to parse the DateTime data into individual arrays, or functions to process the DateTime data in different ways. [8]: LP_catalog . iloc[ 3 ][ 'DateTime' ] . minute [8]: 54 [10]: # create data arrays, it will speed up our loops later # Note .dt provides functional resources to work on the DateTime file. year = LP_catalog[ 'DateTime' ] . dt . year month = LP_catalog[ 'DateTime' ] . dt . month day = LP_catalog[ 'DateTime' ] . dt . day lat = LP_catalog[ 'Latitude' ] . values lon = LP_catalog[ 'Longitude' ] . values mag = LP_catalog[ 'Magnitude' ] . values nevt = len (year) #number of events mag [10]: array([1.26, 1.07, 0.99, …, 0.89, 1.06, 1.1 ]) 1.3 Map the Raw Earthquake Catalog On a map of the Bay Area plot the location of events in the raw catalog. Scale the marker color and size by magnitude. Code for you to write 2

[12]: #Make a Map of the earthquake catalog # Set Corners of Map lat0 =36.0 lat1 =38.5 lon0 =-123.0 lon1 =-121.0 tickstep =0.5 #for axes latticks = np . arange(lat0,lat1 + tickstep,tickstep) lonticks = np . arange(lon0,lon1 + tickstep,tickstep) # coordinates for UC Berkeley Berk_lat = 37.8716 Berk_lon = -122.2727 plt . figure( 1 ,( 10 , 10 )) ax = plt . axes(projection = ccrs . PlateCarree()) ax . set_extent([lon0, lon1, lat0, lat1], crs = ccrs . PlateCarree()) ax . coastlines(resolution = '10m' ,linewidth =1 ) ax . set_xticks(lonticks) ax . set_yticks(latticks) ax . set(xlabel = 'Longitude' , ylabel = 'Latitude' ,title = 'Raw Catalog' ) # Sort by magnitude to plot largest events on top LP_catalog_sorted = LP_catalog . sort_values( 'Magnitude' ) #exponent to scale marker size z = np . exp(LP_catalog_sorted[ 'Magnitude' ]) plt . scatter(LP_catalog_sorted[ 'Longitude' ], LP_catalog_sorted[ 'Latitude' ], s = z, c = LP_catalog_sorted[ 'Magnitude' ], cmap = 'plasma' ,alpha =0.4 ,marker = 'o' ) ␣ ↪ # plot circles on EQs plt . plot(Berk_lon,Berk_lat, 's' ,color = 'limegreen' ,markersize =8 ) # plot red ␣ ↪ square on Berkeley Campus plt . colorbar(label = 'Magnitude' ) plt . show() 3

d0 = datetime.date(year[0], month[0], day[0]) ## Plot Magnitude vs. Time for the Raw Catalog Plot magnitude vs. time for the raw catalog and print out the number of events as we did in-class. Use the alpha = 0.2 argument in plot to make the markers transparent. We have an array of magnitudes but what about time? We need to construct an array of fractions of years from the first entry in the catalog. Code for you to write [14]: d0 = datetime . date(year[ 0 ], month[ 0 ], day[ 0 ]) d0 4

[14]: datetime.date(1989, 1, 1) [21]: #Let's create a time array and then use the plot() function to plot the ␣ ↪ earthquakes #as points. d0 = datetime . date(year[ 0 ], month[ 0 ], day[ 0 ]) time = np . zeros(nevt) # initialize the size of the array days for i in range ( 0 ,nevt, 1 ): d1 = datetime . date(year[i], month[i], day[i]) time[i] = ((d1 - d0) . days) /365.25 + year[ 0 ] #Now the plot fig, ax = plt . subplots(figsize = ( 6 , 6 )) ax . plot(time, mag, 'o' ,alpha =0.2 ,markersize =5 ) ax . set(xlabel = 'Time' , ylabel = 'Magnitude' , title = 'Raw Event Catalog' ) ax . grid() plt . show() 5