In GEOPHYS330, we use python to download, plot and do a regression on a hundred years of sea level measurements in the Auckland harbour:
# these are the modules that contain functions we will use in the script:
import csv # this is to read a "comma separate values" file
from urllib.request import urlopen # to read csv data from a website
import matplotlib.pyplot as plt # plotting functions
import numpy as np # numerical tools
from scipy.stats import linregress # linear regression function
import pandas as pd
# you can find and load the sea level data for Auckland here:
url="http://www.psmsl.org/data/obtaining/rlr.monthly.data/150.rlrdata"
df = pd.read_csv(url,delimiter=';')
# "bad" values are defined by a large negative number, so we only
# accept positive values:
df = df[df.iloc[:,1] >= 0]
# convert the pandas dataframe to numpy matrix:
npdat =df.as_matrix()
time= npdat[:,0]
height = npdat[:,1]
slope, intercept, r_value, p_value, std_err = linregress(time,height)
print("r-squared: %f", r_value**2)
# plot the sea height data:
plt.plot(time,height,'ro')
plt.plot(time,intercept+slope*time,'k')
plt.grid()
plt.xlabel('Date (year)')
plt.ylabel('Sea height (mm)')
plt.title('Sea height at Princes Wharf, Auckland (New Zealand)')
plt.legend(['measurements','best linear fit'],loc=4)
plt.axis('tight')
#plt.savefig('sealevel_auckland.png')
plt.show()
