R help - Feb 2020 - Identification of Turning Points in a Signal

Dear Abby,
Many thanks for your feedback on this.

True, the structure of my question was not clear.  I realized it later
and have been thinking of a better way to re-posting. I am sorry about
that. I will make extra efforts here to make it clearer.

I have just a single data. I have tried to use dput to send the full
data (see attached, Ogbos-dput). This is the raw data with labels
year, month, day and cosmic ray count.

I have two codes for the same data. The first is a combination of
Fourier transform and an R code for pits/trough identification. I am
displaying the script here:
data <- read.table("OULU05", col.names = c("year",
"month", "day", "counts"))

new.century <- data$year < 50

data$year <- ifelse(new.century, data$year + 2000, data$year + 1900)

data$date <- as.Date(ISOdate(data$year, data$month, data$day))
x = data$date
 y = data$counts
 y1<-approx(x,y,xout=x)$y
 RECON = 12
 f = fft(d<-y1)
p<-1:length(f)
 f[(p > RECON/2 + 1) & (p <= length(d) - RECON/2)] = 0
 f<-fft(f, inverse = TRUE)/length(f)
 data$smooth = abs(f)
 data$residual = y1 - data$smooth

 library(pastecs)
 tp<-turnpoints(data$residual)
 res<-(1:length(data$residual))[tp$pit]
 minima<-which(tp$pit & data$residual<= -100)
 dates<-data$date[minima]
 k<-data$residual[minima]

png("2005A.png")

 plot(x, data$residual/100, type = "l")
 points(dates, k/100, lty = 2, col = 4)
dev.off()

The above runs OK as you may see.

The second code is similar but without the Fourier transform part. I
have only standardized the data in the second script. My aim here is
the same, to pick turning points/pits in the standardized (raw) data.
The script I tried to use is:
data <- read.table("OULU05", col.names = c("year",
"month", "day", "counts"))

new.century <- data$year < 50

data$year <- ifelse(new.century, data$year + 2000, data$year + 1900)

data$date <- as.Date(ISOdate(data$year, data$month, data$day))
x = data$date
 y = data$counts
y<-c(y-mean(y))/mean(y)*100
  data$residual = y

 library(pastecs)
 tp<-turnpoints(data$residual)
 res<-(1:length(data$residual))[tp$pit]
 minima<-which(tp$pit & data$residual<= -100)
 dates<-data$date[minima]
 k<-data$residual[minima]

png("2005B.png")

 plot(x, data$residual/100, type = "l")
 points(dates, k/100, lty = 2, col = 4)
dev.off()

This second script is where I have issues. It runs with warning/error.
And this is where I need help.


Thanks again for your kind assistance.

Best regards
Ogbos

On Sun, Feb 16, 2020 at 3:31 AM Abby Spurdle <spurdle.a at gmail.com>
wrote:
>
> Sorry no one replied sooner.
>
> Note that I find your question difficult to follow.
>
> It sounds like you have two datasets, each with different sizes.
> (e.g. Two data.frame objects, each with a different numbers of rows).
>
> Given that tp$pits is a logical index, trying to apply it to a vector
> (or data.frame) of a different size is problematic.
> Assuming that you're dealing with datasets of different sizes, the
> simplest solution is to modify your code, such that they're the same
> size.
>
>
> On Mon, Feb 10, 2020 at 6:40 PM Ogbos Okike <giftedlife2014 at
gmail.com> wrote:
> >
> > Dear Friends,
> > Wishing you the best of the day.
> >
> > I have a data (Cosmic Ray) which exhibit flow patterns of a
> > sine/cosine wave, i.e. decreasing/increasing and registering crests
> > (points maximal increases) and troughs/pits (points maximal
> > decreases). These turning points are of interest to me. With pastecs
> > package and a few lines of code as (the residual is coming from
> > Fourier transformation of the data):
> > library(pastecs)
> >  tp<-turnpoints(data$residual)
> >  res<-(1:length(data$residual))[tp$pit]
> >  minima<-which(tp$pit & data$residual<= -100)
> >  dates<-data$date[minima]
> >  k<-data$residual[minima]
> > I usually pick all the turning points (trough) equal or below -100. If
> > I change the <= to >=, I pick all the crests.
> >
> > Now, without first transforming the data, I wish to pick the same
> > turning points in the raw data. Indeed, the difference between the
> > transformed data and the raw data lies only in the amplitude of the
> > crest or trough, otherwise, the crests and trough are the same in both
> > signals.
> >
> > When I tried the above code in the raw signal, the warning/error
message is:
> > Warning message:
> > In tp$pit & data$residual <= -100 :
> >   longer object length is not a multiple of shorter object length.
> >
> > A sample of the raw data is:
> > 03 10 01 6.20636953199224
> > 03 10 02 6.90829266565563
> > 03 10 03 6.40434785174345
> > 03 10 04 6.33235573547028
> > 03 10 05 5.99039318317273
> > 03 10 06 5.09049172975812
> > 03 10 07 4.35257253795814
> > 03 10 08 4.49655677050448
> > 03 10 09 4.49655677050448
> > 03 10 10 4.4425626832996
> > 03 10 11 5.16248384603129
> > 03 10 12 5.72042274714835
> > 03 10 13 6.26036361919711
> > 03 10 14 5.8284109215581
> > 03 10 15 5.30646807857763
> > 03 10 16 5.32446610764592
> > 03 10 17 5.68442668901176
> > 03 10 18 6.33235573547028
> > 03 10 19 6.80030449124588
> > 03 10 20 7.26825324702148
> > 03 10 21 6.83630054938246
> > 03 10 22 2.53477160206063
> > 03 10 23 2.55276963112892
> > 03 10 24 2.39078736951429
> > 03 10 25 -0.48889728141246
> > 03 10 26 -0.110938670978323
> > 03 10 27 0.303015997592397
> > 03 10 28 1.81485043932894
> > 03 10 29 -8.04806949009518
> > 03 10 30 -16.1471825708267
> > 03 10 31 -17.0470840242413
> > 03 11 01 -13.6094604721975
> > 03 11 02 -8.98396700164638
> > 03 11 03 -6.28426264140255
> > 03 11 04 -5.78031782749036
> > 03 11 05 -3.72854251370505
> > 03 11 06 -2.95462726376849
> > 03 11 07 -4.52045579270991
> > 03 11 08 -3.54856222302213
> > 03 11 09 -0.884853920914888
> > 03 11 10 0.447000230138735
> > 03 11 11 0.0150475324997218
> > 03 11 12 -0.308916990729538
> > 03 11 13 0.0690416197045984
> > 03 11 14 -0.110938670978323
> > 03 11 15 -0.938848008119764
> > 03 11 16 -3.02661938004166
> > 03 11 17 -3.92652083345627
> > 03 11 18 -3.24259572886117
> > 03 11 19 -1.67676719991974
> > 03 11 20 -2.30669821730997
> > 03 11 21 -2.9366292347002
> > 03 11 22 -2.75664894401728
> > 03 11 23 -3.44057404861238
> > 03 11 24 -4.34047550202699
> > 03 11 25 -3.87252674625139
> > 03 11 26 -2.72065288588069
> > 03 11 27 -2.25270413010509
> > 03 11 28 -1.37080070575878
> > 03 11 29 -0.0389465547051547
> > 03 11 30 0.033045561568014
> > the first three columns are year, month and day, the last column % CR
variation.
> >
> > abline (h=0) specifies values below the average. I am interested in
> > picking the time and magnitude of all the turning points below zero.
> >
> > Thank you for assisting me.
> > Best regards
> > Ogbos
> >
> > ______________________________________________
> > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.

Note that your post does not contain a minimal reproducible example.
I, and presumably most other readers, do not have the file "OULU05".

Also, your first post referred to "% CR variation", however, your
second post referred to "counts".

I created a simple simulated data set:

--------
sim.data = function ()
{   year = 1:100
    month = (0:99 %% 12) + 1
    day = (0:99 %% 28) + 1
    counts = sample (1:2000, 100)
    data.frame (year, month, day, counts)
}
data = sim.data ()
--------

After replacing the "data" object (as above), everything worked fine.
(Except that the inequality needed modification based on the value of counts).

Maybe the problem is with your dataset...?
Or maybe there's some step in your code that results in missing
values, given your input?

Also note that the head and tail functions, are useful for both
inspecting data, and describing your data to others.

--------
head (data)
tail (data)
--------

Dear Abby,
Thank you.  I will look at your stimulated data and then run the code with
it.

But since I am dealing with real data and also have volumes of it,  I would
like to send my real data to you.

The OULU05 is attached with dput function.  It is labeled Ogbos_dput.  I
would be surprised if it was stripped off.  Then I will resend it.  If,  on
the other hand,  you don't want it as dput data,  then I will email the
large data through your private box.
Warmest regards
Ogbos

On Sun, Feb 16, 2020, 08:14 Abby Spurdle <spurdle.a at gmail.com> wrote:
> Note that your post does not contain a minimal reproducible example.
> I, and presumably most other readers, do not have the file
"OULU05".
>
> Also, your first post referred to "% CR variation", however, your
> second post referred to "counts".
>
> I created a simple simulated data set:
>
> --------
> sim.data = function ()
> {   year = 1:100
>     month = (0:99 %% 12) + 1
>     day = (0:99 %% 28) + 1
>     counts = sample (1:2000, 100)
>     data.frame (year, month, day, counts)
> }
> data = sim.data ()
> --------
>
> After replacing the "data" object (as above), everything worked
fine.
> (Except that the inequality needed modification based on the value of
> counts).
>
> Maybe the problem is with your dataset...?
> Or maybe there's some step in your code that results in missing
> values, given your input?
>
> Also note that the head and tail functions, are useful for both
> inspecting data, and describing your data to others.
>
> --------
> head (data)
> tail (data)
> --------
>
	[[alternative HTML version deleted]]