R/Examples/ExWind.R

In Splinets: Functional Data Analysis using Splines and Orthogonal Spline Bases

#------------------------------------------------#
#----------- Plotting the Wind data -------------#
#------------------------------------------------#

data(wind) #activating the data

wind1=wind[,-1] #Removing YEAR as irrelevant

#Transforming data to daily with the periodic form, i.e. the arguments in [0,1], 
#which is required in the periodic case.

numbdays=length(wind1[,1])/24 

Days=vector(mode='list', length=numbdays)

for(i in 1:numbdays){
  Days[[i]]=wind1[i*(1:24),]
  Days[[i]][,c(4,6)]=Days[[i]][,c(4,6)]/360 #the direction in [0,1]
}

#Raw discretized data for the first day 

par(mfrow=c(2,2))
hist(Days[[1]][,4],xlim=c(0,1),xlab='Wind direction',main='First day 10[m]')
hist(Days[[1]][,6],xlim=c(0,1),xlab='Wind direction',main='First day 50[m]')

plot(Days[[1]][,4],Days[[1]][,5],xlim=c(0,1),pch='.',cex=4,xlab='Wind direction',ylab='Wind speed')
plot(Days[[1]][,6],Days[[1]][,7],xlim=c(0,1),pch='.',cex=4,xlab='Wind direction',ylab='Wind speed')

#First Day Data:
#Projections of the histograms to the periodic spline form 

FirstDayDataF1=cbind(hist(Days[[1]][,4],xlim=c(0,1),breaks=seq(0,1,by=0.1))$mids,
hist(Days[[1]][,4],xlim=c(0,1),breaks=seq(0,1,by=0.1))$counts)


k=3
N=2 
n_knots=2^N*k-1 #the number of internal knots for the dyadic case
xi = seq(0, 1, length.out = n_knots+2) 
#Note that the range of the argument is assumed to be between 0 and 1

PrF1=project(FirstDayDataF1,xi,periodic = TRUE, graph = TRUE)

F1=PrF1$sp #The first day projection of the direction histogram at 10[m]

#Projections of the scatterplots to the periodic spline form 
#The bivariate sampl
FirstDayDataF1V1=as.matrix(Days[[1]][,4:5]) #we note that wind directions are scaled but not ordered 

#Padding the data with zeros as the sampling frequency is not sufficiently dense over [0,1]
FirstDayDataF1V1=rbind(FirstDayDataF1V1,cbind(seq(0,1,by=1/24),rep(0,25)))

#Another knot selection with more knots but still dyadic case
k=4
N=3 
n_knots2=2^N*k-1 #the number of internal knots for the dyadic case
xi2 = seq(0, 1, length.out = n_knots2+2) 

#For illustration one can plot the B-splines and the corresponding splinet
so = splinet(xi2,smorder = k, periodic = TRUE,norm = TRUE) 

plot(so$bs)
plot(so$bs,type='dyadic') #To facilitate the comparison with the splinet better 
                          #one can choose the dyadic grapph
plot(so$os)


#Projecting direction/wind data onto splines
PrS1=project(FirstDayDataF1V1,xi2,smorder=k,periodic = TRUE, graph = TRUE)

S1=PrS1$sp

#the next 7 days 
days= 7
#Transforming to the periodic data 

#The direction histogram
for(i in 2:days){
  DataF1=cbind(hist(Days[[i]][,4],plot=FALSE,breaks=seq(0,1,by=0.1))$mids,
               hist(Days[[i]][,4],plot=FALSE,breaks=seq(0,1,by=0.1))$counts)
  PrF1=project(DataF1,xi,periodic = TRUE)
  
  F1=gather(F1,PrF1$sp) #Collecting projections of daily wind-direction histograms at 10[m]
  
}

plot(F1) #plot of all daily functional data wind direction distributions


#Wind direction vs speed data at 10[m]
for(i in 2:days){
  DataF1V1=as.matrix(Days[[i]][,4:5]) #we note that wind directions are scaled but not ordered 
  
  #Padding the data with zeros as the sampling frequency is not sufficiently dense over [0,1]
  DataF1V1=rbind(DataF1V1,cbind(seq(0,1,by=1/24),rep(0,25)))
  
  PrS1=project(DataF1V1,xi2,smorder=k,periodic = TRUE)
  
  S1=gather(S1,PrS1$sp) #Collecting projections of daily wind-direction histograms at 10[m]
  
}

plot(S1) #plot of all daily functional data wind speed at wind direction

#Computing means of the data
A=matrix(rep(1/days,days),ncol=days)

MeanF1=lincomb(F1,A)
plot(MeanF1)

MeanS1=lincomb(S1,A)
plot(MeanS1)