<jupyter_start><jupyter_text>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  SetupRun this cell as it is to setup your environment.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>On October 17 at 5:15pm PDT (October 18 1989 at 04:15am UTC) the M6.9 Loma Prieta earthquake occurred in the Santa Cruz mountains approximately 80 km southwest of the Berkeley Campus. We will use this earthquake sequence to investigate the performance of catalog declustering algorithm.https://en.wikipedia.org/wiki/1989_Loma_Prieta_earthquake  Load the Earthquake CatalogLoad 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](http://ncedc.org/anss/catalog-search.html)).<jupyter_code># 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()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>#Lets explore some features of datetime
#  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<jupyter_output><empty_output><jupyter_text>Map the Raw Earthquake CatalogOn 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_**<jupyter_code>#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(...)
#exponent to scale marker size
z=np.exp(...)    

plt.scatter(..., ..., s=..., 
            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()<jupyter_output><empty_output><jupyter_text>Plot Magnitude vs. Time for the Raw CatalogPlot 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_**<jupyter_code>#Let's create a time array and then use the plot() function to plot the earthquakes
#as points.

d0 = datetime.date(..., ..., ...)
time=np.zeros(nevt) # initialize the size of the array days
for i in range(0,nevt,1):
    d1 = datetime.date(..., ..., ...)
    time[i]=((d1 - d0).days)/365.25 + year[0]


#Now the plot
fig, ax = plt.subplots(figsize=(6,6))
ax.plot(..., ...,'o',alpha=0.2,markersize=5)
ax.set(xlabel='Time', ylabel='Magnitude',
       title='Raw Event Catalog')
ax.grid()
plt.show()


#Here is another way to do it using some built in functions. It is as clear but it is 
#powerful using the datetime attributes, and you may want to use it when comparing other datetime ranges.

fig, ax = plt.subplots(figsize=(6,6))

ax.plot(mdates.date2num(LP_catalog['DateTime']), mag,'o',alpha=0.2,markersize=5)
locator = mdates.AutoDateLocator()
formatter = mdates.ConciseDateFormatter(locator)
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(formatter)
ax.set(xlabel='Time', ylabel='Magnitude',
       title='Raw Event Catalog')
ax.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Earthquake Catalog Statistics and AftershocksThere are two primary statistics that describe earthquake catalog data. The first is Omori's Law which describes the distribution of aftershocks following a primary earthquake or mainshock. It can be used to estimate the rate of aftershock production over time after the earthquake.The second is Gutenberg-Richter, which is used to characterize the rates of primary earthquakes over some period of time to determine rates of occurence, which can then be used to assess hazard and make earthquake forecasts.Both statistics ***require a means of identifying earthquakes as aftershocks***. Since aftershocks are related to the redistribution of mainshock stress including them in a statistic that assumes `primary (or mainshock)` event occurrence would lead to biased results. Aftershocks are therefore identified and removed from analysis. This is done based on spatial and temporal proximity of a given earthquake to another earlier, larger earthquake. There are several empirical models developed for this identification and the following cell describes one of these. We will start by examing the Loma Prieta aftershocks, which we need to isolate from the catalog.  Gardner and Knopoff (1974) Aftershock ModelGardner and Knopoff (1974) studied many hundreds of earthquake sequences and empirically determined distance and time filters as a function of mainshock magnitude to help identify aftershocks. The table below shows their model with distance to the event in km, and time after the event in days. Both distance from the mainshock and time afterward increase dramatically with increasing magnitude. Note that aftershocks following M5.5+ earthquakes can last for years.These windows are based on the following equations:To build our algorithm to identify aftershock using these windows we need to convert the year-month-day formate of dates to a timeline in number of days. We'll do this using the function `datetime.date` which for a given year, month, and day returns a datetime class object, which can be used to compute the time interval in days.We will also use the use the `np.power(base, exponent)` function to compute the distance and time interval bounds.   First, Determine the index value and the number of days from the beginning of the catalog to the Loma Prieta QuakeThe Loma Prieta earthquake is the largest in the catalog<jupyter_code>LP_ind = LP_catalog.index[...]
lp_mag=LP_catalog['Magnitude'][LP_ind].values
print(f'Loma Prieta Earthquake M{lp_mag[0]:.1f} Index={LP_ind[0]:d}')
#Some information about using datetime<jupyter_output><empty_output><jupyter_text>To calculate the days since the Loma Prieta earthquake for many entries we can use a loop**_Code for you to write_**<jupyter_code>#Determine the number of days from the Loma Prieta Earthquake
days=np.zeros(nevt) # initialize the size of the array days

for i in range(0,nevt,1):<jupyter_output><empty_output><jupyter_text>Why are some of the days[...] values negative?***Write your answer here***  Second, Define a function to compute the distance between to geographic points on a sphere - we will use the haversine function from last weekGreat-circle distance shown in red between two points on a sphere, P and Q. Source: https://en.wikipedia.org/wiki/Great-circle_distance<jupyter_code>#This function computes the spherical earth distance between to geographic points and is used in the
#declustering algorithm below
def haversine_np(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)

    All args must be of equal length.
    
    The first pair can be singular and the second an array

    """
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2]) # convert degrees lat, lon to radians

    dlon = lon2 - lon1
    dlat = lat2 - lat1

    a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2  # great circle inside sqrt

    c = 2 * np.arcsin(np.sqrt(a))   # great circle angular separation
    km = 6371.0 * c   # great circle distance in km, earth radius = 6371.0 km
    return km<jupyter_output><empty_output><jupyter_text>Next we need to define the code for applying the Gardner & Knopoff (1974) Aftershock Model**_Code for you to write_**<jupyter_code>#code here

print(f'Magnitude={mag[i]:.1f}  Dtest={Dtest:0.2f} km  Ttest={Ttest:.2f} days')<jupyter_output><empty_output><jupyter_text>Put the pieces together to build a aftershock detection algorithm**_Code for us to examine and write_**Note that you need to be careful about array dimensionality (column vs row ordered). The dimensionalty when performing certain (most) operations on arrays needs to be the same for both. We will also use np.array() which has as default column major ordering.<jupyter_code>#Declustering Algorithm Catalog
cnt=LP_ind[0] # initialize a counting variable at the LP index
save=np.zeros((1000000),dtype=int) # initialize a storage array

i=LP_ind[0]
# logical if statements to incorporate definitions of Dtest and Ttest aftershock window bounds


a=days[i+1:nevt]-days[i]    # time interval in days to subsequent earthquakes in catalog
m=mag[i+1:nevt]   # magnitudes of subsequent earthquakes in catalog
b=haversine_np(lon[i],lat[i],lon[i+1:nevt],lat[i+1:nevt]) # distance in km to subsequent EQs in catalog

#First Test find number of events that are within the aftershock time period
icnt=np.count_nonzero(a <= Ttest)   # counts the number of potential aftershocks
    
if icnt > 0:  # if there are potential aftershocks
    itime=np.array(np.nonzero(a <= Ttest)).transpose() + (i+1) # indices of potential aftershocks <= Ttest bound
    for j in range(0,icnt,1):   # loops over the aftershocks         
        if b[j] <= Dtest and m[j] < mag[i]: # test if the event is inside the distance window 
                                            # and that the event is smaller than the current main EQ
            save[cnt]=itime[j]  # index value of the aftershock
            cnt += 1 # increment the counting variable

                
af_ind=np.delete(np.unique(save),0)   # This is an array of indexes of aftershocks<jupyter_output><empty_output><jupyter_text>Create a set of arrays for the aftershock catalog. Use `af_ind` to select the aftershock events from the raw calalog.<jupyter_code># select the aftershock events
aftershock_df = LP_catalog.iloc[af_ind] 
aftershock_days=days[af_ind]  #The aftershocks are selected from the days array 
aftershock_mag=mag[af_ind]    #The aftershocks are selected from the mag array 
aftershock_lon=lon[af_ind]    #The aftershocks are selected from the lon array 
aftershock_lat=lat[af_ind]    #The aftershocks are selected from the lat array 
n2=len(aftershock_days)<jupyter_output><empty_output><jupyter_text>**_Concept question:_** How many aftershock events were there? **Write your answer here.**  Plot Magnitude vs. Time for the Aftershock CatalogPlot magnitude vs. time for the aftershock events and print out the number of events. Use the `alpha = 0.2` argument in `plot` to make the markers transparent. Plot the mainshock in a different color.<jupyter_code>#Plot DeClustered Catalog
fig, ax = plt.subplots(figsize=(6,6))

ax.plot(mdates.date2num(LP_catalog.iloc[LP_ind]['DateTime']), mag[LP_ind[0]],'o',color='red',markersize=10,label='Loma Prieta')
ax.plot(mdates.date2num(aftershock_df['DateTime']), aftershock_mag,'o',alpha=0.2,markersize=5,label='Aftershocks')
locator = mdates.AutoDateLocator()
formatter = mdates.ConciseDateFormatter(locator)
ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(formatter)

ax.set(xlabel='Time', ylabel='Magnitude',
       title='Aftershock Event Catalog')
ax.grid()
ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Plot a histogram of `aftershock_days`<jupyter_code>plt.hist(...)
plt.xlabel('Days after main shock')
plt.ylabel('Number of Events')
plt.title('Aftershocks of the Loma Prieta 1989 Quake')
plt.show()<jupyter_output><empty_output><jupyter_text>**_Discussion question:_** How would you describe the distribution of number of aftershocks with time after the main quake?By how much is the rate of aftershock occurrence reduced after 7 days? After 30 days? **Write your answers here.**  Map the Aftershock EventsOn a map of the Bay Area plot the location of events in the aftershock catalog. Again scale the marker color and size by magnitude, set `vmax=6.9` so the colorscale matches our original map. Add the locations of UC Berkeley campus and the Loma Prieta event epicenter.<jupyter_code>#Make a Map of Aftershock events

#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)

plt.figure(1,(10,10))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([lon0, lon1, lat0, lat1], crs=ccrs.PlateCarree())
ax.set_aspect('auto')
ax.coastlines(resolution='10m',linewidth=1) #downloaded 10m, 50m
ax.set_xticks(lonticks)
ax.set_yticks(latticks, crs=ccrs.PlateCarree())
ax.set(xlabel='Longitude', ylabel='Latitude',
       title='Aftershock Catalog')


#Sort Descending to plot largest events on top
indx=np.argsort(aftershock_mag)   #determine sort index
x=aftershock_lon[indx]            #apply sort index
y=aftershock_lat[indx]
z=np.exp(aftershock_mag[indx])    #exponent to scale size

plt.scatter(x, y, s=z, c=aftershock_mag[indx],vmax=6.9, cmap='plasma',alpha=0.4,marker='o',label='aftershocks')
plt.plot(lon[LP_ind[0]], lat[LP_ind[0]],'*',color='yellow',markersize=20,markeredgecolor='black',label='Loma Prieta Quake')
plt.plot(-122.2727,37.8716,'rs',markersize=8,label='UC Berkeley')
plt.legend()
plt.colorbar(label='Magnitude')
plt.show()<jupyter_output><empty_output><jupyter_text>Map of Bay Area faults. Source: https://pubs.er.usgs.gov/publication/fs20163020 **_Discussion questions:_**  - Were aftershocks only triggered on the same fault as the main quake?- What faults were active?  **Write your answers here.**  Lets decluster the whole catalog and examine the Gutenberg-Richter StatisticGutenberg devised a statistic that counts, for a given time period, the number of events with magnitude equal to or greater than a given magnitude. It is computed by defining a range of test magnitudes (sometimes called bins) and then determining the number of events of that magnitude and larger. Therefore it is a statistic on the rate of earthquake above a given magnitude. When plotted in log10(number(mag)) vs. mag there is a linear trend with an approximate -1 slope. That means that for each unit increase in magnitude there is a 10-fold decrease in the number of earthquakes. Gutenberg and Richter fit a line to the data compiled in this way to determine the intercept and slope (b-value, or rate of earthquake occurrence. $log10(N(m)) = A + B*M$<jupyter_code>#Write Algorithm to Decluster the entire catalog
#Reset variables
year=np.array(LP_catalog['DateTime'].dt.year)
month=np.array(LP_catalog['DateTime'].dt.month)
day=np.array(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<jupyter_output><empty_output><jupyter_text>For the Loma Prieta earthquake we only investigated the catalog from the time of the event to identify the aftershocks to plot the Omori curve.More typically the declustering algorithm is used to remove the aftershocks from a catalog. To apply the declustering algorithm to the entire catalog we need to pass through the entire catalog from beginning to end. To do this we need to add a loop.<jupyter_code>#We need to add a loop over time from what we had before

cnt=0
save=np.zeros(10000000,dtype=int) # initialize a counting variable

for i in range(nevt):                        #You only need this outside loop with proper indentation
    Dtest=np.power(10,0.1238*mag[i]+0.983)   # distance bounds
    if mag[i] >= 6.5:
        Ttest=np.power(10,0.032*mag[i]+2.7389)  # aftershock time bounds for M >= 6.5
    else:
        Ttest=np.power(10,0.5409*mag[i]-0.547)  # aftershock time bounds for M < 6.5
    
    a=days[i+1:nevt]-days[i]    # time interval in days to subsequent earthquakes in catalog
    m=mag[i+1:nevt]   # magnitudes of subsequent earthquakes in catalog
    b=haversine_np(lon[i],lat[i],lon[i+1:nevt],lat[i+1:nevt]) # distance in km to subsequent EQs in catalog
    
    icnt=np.count_nonzero(a <= Ttest)   # counts the number of potential aftershocks, 
                                    # the number of intervals <= Ttest bound
    
    if icnt > 0:  # if there are potential aftershocks
        itime=np.array(np.nonzero(a <= Ttest)).transpose() + (i+1) # indices of potential aftershocks <= Ttest bound
    
    for j in range(0,icnt,1):   # loops over the aftershocks         
        if b[j] <= Dtest and m[j] < mag[i]: # test if the event is inside the distance window 
                                            # and that the event is smaller than the current main EQ
            save[cnt]=itime[j]  # index value of the aftershock
            cnt += 1 # increment the counting variable
            
af_ind=np.delete(np.unique(save),0)   # This is an array of indexes of aftershocks<jupyter_output><empty_output><jupyter_text>The np.delete() function is used to remove aftershocks from the lat, lon, etc arrays. It is good to rename the receiving array (plat, plot, etc - p for 'primary').<jupyter_code>#This time Delete the aftershocks from the catalog and create new arrays
plat=np.delete(lat,af_ind)
plon=np.delete(lon,af_ind)
pmag=np.delete(mag,af_ind)
pyear=np.delete(year,af_ind)
pmonth=np.delete(month,af_ind)
pday=np.delete(day,af_ind)<jupyter_output><empty_output><jupyter_text>How many events were removed from the catalog? What percentage of the events are considered primary, or mainshocks?***Write your answers here***<jupyter_code># RePlot Mag vs Time and compare with what was obtained before

d0 = datetime.date(pyear[0], 1, 1)
time=np.zeros(len(pday)) # initialize the size of the array days
for i in range(0,len(pday),1):
    d1 = datetime.date(pyear[i], pmonth[i], pday[i])
    time[i]=((d1 - d0).days)/365.25 + year[0]


#This loop created a time array in units of years.
fig, ax = plt.subplots(figsize=(6,6))
ax.plot(time, pmag,'o',alpha=0.2,markersize=5)
ax.set(xlabel='Time', ylabel='Magnitude',
       title='Raw Event Catalog')
ax.grid()
plt.show()<jupyter_output><empty_output><jupyter_text>Map a map of the declustered earthquake catalog and compare with the complete catalog<jupyter_code>#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(by=['Magnitude'])
#exponent to scale marker size
z=np.exp(pmag)    

plt.scatter(plon, plat, s=z, 
            c=pmag, 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()<jupyter_output><empty_output><jupyter_text>Compute the Gutenberg Richter Statistic- Compile the annual number of earthquakes above a given magnitude- Examine the distribution- Fit the Gutenberg Richter model to the data- Use the Gutenberg Richter model to forecast recurrence of earthquakes$log10(N(m)) = A + B*M$Let's try it.**_Code for us to examine and write_**<jupyter_code># Compute G&R Statistics
min_mag=0.0
max_mag=max(pmag)

m=np.arange(min_mag,max_mag,0.1)    #The magnitude "bins"
N=np.zeros(len(m))                  #N is the number of events
numyr=1995-1989 + 1                 #number of years in the catalog - note 7 years is too low for this analysis

for i in range(0,len(m),1):
    N[i]=np.log10(...)
 
#plot it
plt.figure()
plt.plot(m,N,'o')
plt.xlabel('Magnitude')
plt.ylabel('log10(N)')
plt.title('Gutenberg Ricter')
plt.show()
#Next Fit the G&R model to the data and plot

y=np.polyfit(...)
intercept=y[1]
slope=y[0]

print(f'The declustered data B value: {slope:0.3f}')

#plot it
plt.figure()
plt.plot(m,N,'o')
plt.plot(m,intercept+slope*m,'-',color='orange')
plt.xlabel('Magnitude')
plt.ylabel('log10(N)')
plt.title('Gutenberg Ricter')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09871685965931953, 'fraction_numerical': 0.021460186882460767, 'mean_word_length': 4.979149959903769, 'pattern_counts': {'":': 0, '<': 93, '<?xml version=': 0, '>': 85, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '482626', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7382, 'n_tokens_neox': 6885, 'n_words': 2449, 'path': 'week04_inclass/W4_inclass_loma_prieta.ipynb', 'repo_name': 'ds-connectors/EPS-88-FA21', 'size': '22367'}