SoilAirSensors/Analysis.R
2023-08-19 23:38:38 +02:00

201 lines
6.9 KiB
R

library(gplots)
setwd("/home/pi/Desktop/SoilAirSensors/data/")
FlankLength = 30
Steps = 30
MinS = -1
MaxS = 1
CorCols = colorRampPalette(c('blue','grey','green'))(n=Steps)
CorBreaks = seq(from=MinS,to=MaxS,length.out=(Steps+1))
Date = paste(format(Sys.Date()-1,c("%Y","%-m","%-d")),collapse="-")
Data = read.table(paste("filter/soilair_",Date,".dat",sep=""),head=T)
pdf(file=paste("pdf/soilair_",Date,".pdf",sep=""))
MeanData = data.frame(
DailySecond = Data$DailySecond,
GPIO_C = rep(NA,nrow(Data)),
GPIO_PY = rep(NA,nrow(Data)),
Humidity = rep(NA,nrow(Data)),
Temperature = rep(NA,nrow(Data)),
Pressure = rep(NA,nrow(Data))
)
MedData = data.frame(
DailySecond = Data$DailySecond,
GPIO_C = rep(NA,nrow(Data)),
GPIO_PY = rep(NA,nrow(Data)),
Humidity = rep(NA,nrow(Data)),
Temperature = rep(NA,nrow(Data)),
Pressure = rep(NA,nrow(Data))
)
for(idx in seq(FlankLength,(nrow(Data)-FlankLength),FlankLength)){
MeanData$GPIO_C[idx] = mean(Data$GPIO_C[ (idx-FlankLength+1) : (idx+FlankLength)])
MeanData$GPIO_PY[idx] = mean(Data$GPIO_PY[ (idx-FlankLength+1) : (idx+FlankLength)])
MeanData$Humidity[idx] = mean(Data$Humidity[ (idx-FlankLength+1) : (idx+FlankLength)])
MeanData$Temperature[idx] = mean(Data$Temperature[ (idx-FlankLength+1) : (idx+FlankLength)])
MeanData$Pressure[idx] = mean(Data$Pressure[ (idx-FlankLength+1) : (idx+FlankLength)])
MedData$GPIO_C[idx] = median(Data$GPIO_C[ (idx-FlankLength+1) : (idx+FlankLength)])
MedData$GPIO_PY[idx] = median(Data$GPIO_PY[ (idx-FlankLength+1) : (idx+FlankLength)])
MedData$Humidity[idx] = median(Data$Humidity[ (idx-FlankLength+1) : (idx+FlankLength)])
MedData$Temperature[idx] = median(Data$Temperature[ (idx-FlankLength+1) : (idx+FlankLength)])
MedData$Pressure[idx] = median(Data$Pressure[ (idx-FlankLength+1) : (idx+FlankLength)])
}
MeanData = na.omit(MeanData)
rownames(MeanData) = 1:nrow(MeanData)
MedData = na.omit(MedData)
rownames(MedData) = 1:nrow(MedData)
write.table(MeanData,file=paste("filter/soilair-mean_",Date,".dat",sep=""),col.names = T,row.names = F)
write.table(MedData,file=paste("filter/soilair-median_",Date,".dat",sep=""),col.names = T,row.names = F)
heatmap.2(
main="Pearson Correlation - Original",
cor(Data,method="pearson"),
Colv=F,
Rowv=F,
dendrogram="none",
breaks=CorBreaks,
col=CorCols
)
heatmap.2(
main=paste("Pearson Correlation - Mean of",FlankLength,"+1"),
cor(MeanData,method="pearson"),
Colv=F,
Rowv=F,
dendrogram="none",
breaks=CorBreaks,
col=CorCols
)
heatmap.2(
main=paste("Pearson Correlation - Median of",FlankLength,"+1"),
cor(MeanData,method="pearson"),
Colv=F,
Rowv=F,
dendrogram="none",
breaks=CorBreaks,
col=CorCols
)
###### PLOT DAILY SECOND; VIEW DATA ######
plot(GPIO_C~DailySecond,Data,main=paste("GPIO_C",Date),xlab="Daily Second [s]",ylab="Frequency [Hz]",pch=19,cex=0.2)
points(GPIO_C~DailySecond,MeanData,col="red",pch=19,cex=0.2)
points(GPIO_C~DailySecond,MedData,col="green",pch=19,cex=0.2)
abline(h=mean(Data$GPIO_C),col="red")
abline(h=median(Data$GPIO_C),col="green")
plot(GPIO_PY~DailySecond,Data,main=paste("GPIO_PY",Date),xlab="Daily Second [s]",ylab="Frequency [Hz]",pch=19,cex=0.2)
points(GPIO_PY~DailySecond,MeanData,col="red",pch=19,cex=0.2)
points(GPIO_PY~DailySecond,MedData,col="green",pch=19,cex=0.2)
abline(h=mean(Data$GPIO_PY),col="red")
abline(h=median(Data$GPIO_PY),col="green")
plot(Humidity~DailySecond,Data,pch=19,cex=0.2,main=Date)
points(Humidity~DailySecond,MeanData,pch=19,cex=0.2,col="red")
points(Humidity~DailySecond,MedData,pch=19,cex=0.2,col="green")
plot(Temperature~DailySecond,Data,pch=19,cex=0.2,main=Date)
points(Temperature~DailySecond,MeanData,pch=19,cex=0.2,col="red")
points(Temperature~DailySecond,MedData,pch=19,cex=0.2,col="green")
plot(Pressure~DailySecond,Data,pch=19,cex=0.2,main=Date)
points(Pressure~DailySecond,MeanData,pch=19,cex=0.2,col="red")
points(Pressure~DailySecond,MedData,pch=19,cex=0.2,col="green")
###### PLOT ENVIRONMENT AGAINST ENVIRONMENT; VIEW DATA ######
plot(Humidity~Temperature,Data,pch=19,cex=0.2,main=Date)
points(Humidity~Temperature,MeanData,pch=19,cex=0.2,col="red")
points(Humidity~Temperature,MedData,pch=19,cex=0.2,col="green")
plot(Temperature~Humidity,Data,pch=19,cex=0.2,main=Date)
points(Temperature~Humidity,MeanData,pch=19,cex=0.2,col="red")
points(Temperature~Humidity,MedData,pch=19,cex=0.2,col="green")
plot(Humidity~Pressure,Data,pch=19,cex=0.2,main=Date)
points(Humidity~Pressure,MeanData,pch=19,cex=0.2,col="red")
points(Humidity~Pressure,MedData,pch=19,cex=0.2,col="green")
plot(Pressure~Humidity,Data,pch=19,cex=0.2,main=Date)
points(Pressure~Humidity,MeanData,pch=19,cex=0.2,col="red")
points(Pressure~Humidity,MedData,pch=19,cex=0.2,col="green")
plot(Pressure~Temperature,Data,pch=19,cex=0.2,main=Date)
points(Pressure~Temperature,MeanData,pch=19,cex=0.2,col="red")
points(Pressure~Temperature,MedData,pch=19,cex=0.2,col="green")
plot(Temperature~Pressure,Data,pch=19,cex=0.2,main=Date)
points(Temperature~Pressure,MeanData,pch=19,cex=0.2,col="red")
points(Temperature~Pressure,MedData,pch=19,cex=0.2,col="green")
###### PLOT GPIO AGAINST ENVIRONMENT; VIEW DATA ######
plot(GPIO_C~Humidity,Data,pch=19,cex=0.2,main=Date)
points(GPIO_C~Humidity,MeanData,pch=19,cex=0.2,col="red")
points(GPIO_C~Humidity,MedData,pch=19,cex=0.2,col="green")
plot(GPIO_PY~Humidity,Data,pch=19,cex=0.2,main=Date)
points(GPIO_PY~Humidity,MeanData,pch=19,cex=0.2,col="red")
points(GPIO_PY~Humidity,MedData,pch=19,cex=0.2,col="green")
plot(GPIO_C~Temperature,Data,pch=19,cex=0.2,main=Date)
points(GPIO_C~Temperature,MeanData,pch=19,cex=0.2,col="red")
points(GPIO_C~Temperature,MedData,pch=19,cex=0.2,col="green")
plot(GPIO_PY~Temperature,Data,pch=19,cex=0.2,main=Date)
points(GPIO_PY~Temperature,MeanData,pch=19,cex=0.2,col="red")
points(GPIO_PY~Temperature,MedData,pch=19,cex=0.2,col="green")
plot(GPIO_C~Pressure,Data,pch=19,cex=0.2,main=Date)
points(GPIO_C~Pressure,MeanData,pch=19,cex=0.2,col="red")
points(GPIO_C~Pressure,MedData,pch=19,cex=0.2,col="green")
plot(GPIO_PY~Pressure,Data,pch=19,cex=0.2,main=Date)
points(GPIO_PY~Pressure,MeanData,pch=19,cex=0.2,col="red")
points(GPIO_PY~Pressure,MedData,pch=19,cex=0.2,col="green")
###### PLOT GPIO AGAINST GPIO; VIEW DATA ######
plot(GPIO_PY~GPIO_C,Data,pch=19,cex=0.2,main=Date)
plot(GPIO_PY~GPIO_C,MedData,pch=19,cex=0.2,main=paste("Median",Date))
abline(lm(GPIO_PY~GPIO_C,MedData),col="red")
###### Missing Data Analysis ######
MissingTimePoint = c()
for(i in 2:nrow(Data)){
if(Data$DailySecond[i] - Data$DailySecond[i-1] != 1){
MissingTimePoint[length(MissingTimePoint)+1] = i
}
}
plot(1:length(MissingTimePoint),MissingTimePoint,pch=19,cex=0.2)
MaxDiff = 0
for( i in 2:length(MissingTimePoint)){
MaxDiff = max(MaxDiff, MissingTimePoint[i] - MissingTimePoint[i-1] )
}
MisDist = rep(0,length(MissingTimePoint)-1)
for( i in 2:length(MissingTimePoint)){
MisDist[i-1] = MissingTimePoint[i] - MissingTimePoint[i-1]
}
hist(MisDist)
plot(1:length(MisDist),MisDist,type="b",pch=19,cex=0.2)
dev.off()