dankelley / oce

R package for oceanographic processing
http://dankelley.github.io/oce/
GNU General Public License v3.0
143 stars 42 forks source link

request : vignette on tidal analysis #1513

Open clayton33 opened 5 years ago

clayton33 commented 5 years ago

Based on some investigative work done with tidal analysis using python with in R outlined in issue #1344 . I'll start this next week, an will probably edit this first comment.

richardsc commented 5 years ago

You and @dankelley and your "I'll just edit the first comment". HOW WILL I KNOW WHEN YOU'VE ACTUALLY STARTED WORKING ON IT IF I GET NO NOTIFICATIONS?!?!?!

(edited to add, I'm not actually angry at you guys. The shouting was just for emphasis)

stale[bot] commented 5 years ago

This issue has been automatically marked as stale because 8 weeks have passed with no activity. (Adding the 'pinned' label will prevent this issue from ever going stale.)

dankelley commented 5 years ago

FYI https://flaterco.com/files/xtide/Tidal_Analysis_and_Predictions_SP_NOS_CO-OPS_3.pdf might be helpful. Its chapter 5 deals with analysis of currents. (I have not read it ... I just ran across the reference whilst looking for something else.)

dankelley commented 4 years ago

I unpinned this, because no user has requested it in a long time, so I think it falls into the "would be nice" category, and not the "must have" category.

richardsc commented 4 years ago

Actually, it has come up a few times at BIO, where currently the standard way of doing TCA is to use Matlab. Basically all our moored processing is moving to R, but it is still a problem for tidal analysis of moored current meters. I know @clayton33 has a lot on her plate, but I still think this is more of a “required” addition than a “would be nice”.

clayton33 commented 4 years ago

I always see this, and I feel bad that I haven't gotten around to it. But I might have some time when I go to sea this September. Does oce have a built in data set that would be suitable for the vignette ?

dankelley commented 4 years ago

No, we don't have a suitable dataset built in to oce. I suppose there might be a dataset going back to Foreman's manual, though. I'll look into that.

dankelley commented 4 years ago

Oh, and don't feel bad. The idea of keeping issues "pinned" is not to make us feel bad, but instead just so we don't forget entirely. Some issues are a lot older than this!

dankelley commented 4 years ago

There are data in appendix 2 of [1], but the copy I have is a scanned version, so the only way to get those particular data would be to use a good OCR on the file, or to get two independent people to type them in by hand (with two people to catch errors).

Possibly T_TIDE has some companion data. A source for T_TIDE material is [2] but I have not checked it for a sample file of currents. (That would be easier than typing in Foremans results from [1]...)

I wonder if @richardsc has some thoughts on a test dataset, since (I think) he has used T_TIDE for currents.

References

  1. Foreman, M. G. G. “Manual for Tidal Currents Analysis and Prediction.” Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean Sciences, Patricia Bay, 1978. https://www.dropbox.com/s/k6ztw6bt9lo1e9b/Foreman%20-%201978%20-%20Manual%20for%20tidal%20currents%20analysis%20and%20prediction.pdf?dl=0
  2. T_TIDE website https://www.eoas.ubc.ca/~rich/#T_Tide
dankelley commented 4 years ago

Oooh, I did some searching and came up with https://www.dfo-mpo.gc.ca/science/data-donnees/tidal-marees/index-eng.html which seems to hold software and data files. For some reason, though, I am unable to download from the links. Maybe that's a temporary thing, though. I'll try again in an hour.

dankelley commented 4 years ago

I was able to download from https://www.dfo-mpo.gc.ca/science/documents/data-donnees/tidal-marees/tidpack.zip and I see therein not just Foreman's manuals, but also data files.

The relevant data files tide8.dat for the north-south component and tide9.dat for the east-west component, both (I think) in cm/s.

These are in a fixed format, e.g. top 6 lines of tide8.dat are

1  7077        8 272999999999999999999999999999999999999999999999999999999999999
2  7077        8 27299999  166  262  141  -20 -176  -26   28  128  288  133 -144
1  7077        9 272 -108  -14   60   38  -95  -84   41   83   91   84   30   14
2  7077        9 272   23  110  165   95   25 -107 -143   14  -39   -1  -21 -141
1  7077       10 272 -188 -148   50  124   69  -82 -196 -187  -85  -24   19   17
2  7077       10 272  -84  -95    0   72   48  -56 -115  -74  -45 -159  -73 -116

and the following code seems to read them well:

f <- "tide8.dat"
## The data format is described on page 5 (PDF page 12) of Foreman (1978).
##
## Foreman, M. G. G. “Manual for Tidal Currents Analysis and Prediction.”
## Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
## Sciences, Patricia Bay, 1978.
widths <- c(2, 5, 9, 2, 2, rep(5, 12))
names <- c("KOLI", "JSTN", "ID", "IM", "IY", paste("KARD", 1:12, sep=""))
d <- utils::read.fwf(f, widths=widths, col.names=names)
head(d)

since the output is

  KOLI JSTN ID IM IY KARD1 KARD2 KARD3 KARD4 KARD5 KARD6 KARD7 KARD8 KARD9 KARD10 KARD11 KARD12
1    1 7077  8  2 72 99999 99999 99999 99999 99999 99999 99999 99999 99999  99999  99999  99999
2    2 7077  8  2 72 99999   166   262   141   -20  -176   -26    28   128    288    133   -144
3    1 7077  9  2 72  -108   -14    60    38   -95   -84    41    83    91     84     30     14
4    2 7077  9  2 72    23   110   165    95    25  -107  -143    14   -39     -1    -21   -141
5    1 7077 10  2 72  -188  -148    50   124    69   -82  -196  -187   -85    -24     19     17
6    2 7077 10  2 72   -84   -95     0    72    48   -56  -115   -74   -45   -159    -73   -116  

If I save that, the result is only 2K for this one component of velocity. I think that's certainly small enough to be added to oce, and I see a good benefit in doing so, and not just for a vignette ... this could be good for my classes, too).

Therefore, I plan to continue this work, connecting the KOLI=1 and =2 lines to get a 24 hour signal, and then string them together and using the ID (day), IM (month)andIY` (year) to get a timestamp. My goal will be to make a single item, perhaps retrieved with

data(tidecurrendata)

for convenience in testing.

Since I don't see any real negative in proceeding with this, I won't ask for coauthor votes ... I can always just delete the file, if others don't want to burn up 4K for it.

dankelley commented 4 years ago

I've added a new dataset to commit 2600fb282d28e246a8a580348603cd5391d261d8 of the "develop" branch. To see the docs, type

?tidalCurrent

and to get a sample graph like below, type

example(tidalCurrent)
Screen Shot 2020-08-05 at 11 18 07 AM
dankelley commented 3 years ago

@clayton33 unless you object, I would like to close this issue. We are going to release to CRAN in a few days (#1795) and I want to clean up the issue list. This one has not been active in a while. Since closed issues are still available, I propose that you just reopen this whenever you feel you want to work on it and seek new input.

OK?

clayton33 commented 3 years ago

Sounds good.

richardsc commented 2 months ago

I'm reopening this, because:

a) it was never actually completed (just lapsed 😄), and b) I've been wanting this for some work, and have made some progress using the original example @clayton33 posted in #1344, but now using the tidalCurrent dataset that has been added to oce.

A vignette might not be the right place for this, as the system dependencies are a little tricky (e.g. python, utide package with dependencies, reticulate package, etc), but at the very least it could be documented on the wiki or the webpage. I will likely write it up first for my blog.

As a teaser, see the below example. Note that I have pre-configured a python "conda" environment called tides that contains the required python libraries:

library(oce)
#> Loading required package: gsw
library(reticulate)
use_condaenv("tides")

data(tidalCurrent)
utide <- import("utide")
pandas <- import("pandas")
np <- import("numpy")
mpl <- import("matplotlib")

t <- tidalCurrent$time
u <- np$array(tidalCurrent$u)
v <- np$array(tidalCurrent$v)

tpy <- pandas$to_datetime(as.numeric(t), unit = "s", utc = TRUE)
coef <- utide$solve(tpy, u, v,
    lat = 45, nodal = FALSE, trend = FALSE, method = "ols",
    conf_int = "linear", Rayleigh_min = 0.95
)
#> solve: matrix prep ... solution ... done.
tide <- utide$reconstruct(t = tpy, coef = coef)
#> prep/calcs ... done.

if (!interactive()) pdf("01_utide_example.pdf")

par(mfrow=c(2, 1))
oce.plot.ts(t, u)
lines(t, tide["u"], col=2)
oce.plot.ts(t, v)
lines(t, tide["v"], col=2)

if (!interactive()) dev.off()
#> quartz_off_screen 
#>                 2
image

Created on 2024-09-04 with reprex v2.1.1

dankelley commented 2 months ago

I like this and think it would be good for a vignette. Maybe that's just me, but I would never look for a wiki on R packages. I'd look in the docs, and then in the vignettes. But not in a wiki.

Note that you can always do

stuff

and it won't try to run it. Basically, it just pretty-prints the code. You make up the PNGs separately and include them -- like in a blog posting...

Readers would also like to know a bit about how you set up the conda.

I think this is a great addition.

Thanks!!

dankelley commented 2 months ago

Sorry, for the 'stuff' I meant, of course, to prefix it with the normal R material, but put eval=FALSE in there. It's hard to enter triple-backtick stuff on gh comments.

richardsc commented 2 months ago

First step, blog post!

https://www.clarkrichards.org/2024/09/05/tidal-current-analysis-utide/