-
Notifications
You must be signed in to change notification settings - Fork 71
Example Pytides Usage
This short example will demonstrate how Pytides can be used for prediction. To predict the tide, we need some water level data. For example let's use data from NOAA about King's Point, New York - which is identified by the station id 8516945. The data can be accessed here. It is a fairly simple matter to scrape this information. For this example, I collected two years of readings (in GMT and metres) from 2000-01-01 to 2001-12-31. NOAA provide readings every 6-minutes, but I will use hourly readings for my analysis. This speeds up the process significantly without sacrificing much accuracy. I stored the data in a file called data/8516945 - the first line of which reads:
2000-01-01 00:00 1.887
Once we've fit our model, we want to use it to predict some tides. Let's predict the tide for the first week of 2013 - bear in mind that's over a decade after any of our original readings. For comparison, I scraped NOAA predictions and verified readings for this same period in a file (data/8516945_noaa). The first line of this file reads:
2013-01-01 00:00 -0.276 -0.081
where -0.276 is the actual height observed and -0.081 is the height that NOAA predicted. I will use matplotlib to plot my prediction, the NOAA prediction, and the actual tide for comparison. Now onto the code:
from datetime import datetime
from pytides.tide import Tide
import numpy as np
import matplotlib.pyplot as plt
##Prepare our tide data
station_id = '8516945'
heights = []
t = []
f = open('data/'+station_id, 'r')
for i, line in enumerate(f):
t.append(datetime.strptime(" ".join(line.split()[:2]), "%Y-%m-%d %H:%M"))
heights.append(float(line.split()[2]))
f.close()
#For a quicker decomposition, we'll only use hourly readings rather than 6-minutely readings.
heights = np.array(heights[::10])
t = np.array(t[::10])
##Prepare a list of datetimes, each 6 minutes apart, for a week.
prediction_t0 = datetime(2013,1,1)
hours = 0.1*np.arange(7 * 24 * 10)
times = Tide._times(prediction_t0, hours)
##Fit the tidal data to the harmonic model using Pytides
my_tide = Tide.decompose(heights, t)
##Predict the tides using the Pytides model.
my_prediction = my_tide.at(times)
##Prepare NOAA's results
noaa_verified = []
noaa_predicted = []
f = open('data/'+station_id+'_noaa', 'r')
for line in f:
noaa_verified.append(line.split()[2])
noaa_predicted.append(line.split()[3])
f.close()
##Plot the results
plt.plot(hours, my_prediction, label="Pytides")
plt.plot(hours, noaa_predicted, label="NOAA Prediction")
plt.plot(hours, noaa_verified, label="NOAA Verified")
plt.legend()
plt.title('Comparison of Pytides and NOAA predictions for Station: ' + str(station_id))
plt.xlabel('Hours since ' + str(prediction_t0) + '(GMT)')
plt.ylabel('Metres')
plt.show()
A few seconds later, and we have our result:
Not bad! The difference in the peaks between the Pytides prediction and the NOAA prediction could be due to the mean water level increasing in the last decade, or possibly the NOAA increase the amplitude of their high predictions as a safety precaution. Of course, neither prediction perfectly matches reality - partly because we cannot predict the meteorological conditions which affect the tide.