{"id":1549,"date":"2013-05-19T15:32:13","date_gmt":"2013-05-19T19:32:13","guid":{"rendered":"http:\/\/lukemiller.org\/?p=1549"},"modified":"2016-09-06T17:31:45","modified_gmt":"2016-09-07T00:31:45","slug":"more-tide-prediction-with-r","status":"publish","type":"post","link":"https:\/\/lukemiller.org\/index.php\/2013\/05\/more-tide-prediction-with-r\/","title":{"rendered":"More tide prediction with R"},"content":{"rendered":"

Edit: <\/strong>There is now a full-fledged R package, rtide<\/a>, to accomplish the same basic task of generating time series of tide predictions that is outlined here. See this more recent post<\/a> for information.<\/p>\n

In the previous post<\/a> I outlined how to query the XTide<\/a> software with R and parse the results into a handy-dandy data frame. The biggest hurdle with that method is getting XTide up and running on your computer. The code outlined here works entirely within R, so you don’t need XTide installed on your computer. The trade-off is that it provides substantially less information. This script only produces tide height predictions at a user-specified interval (1 to 60 minutes), and doesn’t pick out high tide, low tide, sunrise, sunset, or any of the other fun stuff that XTide can produce. Additionally, this version only works for a subset of 637 of the NOAA tide stations rather than the 1000s of stations that XTide can output. The available stations are found on this page<\/a> listed with the type ‘Ref’, rather than the ones labeled ‘Sub’. If all you want is a series of continuous tide height predictions for a NOAA reference tide station, this script and the associated Rdata file will do it for you. Fair warning: the actual prediction routine is not vectorized, so if you ask for many weeks of high-frequency data, the script will take a few minutes to run.<\/p>\n

 <\/p>\n

\"An<\/a>
An illustration of tide height predictions from the script.<\/figcaption><\/figure>\n

The script tide_prediction.R<\/a><\/code> and the harmonics file Harmonics-20120302.Rdata<\/a><\/code> can be found in my GitHub repository<\/a>. The tide_prediction.R script is shown below.<\/p>\n

\n
\n
# Filename: tide_prediction.R<\/span>\r\n\r\n# NOT FOR NAVIGATIONAL USE. DON'T WRECK YOUR BOAT BASED ON MY BAD CODING.<\/span>\r\n# Generate tide height predictions for US NOAA *reference* tide stations<\/span>\r\n# entirely within R.<\/span>\r\n# This requires a previously-produced file of harmonics for the sites, stored as<\/span>\r\n# a Rdata file. The current file is Harmonics-20120302.Rdata. The harmonics file<\/span>\r\n# was originally derived from XTide, found at http:\/\/www.flaterco.com\/xtide\/<\/span>\r\n# The available reference tide stations are listed here:<\/span>\r\n# http:\/\/www.flaterco.com\/xtide\/locations.html<\/span>\r\n# Note that only stations listed as type 'Ref' on that page can be used<\/span>\r\n# with this script.<\/span>\r\n# The script assumes that all input times (and corresponding tide<\/span>\r\n# time outputs) are in UTC\/GMT time zone, not your local time zone. Check the<\/span>\r\n# output against NOAA's tide predictions at<\/span>\r\n# http:\/\/tidesandcurrents.noaa.gov\/tide_predictions.shtml<\/span>\r\n# This can produce predictions for years between 1700 and 2100, but with<\/span>\r\n# decreasing accuracy as you move further from 2012 (say beyond 1990 or 2030)<\/span>\r\n# since the harmonics file is derived from the current tidal epoch, and things<\/span>\r\n# like sea level have a tendency to change over time.<\/span>\r\n\r\n# This script may take a few minutes to run if you ask for weeks of<\/span>\r\n# high-frequency tide predictions.<\/span>\r\n# The output is stored in a data frame called 'results'.<\/span>\r\n\r\n# Author: Luke Miller May 19, 2013<\/span>\r\n###############################################################################<\/span>\r\n# Sys.setenv(TZ='UTC') # Change R session time zone to UTC (optional)<\/span>\r\n\r\n# User-friendly version for loading harmonics Rdata file:<\/span>\r\ncat<\/span><\/a>(<\/span>'Please load the harmonics Rdata file.\\n<\/span>'<\/span>)<\/span>\r\nload<\/span><\/a>(<\/span>file.choose<\/span><\/a>(<\/span>)<\/span>)<\/span>\r\n\r\n# Hard-coded version, if your harmonics file is always in the same place:<\/span>\r\n#load('D:\/R\/Weather_data\/Harmonics_20120302.Rdata')<\/span>\r\n\r\n# Specify tide station, Site names: http:\/\/www.flaterco.com\/xtide\/locations.html<\/span>\r\n# Only the stations labeled 'Ref' will work with this script.<\/span>\r\ncat<\/span><\/a>(<\/span>'Please enter the station name: \\n<\/span>'<\/span>)<\/span>\r\nstationID = scan<\/span><\/a>(<\/span>file<\/span><\/a> = ''<\/span>,<\/span> character<\/span><\/a>(<\/span>)<\/span>,<\/span> quiet = TRUE<\/span>,<\/span> n = 1<\/span>,<\/span> sep = '#'<\/span>)<\/span>\r\n# Hard-coded version:<\/span>\r\n#stationID = 'Monterey Harbor'<\/span>\r\n#stationID = 'Boston, Boston Harbor'<\/span>\r\n#stationID = 'San Diego, San Diego Bay'<\/span>\r\n#stationID = 'Boston'<\/span>\r\n\r\n# Find row index of the desired station in the harms list that was loaded from<\/span>\r\n# the harmonics Rdata file.<\/span>\r\nstInd = grep<\/span><\/a>(<\/span>stationID,<\/span> harms$station)<\/span>\r\n\r\n# If no stations are found, or multiple possible stations are found, notify the<\/span>\r\n# user.<\/span>\r\nif<\/span> (<\/span>length<\/span><\/a>(<\/span>stInd)<\/span> == 0<\/span>)<\/span>{<\/span>\r\ncat<\/span><\/a>(<\/span>'No station found.\\n<\/span>'<\/span>)<\/span>\r\n}<\/span> else<\/span> if<\/span> (<\/span>length<\/span><\/a>(<\/span>stInd)<\/span> > 1<\/span>)<\/span> {<\/span>\r\ncat<\/span><\/a>(<\/span>'Multiple stations found, please choose one (enter number): \\n<\/span>'<\/span>)<\/span>\r\nfor<\/span> (<\/span>i in<\/span> 1<\/span>:length<\/span><\/a>(<\/span>stInd)<\/span>)<\/span> {<\/span>\r\ncat<\/span><\/a>(<\/span>'['<\/span>,<\/span> i,<\/span> ']'<\/span>,<\/span> harms$station[<\/span>stInd[<\/span>i]<\/span>]<\/span>,<\/span>'\\n<\/span>'<\/span>)<\/span>\r\n}<\/span>\r\nchoice = scan<\/span><\/a>(<\/span>file<\/span><\/a> = ''<\/span>,<\/span> what = numeric<\/span><\/a>(<\/span>)<\/span>,<\/span> n = 1<\/span>,<\/span> quiet = TRUE<\/span>)<\/span>\r\nstInd = stInd[<\/span>choice]<\/span>\r\n}<\/span>\r\n\r\n# Show the name of the station being used<\/span>\r\ncat<\/span><\/a>(<\/span>'Using station: '<\/span>,<\/span> harms$station[<\/span>stInd]<\/span>,<\/span>'\\n<\/span>'<\/span>)<\/span>\r\n\r\n# Specify number of minutes between tide predictions<\/span>\r\ncat<\/span><\/a>(<\/span>'Enter desired frequency of predictions in minutes (1-60): \\n<\/span>'<\/span>)<\/span>\r\nfreq = scan<\/span><\/a>(<\/span>file<\/span><\/a> = ''<\/span>,<\/span> numeric<\/span><\/a>(<\/span>)<\/span>,<\/span> quiet = TRUE<\/span>,<\/span> n = 1<\/span>)<\/span>\r\n#freq = 6 # units of minutes, hard coded version<\/span>\r\n\r\nfreq = freq * 60<\/span> # convert minutes to seconds<\/span>\r\n\r\n# Specify start and end times for the tide predictions<\/span>\r\n# If you cross a year boundary, the predictions must use the<\/span>\r\n# correct equilibrium and nodefactors for each year.<\/span>\r\ncat<\/span><\/a>(<\/span>'Enter starting time (YYYY-MM-DD HH:MM) in UTC\/GMT time zone: \\n<\/span>'<\/span>)<\/span>\r\nt1 = scan<\/span><\/a>(<\/span>file<\/span><\/a> = ''<\/span>,<\/span> character<\/span><\/a>(<\/span>)<\/span>,<\/span> quiet = TRUE<\/span>,<\/span> n = 1<\/span>)<\/span>\r\ncat<\/span><\/a>(<\/span>'Enter ending time (YYYY-MM-DD HH:MM) in UTC\/GMT time zone: \\n<\/span>'<\/span>)<\/span>\r\nt2 = scan<\/span><\/a>(<\/span>file<\/span><\/a> = ''<\/span>,<\/span> character<\/span><\/a>(<\/span>)<\/span>,<\/span> quiet = TRUE<\/span>,<\/span> n = 1<\/span>)<\/span>\r\n# Hard coded version:<\/span>\r\n#t1 = '2012-12-31 23:00'<\/span>\r\n#t2 = '2013-01-01 09:00'<\/span>\r\n\r\n# Convert to POSIXct time format, with UTC time zone<\/span>\r\nt1 = as.POSIXct<\/span><\/a>(<\/span>t1,<\/span> tz = 'UTC'<\/span>)<\/span>\r\nt2 = as.POSIXct<\/span><\/a>(<\/span>t2,<\/span> tz = 'UTC'<\/span>)<\/span>\r\n\r\n# Create data frame with sequence of times at desired frequency<\/span>\r\ntimes = data.frame<\/span><\/a>(<\/span>POStime = seq<\/span><\/a>(<\/span>t1,<\/span>t2,<\/span>freq)<\/span>)<\/span>\r\n# Extract the year for each of those time points<\/span>\r\ntimes$yr = as.POSIXlt<\/span><\/a>(<\/span>times$POStime)<\/span>$year + 1900<\/span>\r\n# Use the year to make a character value representing the start of the year<\/span>\r\ntimes$yrstart = paste<\/span><\/a>(<\/span>times$yr,<\/span> '01'<\/span>,<\/span> '01'<\/span>,<\/span> sep = '-'<\/span>)<\/span>\r\n# Calculate hours since the start of the year for each time point<\/span>\r\ntimes$hours = as.numeric<\/span><\/a>(<\/span>difftime<\/span><\/a>(<\/span>times$POStime,<\/span>\r\nas.POSIXct<\/span><\/a>(<\/span>times$yrstart,<\/span> tz = 'UTC'<\/span>)<\/span>,<\/span> units<\/span><\/a> = 'hours'<\/span>)<\/span>)<\/span>\r\n# Make a year index to pull correct starting values from harms for the year<\/span>\r\ntimes$yrindx = times$yr - harms$startyear + 1<\/span>\r\n\r\n########################################<\/span>\r\n# Tide height calculation - Magic happens here<\/span>\r\n# Input times behave as if they are given in UTC, so output tides correspond to<\/span>\r\n# the UTC times.<\/span>\r\npred = vector<\/span><\/a>(<\/span>'numeric'<\/span>,<\/span>nrow<\/span><\/a>(<\/span>times)<\/span>)<\/span>\r\npred[<\/span>1<\/span>:nrow<\/span><\/a>(<\/span>times)<\/span>]<\/span> = harms$datum[<\/span>stInd]<\/span>\r\nfor<\/span> (<\/span>j in<\/span> 1<\/span>:nrow<\/span><\/a>(<\/span>times)<\/span>)<\/span> {<\/span>\r\nfor<\/span> (<\/span>i in<\/span> 1<\/span>:length<\/span><\/a>(<\/span>harms$name)<\/span>)<\/span> {<\/span>\r\npred[<\/span>j]<\/span> = pred[<\/span>j]<\/span> + (<\/span>(<\/span>harms$nodefactor[<\/span>i,<\/span>times$yrindx[<\/span>j]<\/span>]<\/span> *\r\nt<\/span><\/a>(<\/span>harms$A[<\/span>stInd,<\/span>i]<\/span>)<\/span> *\r\n(<\/span>cos<\/span><\/a>(<\/span>(<\/span>harms$speed[<\/span>i]<\/span> *\r\ntimes$hours[<\/span>j]<\/span> +\r\nharms$equilarg[<\/span>i,<\/span> times$yrindx[<\/span>j]<\/span>]<\/span> -\r\nt<\/span><\/a>(<\/span>harms$kappa[<\/span>stInd,<\/span>i]<\/span>)<\/span>)<\/span> *\r\npi\/180<\/span>)<\/span>)<\/span>)<\/span>)<\/span>\t\r\n}<\/span>\r\n}<\/span>\r\n# Finished tide height calculations, units of feet (sorry, that's the way it is)<\/span>\r\n####################################<\/span>\r\n\r\ntimesout = as.POSIXct<\/span><\/a>(<\/span>times$hours * 3600<\/span>,<\/span> origin = times$yrstart,<\/span> tz = 'UTC'<\/span>)<\/span>\r\n\r\n# Make an output data frame<\/span>\r\nresults = data.frame<\/span><\/a>(<\/span>TimeUTC = timesout,<\/span> TideHt.ft = pred)<\/span>\r\n\r\n# This will produce a column of times in the current R session's time zone,<\/span>\r\n# which isn't necessarily the time zone of your tide station.<\/span>\r\nresults$TimeLocal = c<\/span><\/a>(<\/span>results$TimeUTC)<\/span>\r\n\r\ncat<\/span><\/a>(<\/span>'Finished. See \"results\" data frame for output.\\n<\/span>'<\/span>)<\/span>\r\n# Optionally print the results<\/span>\r\nplot<\/span><\/a>(<\/span>results$TimeUTC,<\/span>results$TideHt.ft,<\/span> type = 'l'<\/span>,<\/span> col<\/span><\/a> = 'blue'<\/span>,<\/span> las = 1<\/span>,<\/span>\r\nxlab = 'Time, UTC time zone'<\/span>,<\/span> ylab = 'Tide Height, ft'<\/span>,<\/span>\r\nmain = harms$station[<\/span>stInd]<\/span>)<\/span><\/pre>\n<\/div>\n<\/div>\n
In addition to these tide prediction routines, I have also previously provided R scripts to query the NOAA servers to retrieve measured, verified historical tide height data. That process is shown in this older post.<\/a><\/p>\n","protected":false},"excerpt":{"rendered":"
Edit: There is now a full-fledged R package, rtide, to accomplish the same basic task of generating time series of tide predictions that is outlined here. See this more recent post for information. In the previous post I outlined how to query the XTide software with R and parse the results into a handy-dandy data […]<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[218,240],"tags":[58,62,157,153],"class_list":["post-1549","post","type-post","status-publish","format-standard","hentry","category-r-project","category-tide-prediction","tag-r-project","tag-tide-height","tag-tide-prediction","tag-xtide"],"_links":{"self":[{"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/posts\/1549","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/comments?post=1549"}],"version-history":[{"count":5,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/posts\/1549\/revisions"}],"predecessor-version":[{"id":2317,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/posts\/1549\/revisions\/2317"}],"wp:attachment":[{"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/media?parent=1549"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/categories?post=1549"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/tags?post=1549"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}