<jupyter_start><jupyter_text>Scientific Research Project 7 - Practical 2  Spatial and temporal analysis of the earthquake datasetIn the first week, we downloaded a CSV file of earthquake information from the USGS database and used Python to make some map plots of the dataset. In this notebook, we're going to use some more Python to perform some analysis on the spatial and temporal distribution of earthquakes in the San Francisco region.**** **Firstly, copy the notebook to your Google Drive using the COPY TO DRIVE button.**

Now we'll import the packages that we need again. This should be familiar from last time.<jupyter_code># import several helpful packages that we'll use
import numpy as np # linear algebra
import pandas as pd # data processing. We can read in csv files using pd.read_csv('path/to/csv-file')
import matplotlib.pyplot as plt # package for making figures
import matplotlib.cm as cm # for getting colourmaps<jupyter_output><empty_output><jupyter_text>Now we have to load our earthquake data into Google Colab. This is the exact same as we did last week using the USGS CSV file that you downloaded. To do this, click the small folder icon on the left hand bar and then click `Upload to Session Storage`. You can then navigate to the CSV file we downloaded and add it to Google Colab.**NOTE** - you will have to re-upload the data to Google Colab whenever you want to run this notebook. ![practical1_fig4.PNG](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAvAAAAFiCAYAAACQ4OyvAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAKrrSURBVHhe7f1p0GZFme8LGyfOBz92xPnw1pfutiK69ynDji3u6NcXDLvbXUf6lBOhYCsbQ6Wk28bGBruxBBopFHkoFRRElLalRJDBomRuBm0GmQq0rGKLDGo5MCmCQkmhVYUM+eY/h7Uyr8xcwz09636e/4/I4LnXynXllVfmqvyvvHPl/RJF<KEY>//<KEY>ZBlDAS+RA579XBQNAAOmG4wxmEZiJRqYM4N0cG2OSBQb4RNfH533ZAVSWrYZ+EuCQoqGDmR90aRiRdNS75pmv237xEIs8iMpJ9MGIVHsZF8QuBhtgj89YzUbCnWVItv1q5zAyhHHH8RtkPZZEZts/7QktsO8uT4pz0dtnfYNYx917VFfS0tfSEhjEvmWxKDUL+1xeV+FMW5rj7b7pYzPa33wwr10nzdh6i/jHcWgrr9WaWRZgPXxWAqDX2/FHu5emIfp9NPtrDyWzhBCIijgJXagrgZ5neLBqjTQWsxAJa5Hsvkz10pRYj7H1zaKoRzZfGnZ0cBuyhWCo5fAKQ/sqcDQJGKrRIvfpr1iQRL5kZRj7dV+2s9xzF3sTBwbxE7YVj1jNRMaYmxi6H3Xebbpz2HbyX4cnkvb094zVRsF5SJv0ica+nFiO8ybq0/j+WaxasoK8xtbYb0DH9v6gia99+P7KYlDRObfBkP+eGivtT1kjDr31frfwrr8kp/NmNhE5Uo79WcKeEII6QAFvEQMfgnNA1jzQJ25NhQdTkCE5yN74UDcRDZfWnY08Luya/HRoRxBqe6JUNLkjuVp8du0V38B7+0ZYZGIOF/3ZgFobJu8ts+U231xaO6LMcgb9csG8oIxjJP/jFhn4tfQjxPbYd4q3gHhsaStW9qvwY+UZlvG78hWnL+9LTL/Nhjs8fja+Fh7e9S2+7RzvuyWmBYw91nSruG/Nz6tVP9/yxgejbNXnf4AAAAASUVORK5CYII=) Once you have uploaded the file and you can see it in the left hand bar, we can then load it into Python again using the code in the cell below.<jupyter_code>df = pd.read_csv('San_Francisco_earthquakes_1900-2020.csv') # This needs to correspond to the name of the file you uploaded. You can change it to represent your own file<jupyter_output><empty_output><jupyter_text>Spatial analysisLast week, we made some plots of how earthquakes were distributed, but we didn't think about how to **quantitatively measure** that distribution. One way of doing this is to do a _cluster analysis_: we will try to identify groups of earthquakes which have a similar spatial location.  One of the most common types of cluster analysis is called `K-means` clustering: this analysis will split the dataset into a user-defined number of clusters, or $K$ clusters (that's why it's called K-means). K-means clustering is a type of unsupervised machine learning algorithm. You can read this tutorial for details about how K-means clustering works: https://stanford.edu/~cpiech/cs221/handouts/kmeans.htmlYou can also look at the Python user guide for K-means which provides a lot of information: https://scikit-learn.org/stable/modules/clustering.htmlk-meansThe good thing about this algorithm is that it's very simple and easy to understand. However, one of the downsides is that we have to decide the number of clusters, or groups, that we want to end up with in advance. The result will be dependent on the number of clusters, or value of $K$, that you choose. Here we start with clustering into 5 groups ($K = 5$) but you should experiment and find a value of $K$ that looks appropriate for your dataset. We will start by importing the K-means package:<jupyter_code># import K means clustering package from scikit-learn: https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html
from sklearn.cluster import KMeans<jupyter_output><empty_output><jupyter_text>Now that we've got our dataframe, let's run the K-means clustering algorithm. We'll define a variable $K=5$ which sets the number of clusters that we want.<jupyter_code># set K = number of clusters you want
K = 5<jupyter_output><empty_output><jupyter_text>The next step is to get the X and Y values for clustering from our data. In our case, the X value is the longitude and the Y value is the latitude. We then run the clustering on the data.<jupyter_code># get the longitude and latitude
data = df[['longitude', 'latitude']]
# run the K-means clustering on the data with K number of clusters
kmeans = KMeans(n_clusters=K)
clusters = kmeans.fit_predict(data)<jupyter_output><empty_output><jupyter_text>The `clusters` variable is now a list of the cluster number that was assigned to each earthquake. Python always starts counting at 0, so we have 5 clusters: cluster 0, cluster 1, cluster 2, cluster 3, and cluster 4:<jupyter_code>clusters<jupyter_output><empty_output><jupyter_text>We can now make a simple scatter plot of the earthquakes coloured by the cluster number.<jupyter_code># make a simple X and Y plot with the points coloured by the cluster number
# the cmap variable determintes the colour scale. You can change it to any of the named matplotlib colour scales to change the colour: https://matplotlib.org/tutorials/colors/colormaps.html
cmap = cm.viridis

# plot the points using a scatter plot. Some important parameters that you might want to change are:
# s = controls the size of the markers
# marker = controls the style of the markers. You can change it to 'D' for diamonds '.' for points, '+' for crosses, etc.
# edgecolors = controls the colour of the outline of the marker
cbar = plt.scatter(data['longitude'], data['latitude'], c=clusters, s=100, cmap=cmap, marker='o', edgecolors='black')

# add axis labels
plt.xlabel('Longitude ($^\circ$)')
plt.ylabel('Latitude ($^\circ$)')

# add a colour bar
plt.colorbar(cbar, label='Cluster number')<jupyter_output><empty_output><jupyter_text>We can add the cluster number back in to our original dataframe, and then use pandas to save it as a new CSV file:<jupyter_code>df['cluster_number'] = clusters
df.to_csv('San_Francisco_earthquakes_plus_clustering.csv', index=False)<jupyter_output><empty_output><jupyter_text>Testing if points are randomly distributedAnother test, which might be more useful, would be to check whether the earthquakes are randomly distributed, often called _Complete Spatial Randomness (CSR)_, or whether they are clustered into groups. To do that, we will install and import some libraries for analysing point patterns:<jupyter_code># download and install the point pattern library
!pip install libpysal pointpats
import libpysal as ps
from pointpats import PointPattern
import pointpats.quadrat_statistics as qs
# make a point pattern from the latitude and longitude and print the summary
pp = PointPattern(data)
pp.summary()<jupyter_output><empty_output><jupyter_text>To test whether points in the pattern are randomly distributed, we will use what's called a _Quadrat Statistic_. This comes from ecology, where a common way of counting the distribution of, say, beetles or plants, is by dividing an area up in to rectangles and counting the number of organisms in each rectangle. If the organisms are randomly distributed, then there should be an even number of them in each rectangle. You can see more explanation in this document: https://www.seas.upenn.edu/~ese502/NOTEBOOK/Part_I/3_Testing_Spatial_Randomness.pdfHere we will set up our quadrats to be rectangle shapes with 3 rows (`ny`)and 3 columns (`nx`). One disadvantage of this method is that it can be sensitive to the number of rectangles that you choose. You can test the sensitivity by changing `nx` and `ny` in the function below.<jupyter_code>q_r = qs.QStatistic(pp,shape= "rectangle",nx = 3, ny = 3)<jupyter_output><empty_output><jupyter_text>We can plot the quadrats, which will show you how many earthquakes are located in each one. Just from looking at the numbers, we can see that there is a very uneven distribution of earthquakes across the region.<jupyter_code>q_r.plot()<jupyter_output><empty_output><jupyter_text>But how can we say that statistically? The easiest way to test the CSR hypothesis is to use the Pearson $\chi^2$ goodness-of-fit test. This works by comparing the _expected_ number of points in each rectangle if the dataset was random, to the _actual_ number of points in each rectangle that we observe.You can read this notebook for a very nice explanation of how this works (focus on Section 3.1): https://joparga3.github.io/spatial_point_pattern/Our null hypothesis, that the points are randomly distributed, would occur if we have a $\chi^2$ statistic with a $p$-value greater than 0.05. If $p > 0.05$ it means we can't reject the null hypothesis at a 95% confidence level.Our alternative hypothesis, that the points are clustered, would occur if $p < 0.05$. So let's look at our results:<jupyter_code>q_r.chi2 
q_r.chi2_pvalue<jupyter_output><empty_output><jupyter_text>We can see that we have a $p$ value _much_ smaller than 0.05! This means we can reject the null hypothesis and state that the earthquakes are not randomly located.  Exercise 1: Spatial analysisIn this section, we've explored a couple of statistical techniques that you could use to analyse the spatial distribution of earthquakes. Try the following:1. Think about what the K-means clusters tell you. Does the machine learning algorithm relate to the tectonic setting around San Francisco?2. Re-run the K-means clustering with different numbers of clusters (change the value of $K$). How does this affect the results?3. Practice your Python by modifying the scatter plot of earthquakes coloured by cluster number. You can try changing the colourscale, the size of the points, the marker style, etc. Can you work out how to save this plot as an image? (HINT: see notebook from last week for example of how to do this).4. The quadrat method for testing whether points are randomly located can be dependent on the number and size of rectangles that are used. Try varying the number of rectangles, and look at the difference in the $\chi^2$ statistic and $p$ value. Does it have a significant impact on our hypothesis testing? If it does, this would be quite worrying, as it means our statistic wouldn't be very reliable.**BONUS QUESTION** for those who want to mess around with Python more:5. At the end of the K-means clustering, we saved the output to a new CSV file. In a new cell, try to read this to a dataframe using pandas. Then use the code from the first notebook to make a `cartopy` map plot of the earthquakes coloured by cluster number. (HINT - modify the last code block in Practical 1).<jupyter_code><jupyter_output><empty_output><jupyter_text>Temporal analysisNow we're going to move on and think about analysing the distribution of earthquakes through _time_. Time series analysis is essential for many, many things in geographical research: just a few examples of time series data could be precipitation or temperature records, sea level, ice sheet volume, river discharge, carbon emissions, air pollution levels, etc., etc. Alongside geographical research, time series are widely used in finance, economics, business, or pretty much any area where a variable changes through time.I would argue that dealing with time series is a fundamental and employable skill that you can learn from your degree!We're going to learn how to make some simple plots and analyses of the earthquake time series from our pandas dataframe. One of the first things we could do is to simply plot the earthquakes through time according to their magnitude. `pandas` has a very nice function called `datetime` that is specially designed for dealing with timeseries. Our first step is to convert the index (the row numbers) of the dataframe to a datetime:<jupyter_code># Convert the dataframe index to a pandas datetime
df.index = pd.to_datetime(df['time'])<jupyter_output><empty_output><jupyter_text>If you print out the dataframe, you can see that the row numbers in the very left hand column have been replaced with the date. This means that we can use the date to refer to each row.<jupyter_code>df<jupyter_output><empty_output><jupyter_text>We'll now make what's called a "stem plot" to show the earthquakes through time. This is a plot where there is a point for each earthquake, with a vertical line connecting it to the X axis. It makes it easier to see what year each earthquake occurred in.<jupyter_code># set up a figure
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(10,4))
# add the grid lines
ax.grid(color='0.8', linestyle='--', which='both')

# make a stem plot
markerline, stemlines, baseline = ax.stem(df.index, df.mag, linefmt='black', use_line_collection=True)

# adjust the markers and lines
markerline.set_markeredgecolor('k')
markerline.set_markerfacecolor('red')
stemlines.set_linewidth(0.5)
baseline.set_color('white')
markerline.set_markersize(10)

# add x and y labels
plt.xlabel('Year', fontsize=14)
plt.ylabel('Magnitude ($M$)', fontsize=14, labelpad=20)
# set the Y axis to start at magnitude 3
plt.ylim(3,)

# save as a figure
plt.savefig('San_Francisco_earthquakes_stem_plot.png')<jupyter_output><empty_output><jupyter_text>Number of earthquakes through timeAnother analysis that we might want to do is to analyse how the number of earthquakes has changed through time. To do that, we can count the number of earthquakes in a given year. Our `pandas datetime` object can help us with this: we can simply group times by year:<jupyter_code># get the years from the dataframe
y = df.index.year
# group the dataframe by the year, and then get the count
freq_df = df.groupby(y)['id'].count().reset_index()
plt.plot(freq_df['time'], freq_df['id'], c='blue')
plt.ylabel('Number of earthquakes')
plt.xlabel('Time')<jupyter_output><empty_output><jupyter_text>We could also make this a bar plot, rather than a line plot:<jupyter_code>plt.bar(freq_df['time'], freq_df['id'], color='red')
plt.ylabel('Number of earthquakes')
plt.xlabel('Time')<jupyter_output><empty_output><jupyter_text>Summary statistics through timeWe might also want to look at how the mean or median of a parameter changes through time. To show an example, here I'll plot the mean depth of earthquakes through time:<jupyter_code># get the years from the dataframe
y = df.index.year
# group the dataframe by the year, and then get the mean depth in each year
depth_df = df.groupby(y)['depth'].mean().reset_index()

# group the dataframe by the year, and then get the mean magnitude in each year.
mag_df = df.groupby(y)['mag'].mean().reset_index()
depth_df
# plot the points using a scatter plot
plt.scatter(depth_df['time'], depth_df['depth'], c='cyan', edgecolor='black', s=60)

# axis labels
plt.xlabel('Time')
plt.ylabel('Mean depth (km)')<jupyter_output><empty_output><jupyter_text>We can also colour this by earthquake magnitude:<jupyter_code># define the colour map - you can change this to be any colour scale from the matplotlib options
cmap = cm.Reds
# plot the points using a scatter plot
cbar = plt.scatter(depth_df['time'], depth_df['depth'], c=mag_df['mag'], cmap=cmap, edgecolors='black', s=100)

# make the colourbar
plt.colorbar(cbar, label='Mean magnitude')

# axis labels
plt.xlabel('Time')
plt.ylabel('Mean depth (km)')<jupyter_output><empty_output><jupyter_text>Exercise 2: Temporal analysisIn this section, we've explored some plotting techniques to analyse the temporal distribution of earthquakes. Try the following:1. Modify the plot of the number of earthquakes through time to make it look nicer. You could change the line/bar colours, change the font size, etc. Save it as a PNG file. * You can see the matplotlib documentation on [line plots](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.plot.html) or [bar plots](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.bar.html) for help. Remember to look at the examples at the bottom of these pages: they show you nice examples of different plotting styles with code.2. Try to plot some other summary statistics through time. For example, you could:  * Make a plot of the median depth of each year through time.  * Make a plot of the mean magnitude of each year through time.  * Make a plot of the maximum magnitude of each year through time.   * Make sure to save these plots as PNG files. You can make these same figures to explore your own dataset for your project.3. Think about what these data are telling you. Are there particular years in which there are more earthquakes? Can you identify which year had the highest magnitude earthquake? 4. From the plots you've made so far, choose two of the largest earthquakes. Do some research and find out some more information about these earthquakes. If they are large events, you should be able to find reports on them by Googling the year in which they occurred.* Were there any fatalities resulting from these earthquakes?* What was the economic damage?5. As well as magnitude, earthquake depth is important for seismic hazard. Shallow earthquakes are often the most damaging, as the seismic waves have less distance through the Earth's crust to travel before reaching infrastructure, and therefore dissipate less energy. Looking at the plots of earthquake depth, do you see any trend in earthquake depth through time? 6. Modify your plot of the mean magnitude through time so that the points are coloured by earthquake depth. From this plot, which year would you expect the most damage from earthquake events to occur in?<jupyter_code><jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06476484194294525, 'fraction_numerical': 0.014869479017512942, 'mean_word_length': 5.51794687724336, 'pattern_counts': {'":': 0, '<': 88, '<?xml version=': 0, '>': 88, 'https://': 8, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '325097', 'license': 'MIT', 'n_tokens_mistral': 5679, 'n_tokens_neox': 5310, 'n_words': 2380, 'path': 'Practical_2_spatial_temporal.ipynb', 'repo_name': 'fclubb/SciRes-Earthquakes', 'size': '75566'}