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 <spurdl...@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 <giftedlife2...@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@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.
______________________________________________
R-help@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.
