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()
Comments are closed.