Skip to content
Snippets Groups Projects
Commit bc6c28a1 authored by Rikkers, Roxann1's avatar Rikkers, Roxann1
Browse files

removed empty folder

parent 154c6aec
Branches andor
No related tags found
No related merge requests found
######################################################
# ANALYSIS FOR BEC-3 IN 2022
# Author: Roxann Rikkers, WLR-ABG
# Date: 2022
# Transform UBW data records to a percentage activity per area per day per cow
######################################################
rm(list = ls())
## HPC locations #############################################################################
# HPC work directory
setwd("/lustre/nobackup/WUR/ABGC/rikke007/KB_BEC3/")
# HPC location data Ines
rawdata <- "/lustre/shared/BEC3_Rox/Data_Ines/"
# get list of filenames in directory with unfiltered dayfiles
flistall <- list.files(path = rawdata, full.names = FALSE)
# keep only dayfiles in list
flist <- flistall[!startsWith(flistall, "treatment") &
!startsWith(flistall, "diagnosis") &
!startsWith(flistall, "milk")]
# Create new directory to save results and overview files
dataloc <- "/lustre/nobackup/WUR/ABGC/rikke007/KB_BEC3/Data/"
### Onedrive locations ########################################################################
# working directory
#setwd("C:/Users/rikke007/OneDrive - WageningenUR/Documents/WLR/KB_BEC3/")
# practice files location
#rawdata <- "C:/Users/rikke007/OneDrive - WageningenUR/Documents/WLR/KB_BEC3/Data/Data_Ines_voorbeeldje/"
# get list of filenames in directory with unfiltered dayfiles
#flistall <- list.files(path = rawdata, full.names = FALSE)
# keep only dayfiles in list
#flist <- flistall[!startsWith(flistall, "treatment") &
# !startsWith(flistall, "diagnosis") &
# !startsWith(flistall, "milk")]
# Create new directory to save filtered dayfiles and overview files
#dataloc <- "C:/Users/rikke007/OneDrive - WageningenUR/Documents/WLR/KB_BEC3/"
# for HPC and Onedrive the same
dir.create(paste0(dataloc, "Data/Data_filtered"))
#### Packages #############################################################################
#library(dplyr)
#library(reshape2)
#library(ggplot2)
for(f in flist){ # f = flist[1]
### Loop to go trough all location files per day ###
print(f)
locdata <- read.table(paste0(rawdata, "/", f), sep = ",", header = TRUE) # import file
day <- locdata[1, 2] # get day as a variable to use later
### Filtered Dataset information ########################################################
#
# Each second of each day the X and Y coordinates of each cow
# Areas in the barn are known as feed, drink, lying cubicles (A, B, C), concentrate feed and
# no specific area (walking areas?)
#
#########################################################################################
# add barn layout
locdata$area <- as.numeric(locdata$area)
locdata$area_exp <- ifelse(locdata$area == 1, "Cubicle_A",
ifelse(locdata$area == 2, "Cubicle_B",
ifelse(locdata$area == 3, "Cubicle_C",
ifelse(locdata$area == 4, "Feeding_rack",
ifelse(locdata$area == 5, "Drinking_through",
ifelse(locdata$area == 6, "Concentrate_feeder",
ifelse(locdata$area == 0, "No_specific_area", NA)))))))
area_usage <- data.frame(cowid = NULL, date = NULL, concfeed = NULL, cubA = NULL, cubB = NULL, cubC = NULL,
drink = NULL, feed = NULL, nospec = NULL, cuball = NULL, totalrec = NULL, missrec = NULL)
for(id in unique(locdata$cowid)) { # id = 280 id = 419
### Loop to get each cow individually ###
tmp <- locdata[locdata$cowid == id, ] # create temporary dataframe for ith cow
# Filtering of Ines changed area numbers for missing values to 0 and thus no specific area
# thus, if xnew or ynew are missing then area should be NA
tmp$areanew <- ifelse(is.na(tmp$xnew), yes = NA, no = tmp$area)
# get number of missing records per cow per day
nummissing <- nrow(tmp[is.na(tmp$areanew), ])
# calculate percentages of areas usage (e.g. number of feed / number of records)
tmp <- tmp[!is.na(tmp$areanew), ] # df without the missing values
concfeed <- sum(tmp$areanew == 'Concentrate_feeder') / nrow(tmp) * 100
cubA <- sum(tmp$area_exp == 'Cubicle_A') / nrow(tmp) * 100
cubB <- sum(tmp$area_exp == 'Cubicle_B')/ nrow(tmp) * 100
cubC <- sum(tmp$area_exp == 'Cubicle_C') / nrow(tmp) * 100
drink <- sum(tmp$area_exp == 'Drinking_through') / nrow(tmp) * 100
feed <- sum(tmp$area_exp == 'Feeding_rack') / nrow(tmp) * 100
nospec <- sum(tmp$area_exp == 'No_specific_area') / nrow(tmp) * 100
cuball <- (cubA + cubB + cubC)
# create temporary data frame with all previously calcualted percentages
tmparea <- data.frame(data.frame(cowid = id, date = day, concfeed = concfeed, cubA = cubA, cubB = cubB,
cubC = cubC, drink = drink, feed = feed, nospec = nospec, cuball = cuball))
# add total number of records per day
tmparea$totalrec <- nrow(tmp)
tmparea$missrec <- nummissing
# row bind to new data frame for each cow
area_usage <- rbind(area_usage, tmparea)
} # end loop individual cows
rm(tmp) ; rm(id) ; rm(tmparea) ; rm(concfeed, cubA, cubB, cubC, drink, feed, nospec, cuball)
# merge data frames of each day with number of records and missing records
if (f == flist[1]){
farea_usage <- area_usage
} else {
farea_usage <- rbind(farea_usage, area_usage) # rowbind to new data frame for each day
}
write.csv(farea_usage, paste0(dataloc, "Data/Data_filtered/area_usage_per_cow.txt"), row.names = FALSE)
} # end loop of each dayfile
### get means and SD for all cows per day
avg <- aggregate(farea_usage[, 3:10], list(farea_usage$date), FUN = mean)
colnames(avg)[1] <- "date"
colnames(avg)[2:9] <- paste0(colnames(avg)[2:9], "_mean")
stds <- aggregate(farea_usage[, 3:10], list(farea_usage$date), FUN = sd)
colnames(stds)[1] <- "date"
colnames(stds)[2:9] <- paste0(colnames(stds)[2:9], "_sd")
# combine means and standard deviation to 1 data frame and save this data frame
df <- merge(x = avg, y = stds, by = "date")
write.csv(df, paste0(dataloc, "Data/Data_filtered/area_usage_mean_sd.txt"), row.names = FALSE)
######################################################
# ANALYSIS FOR BEC-3 IN 2022
# Author: Roxann Rikkers, WLR-ABG
# Date: 2022
# Combine % area usage file (script Area_usage.R) and transformed
# weather data (script Transform_weatherdata.R)
#
# additional a fix in calculations of activity per day, however this variable
# is not used in downstream analysis. Parameter is not calculated correclty, not
# corrected for missing data/gaps.
######################################################
rm(list = ls())
### Onedrive locations ########################################################################
# working directory
setwd("C:/Users/rikke007/OneDrive - WageningenUR/Documents/WLR/KB_BEC3/")
# import files
areausage <- read.csv("Data/HPC_Area_usage/area_usage_per_cow.txt")
activity <- read.csv("Data/HPC_Activity/activity.txt")
wdata <- read.csv("Data/Weather_data/daily_THI_01012021_24092022.txt")
## added dec 2021 till march 2022 to dataset
area_extra <- read.csv("Data/HPC_Area_usage/extra_dates/area_usage_per_cow.txt")
act_extra <- read.csv("Data/HPC_Activity/extra_dates/activity.txt")
# add extra dates to original datasets
area_all <- rbind(area_extra, areausage)
min(area_all$date)
max(area_all$date)
act_all <- rbind(act_extra, activity)
min(act_all$date)
max(act_all$date)
### merge area results & activity results
tmp <- merge(x = area_all, y = act_all, by = c("date", "cowid"))
### include weatherdata per day
wdata$date <- as.character(wdata$date)
tmp$date <- gsub(pattern = "-", replacement = "", x = tmp$date)
THIs <- wdata[wdata$date %in% tmp$date, ]
length(unique(tmp$date)) # same number of days, so all days accounted for
# combine THI per day and area+activity data
fdata <- merge(x = tmp, y = THIs, by = "date")
### percentage activity niet goed berekend, hier aanpassen
fdata$perc_active_10sec <- fdata$num_active_10sec / (fdata$num_active_10sec + fdata$num_inactive_10sec)
fdata$perc_active_1min <- fdata$num_active_1min / (fdata$num_active_1min + fdata$num_inactive_1min)
colnames(fdata) <- c("date", "cowid",
"ConcFeed", "CubA", "CubB", "CubC", "Drink", "Feed", "Nospec", "CubAll", "Area_avail_rec", "Area_miss_rec",
"Tot_timepoints",
"Miss_1min", "Num_active_1min", "Num_inactive_1min", "perc_active_1min",
"Miss_10sec", "Num_active_10sec", "Num_inactive_10sec", "perc_active_10sec", "dist_moved",
"meanT", "meanRH", "meanTHI" )
fdata <- fdata[, c(1, 2, 11, 12, 3:10, 13, 14:25)]
write.csv(fdata, paste0("Data/Final_data/area_activity_weather_data_BEC3_", Sys.Date(), ".txt"), row.names = FALSE)
hist(fdata$meanTHI)
This diff is collapsed.
######################################################
# ANALYSIS FOR BEC-3 IN 2022
# Author: Roxann Rikkers, WLR-ABG
# Date: 2022
# Create meanTHI variable from hourly weather data
######################################################
rm(list = ls())
setwd("C:/Users/rikke007/OneDrive - WageningenUR/Documents/WLR/KB_BEC3/Scripts/")
#### Packages #############################################################################
library(dplyr)
library(reshape2)
library(ggplot2)
dir.create("C:/Users/rikke007/OneDrive - WageningenUR/Documents/WLR/KB_BEC3/Data/Weather_data")
#### weather data #########################################################################
wdata <- read.table("C:/Users/rikke007/OneDrive - WageningenUR/Documents/WLR/KB_BEC3/Data/Weather_data/Leeuwarden_weather_station_hourly_data_20210101_20220924.txt",
sep = ",", header = TRUE)
colnames(wdata)
# THI = (1.8 × Tmean + 32) − [(0.55 − 0.0055 × RHmean) × (1.8 × Tmean − 26)],
# where Tmean is the daily average T (°C) and RHmean is daily mean RH (%).
# calculate mean temp and mean RH per day
meanT <- aggregate((wdata[, 8] * 0.1), list(wdata$YYYYMMDD), FUN = mean)
colnames(meanT) <- c("date", "meanT")
meanRH <- aggregate(wdata[, 18], list(wdata$YYYYMMDD), FUN = mean)
colnames(meanRH) <- c("date", "meanRH")
# create meanTHI
meanTHI <- merge(meanT, meanRH, by = "date")
meanTHI$meanTHI <- (1.8 * meanTHI$meanT + 32) - (0.55 - 0.0055 * meanTHI$meanRH) * (1.8 * meanTHI$meanT - 26)
write.csv(meanTHI, "C:/Users/rikke007/OneDrive - WageningenUR/Documents/WLR/KB_BEC3/Data/Weather_data/daily_THI_01012021_24092022.txt",
row.names = FALSE)
######################################################
# ANALYSIS FOR BEC-3 IN 2022
# Author: Roxann Rikkers, WLR-ABG
# Date: 2022
# Descriptive statistics of % area usage and transformed weather data
# datafile from script Make_table_all_results.R
# distance moved and percentage active per day not used in downstream analysis,
# not corrected for missing data/gaps in the data.
######################################################
rm(list = ls())
setwd("C:/Users/rikke007/OneDrive - WageningenUR/Documents/WLR/KB_BEC3/")
#### Packages #############################################################################
library(dplyr)
library(reshape2)
library(ggplot2)
### Visualize data #######################################################
fdata <- read.csv("Data/Final_data/area_activity_weather_data_BEC3_2022-09-27.txt")
#######################################################
#### Plot area usage against mean THI #################
### remove cows with > 50% missing records, this removes 5100 - 4606 = 494 rows
datawide <- fdata[fdata$Missrec < 43200, ]
# get needed datavariables
datawide <- datawide[, c(1:2, 5, 9, 10, 11, 12, 25)]
# round THI's to integers (no decimal)
datawide$meanTHI <- round(datawide$meanTHI)
# aggregate results over rounded THIs per day per area
avg <- aggregate(datawide[, 3:8], list(datawide$meanTHI), FUN = mean)
stds <- aggregate(datawide[, 3:8], list(datawide$meanTHI), FUN = sd)
# datalong of averages
datalong <- melt(avg[1:6], id.vars = c("Group.1"))
# datalong of STDs
datalong2 <- melt(stds[1:6], id.vars = c("Group.1"))
## merge 2 datalong files
datalongfinal <- cbind(datalong, datalong2[, 3])
colnames(datalongfinal) <- c("meanTHI", "variable", "means", "stds")
# use colorblind palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
plot <- ggplot(data = datalongfinal, aes(x = meanTHI,
y = means,
group = variable,
colour = variable)) +
geom_line(size = 1) +
ylab("% of records per day") +
xlab("meanTHI") +
theme(axis.text.x = element_text(angle = 90)) +
scale_colour_manual(name = "Area", labels = c("Concentrate Feed", "Drinking Through", "Feeding rack", "No specific area",
"Cubicles"),
values = cbPalette)
plot
# save plot
ggsave(plot = plot, filename = "Results/Area_usage/Area_usage_meanTHI_nobars.png",
width = 8, height = 4)
plot <- ggplot(data = datalongfinal, aes(x = meanTHI,
y = means,
group = variable,
colour = variable)) +
geom_errorbar(aes(ymin = means - stds, ymax = means + stds), width = 0.5) +
geom_line(size = 1) +
ylab("% of records per day") +
xlab("meanTHI") +
theme(axis.text.x = element_text(angle = 90)) +
scale_colour_manual(name = "Area", labels = c("Concentrate Feed", "Drinking Through", "Feeding rack", "No specific area",
"Cubicles"),
values = cbPalette)
plot
# save plot
ggsave(plot = plot, filename = "Results/Area_usage/Area_usage_meanTHI_withbars.png",
width = 8, height = 4)
#############################################################################
#### Plot area usage per date (mean THI included as second y axis) ##########
### remove cows with > 50% missing records, this removes 5100 - 4606 = 494 rows
datawide <- fdata[fdata$Missrec < 43200, ]
# get needed datavariables
datawide <- datawide[, c(1:2, 5, 9, 10, 11, 12, 25)]
# round THI's to integers (no decimal)
datawide$meanTHI <- round(datawide$meanTHI)
# aggregate results over rounded THIs per day per area
avg <- aggregate(datawide[, 3:8], list(datawide$date), FUN = mean)
stds <- aggregate(datawide[, 3:8], list(datawide$date), FUN = sd)
# datalong of averages
datalong <- melt(avg[1:7], id.vars = c("Group.1"))
# datalong of STDs
datalong2 <- melt(stds[1:7], id.vars = c("Group.1"))
## merge 2 datalong files
datalongfinal <- cbind(datalong, datalong2[, 3])
colnames(datalongfinal) <- c("date", "variable", "means", "stds")
datalongfinal$date <- as.Date(x = as.character(datalongfinal$date), "%Y%m%d")
### plot per date
plot <- ggplot(data = datalongfinal, aes(x = date,
y = means,
group = variable,
colour = variable)) +
geom_line(size = 1) +
xlab("") +
theme(axis.text.x = element_text(angle = 90)) +
scale_colour_manual(name = "", labels = c("Concentrate Feed", "Drinking Through", "Feeding rack", "No specific area",
"Cubicles", "Mean THI"),
values = cbPalette) +
scale_y_continuous("% of records per day", sec.axis = sec_axis(~ . * 1, name = "mean THI")
)
plot
# save plot
ggsave(plot = plot, filename = "Results/Area_usage/Area_usage_DATES_nobars.png",
width = 16, height = 4)
plot <- ggplot(data = datalongfinal, aes(x = date,
y = means,
group = variable,
colour = variable)) +
geom_errorbar(aes(ymin = means - stds, ymax = means + stds), width = 0.5) +
geom_line(size = 1) +
xlab("") +
theme(axis.text.x = element_text(angle = 90)) +
scale_colour_manual(name = "", labels = c("Concentrate Feed", "Drinking Through", "Feeding rack", "No specific area",
"Cubicles", "Mean THI"),
values = cbPalette) +
scale_y_continuous("% of records per day", sec.axis = sec_axis(~ . * 1, name = "mean THI")
)
plot
# save plot
ggsave(plot = plot, filename = "Results/Area_usage/Area_usage_DATES_withbars.png",
width = 16, height = 4)
###############################################
#### plot individual cows #####################
dir.create(path = "Results/Area_usage/Individual_cows")
for(cow in unique(fdata$cowid)) { # cow = 443
print(cow)
tmpcow <- fdata[fdata$cowid == cow, ]
tmpcow$percmissrec <- (tmpcow$Missrec / 86399) * 100
datawide <- tmpcow[, c(1:2, 5, 9, 10, 11, 12, 26)]
datalong <- melt(datawide, id.vars = c("date", "cowid"))
datalong$date <- as.Date(x = as.character(datalong$date), "%Y%m%d")
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
plot <- ggplot(data = datalong, aes(x = date,
y = value,
group = variable,
colour = variable)) +
geom_line(size = 1) +
ylab("% of records per day") +
xlab("") +
theme(axis.text.x = element_text(angle = 90)) +
scale_colour_manual(name = "Area", labels = c("Concentrate Feed", "Drinking through", "Feeding rack",
"No specific area", "Cubicles", "% missing records"),
values = cbPalette) +
ylab("% records per day") +
ggtitle(paste0("Area usage cow ", cow))
plot
# save plot
ggsave(plot = plot, filename = paste0("Results/Area_usage/Individual_cows/Area_usage_missrec_cow", cow, ".png"),
width = 8, height = 4)
tmpcow <- fdata[fdata$cowid == cow, ]
datawide <- tmpcow[, c(1:2, 5, 9, 10, 11, 12, 25)]
datalong <- melt(datawide, id.vars = c("date", "cowid"))
datalong$date <- as.Date(x = as.character(datalong$date), "%Y%m%d")
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
plot <- ggplot(data = datalong, aes(x = date,
y = value,
group = variable,
colour = variable)) +
geom_line(size = 1) +
ylab("% of records per day") +
xlab("") +
theme(axis.text.x = element_text(angle = 90)) +
scale_colour_manual(name = "Area", labels = c("Concentrate Feed", "Drinking through", "Feeding rack",
"No specific area", "Cubicles", "meanTHI"),
values = cbPalette) +
ylab("% records per day") +
ggtitle(paste0("Area usage cow ", cow))
plot
# save plot
ggsave(plot = plot, filename = paste0("Results/Area_usage/Individual_cows/Area_usage_meanTHI_cow", cow, ".png"),
width = 8, height = 4)
}
####################################
#### Plot activity data ############
rm(list = ls())
fdata <- read.csv("Data/Final_data/area_activity_weather_data_BEC3_2022-09-27.txt")
colnames(fdata)
######################################################
#### Plot percentage active againt meanTHI ###########
# get needed datavariables
datawide <- fdata[, c(1:2, 17, 21, 22, 25)]
# may 30th is missing, so for now remove this data
#datawide <- datawide[!datawide$date == "20220530",]
# round THI's to integers (no decimal)
datawide$meanTHI <- round(datawide$meanTHI)
# aggregate resilts over rounded THIs per day per area
avg <- aggregate(datawide[, 3:5], list(datawide$meanTHI), FUN = mean)
stds <- aggregate(datawide[, 3:5], list(datawide$meanTHI), FUN = sd)
# datalong of averages
datalong <- melt(avg[1:2], id.vars = c("Group.1"))
# datalong of STDs
datalong2 <- melt(stds[1:2], id.vars = c("Group.1"))
## merge 2 datalong files
datalongfinal <- cbind(datalong, datalong2[, 3])
colnames(datalongfinal) <- c("meanTHI", "variable", "means", "stds")
# use colorblind palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
plot <- ggplot(data = datalongfinal, aes(x = meanTHI,
y = means,
group = variable,
colour = variable)) +
geom_line(size = 1) +
ylab("% active") +
xlab("meanTHI") +
theme(axis.text.x = element_text(angle = 90))
plot
# save plot
ggsave(plot = plot, filename = "Results/Activity/PercActive_meanTHI_nobars.png",
width = 8, height = 4)
plot <- ggplot(data = datalongfinal, aes(x = meanTHI,
y = means,
group = variable,
colour = variable)) +
geom_errorbar(aes(ymin = means - stds, ymax = means + stds), width = 0.5) +
geom_line(size = 1) +
ylab("% of records per day") +
xlab("meanTHI") +
theme(axis.text.x = element_text(angle = 90))
plot
# save plot
ggsave(plot = plot, filename = "Results/Activity/PercActive_meanTHI_withbars.png",
width = 8, height = 4)
######################################################
#### Plot distance moved against meanTHI #############
# datalong of averages
datalong <- melt(avg[c(1,4)], id.vars = c("Group.1"))
# datalong of STDs
datalong2 <- melt(stds[c(1,4)], id.vars = c("Group.1"))
## merge 2 datalong files
datalongfinal <- cbind(datalong, datalong2[, 3])
colnames(datalongfinal) <- c("meanTHI", "variable", "means", "stds")
# use colorblind palette
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
plot <- ggplot(data = datalongfinal, aes(x = meanTHI,
y = means,
group = variable,
colour = variable)) +
geom_line(size = 1) +
ylab("Distance Moved") +
xlab("meanTHI") +
theme(axis.text.x = element_text(angle = 90))
plot
# save plot
ggsave(plot = plot, filename = "Results/Activity/DistanceMoved_meanTHI_nobars.png",
width = 8, height = 4)
plot <- ggplot(data = datalongfinal, aes(x = meanTHI,
y = means,
group = variable,
colour = variable)) +
geom_errorbar(aes(ymin = means - stds, ymax = means + stds), width = 0.5) +
geom_line(size = 1) +
ylab("Distance Moved") +
xlab("meanTHI") +
theme(axis.text.x = element_text(angle = 90))
plot
# save plot
ggsave(plot = plot, filename = "Results/Activity/DistanceMoved_meanTHI_withbars.png",
width = 8, height = 4)
####################################################################################
#### Plot percentage active per date (mean THI included as second y axis) ##########
rm(list = ls())
fdata <- read.csv("Data/Final_data/area_activity_weather_data_BEC3_2022-09-27.txt")
colnames(fdata)
# get needed datavariables
datawide <- fdata[, c(1:2, 17, 21, 22, 25)]
# may 30th is missing, so for now remove this data
#datawide <- datawide[!datawide$date == "20220530",]
# round THI's to integers (no decimal)
datawide$meanTHI <- round(datawide$meanTHI)
# aggregate resilts over rounded THIs per day per area
avg <- aggregate(datawide[, 3:5], list(datawide$date), FUN = mean)
stds <- aggregate(datawide[, 3:5], list(datawide$date), FUN = sd)
# some meanTHI only 1 record, so no SD and gives NA values...change NA values to 0
#stds[is.na(stds)] <- 0
# datalong of averages
datalong <- melt(avg[1:2], id.vars = c("Group.1"))
# datalong of STDs
datalong2 <- melt(stds[1:2], id.vars = c("Group.1"))
## merge 2 datalong files
datalongfinal <- cbind(datalong, datalong2[, 3])
colnames(datalongfinal) <- c("date", "variable", "means", "stds")
datalongfinal$date <- as.Date(x = as.character(datalongfinal$date), "%Y%m%d")
### plot per date
plot <- ggplot(data = datalongfinal, aes(x = date,
y = means,
group = variable,
colour = variable)) +
geom_line(size = 1) +
xlab("") + ylab("% active") +
theme(axis.text.x = element_text(angle = 90))
plot
# save plot
ggsave(plot = plot, filename = "Results/Activity/PercActive_DATES_nobars.png",
width = 16, height = 4)
plot <- ggplot(data = datalongfinal, aes(x = date,
y = means,
group = variable,
colour = variable)) +
geom_errorbar(aes(ymin = means - stds, ymax = means + stds), width = 0.5) +
geom_line(size = 1) +
xlab("") + ylab("% active") +
theme(axis.text.x = element_text(angle = 90))
plot
# save plot
ggsave(plot = plot, filename = "Results/Activity/PercActive_DATES_withbars.png",
width = 16, height = 4)
#################################################################################
#### Plot distance moved per date (mean THI included as second y axis) ##########
rm(list = ls())
fdata <- read.csv("Data/Final_data/area_activity_weather_data_BEC3_2022-09-27.txt")
colnames(fdata)
# get needed datavariables
datawide <- fdata[, c(1:2, 17, 21, 22, 25)]
# may 30th is missing, so for now remove this data
#datawide <- datawide[!datawide$date == "20220530",]
# round THI's to integers (no decimal)
datawide$meanTHI <- round(datawide$meanTHI)
# aggregate resilts over rounded THIs per day per area
avg <- aggregate(datawide[, 3:5], list(datawide$date), FUN = mean)
stds <- aggregate(datawide[, 3:5], list(datawide$date), FUN = sd)
# some meanTHI only 1 record, so no SD and gives NA values...change NA values to 0
#stds[is.na(stds)] <- 0
# datalong of averages
datalong <- melt(avg[c(1,4)], id.vars = c("Group.1"))
# datalong of STDs
datalong2 <- melt(stds[c(1,4)], id.vars = c("Group.1"))
## merge 2 datalong files
datalongfinal <- cbind(datalong, datalong2[, 3])
colnames(datalongfinal) <- c("date", "variable", "means", "stds")
datalongfinal$date <- as.Date(x = as.character(datalongfinal$date), "%Y%m%d")
### plot per date
plot <- ggplot(data = datalongfinal, aes(x = date,
y = means,
group = variable,
colour = variable)) +
geom_line(size = 1) +
xlab("") + ylab("Distance Moved") +
theme(axis.text.x = element_text(angle = 90))
plot
# save plot
ggsave(plot = plot, filename = "Results/Activity/DistanceMoved_DATES_nobars.png",
width = 16, height = 4)
plot <- ggplot(data = datalongfinal, aes(x = date,
y = means,
group = variable,
colour = variable)) +
geom_errorbar(aes(ymin = means - stds, ymax = means + stds), width = 0.5) +
geom_line(size = 1) +
xlab("") + ylab("Distance Moved") +
theme(axis.text.x = element_text(angle = 90))
plot
# save plot
ggsave(plot = plot, filename = "Results/Activity/DistanceMoved_DATES_withbars.png",
width = 16, height = 4)
###############################################
#### plot individual cows #####################
rm(list = ls())
dir.create("Results/Activity/Individual_cows")
inddata <- read.csv("Data/Final_data/area_activity_weather_data_BEC3_2022-09-27.txt")
inddata$meanTHI <- round(inddata$meanTHI)
for(cow in unique(inddata$cowid)) { # cow = 280
print(cow)
tmpcow <- inddata[inddata$cowid == cow, ]
tmpcow$date <- as.Date(x = as.character(tmpcow$date), "%Y%m%d")
plot <- ggplot(data = tmpcow, aes(x = date,
y = perc_active_1min),
group = 1) +
geom_line(size = 1) +
ylab("% active") +
xlab("") +
theme(axis.text.x = element_text(angle = 90)) +
ggtitle(paste0("% active cow ", cow))
plot
# save plot
ggsave(plot = plot, filename = paste0("Results/Activity/Individual_cows/PercActive_cow", cow, ".png"),
width = 8, height = 4)
plot <- ggplot(data = tmpcow, aes(x = date,
y = dist_moved)) +
geom_line(size = 1) +
ylab("Distance Moved") +
xlab("") +
ggtitle(paste0("Distance Moved cow ", cow))
plot
# save plot
ggsave(plot = plot, filename = paste0("Results/Activity/Individual_cows/DistanceMoved_cow", cow, ".png"),
width = 8, height = 4)
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment