library(MASS) library(ggplot2) require(MESS) require(Bolstad2) library(segmented) library(grDevices) #0 data loading #---- # @Yangkou 2016 wader_YK <- read.csv("/Users/tong/Desktop/Chapter 1 local distribution/Yangkou_2016F.csv"); plottitle <- "Yangkou"; wader_YK[is.na(wader_YK)] <- 0 # @Nanpu 2017 wader_NP <- read.csv("/Users/tong/Desktop/Chapter 1 local distribution/Nanpu_2017S.csv"); plottitle <- "Nanpu"; wader_NP[is.na(wader_NP)] <- 0 #matrix converting plot numbers to median distances for each transect plot_distance <- read.csv("/Users/tong/Desktop/Chapter 1 local distribution/plot_to_distance.csv") #1 get the proportion of exposed mudflat estimated from observation time#---- #1.1 Get the moment of exposure and submersion exposure_momentYK <- data.frame(XF = c(1:9), XF_max = c(1:9), SP = c(1:9), SP_max = c(1:9)) exposure_momentNP <- data.frame(NPT1 = c(1:9), NPT1_max = c(1:9), NPT2 = c(1:9), NPT2_max = c(1:9), NPT3 = c(1:9), NPT3_max = c(1:9), NPT4 = c(1:9), NPT4_max = c(1:9)) wader <- wader_YK for (n in 1:(ncol(exposure_momentYK)/2)) { #exposure time approximated by observation time for (i in 1:9) { exposure_momentYK[i, 2*n-1] <- min(wader$TimeSinceHighTide[wader["plot"]==i&wader["transect"]==colnames(exposure_momentYK)[2*n-1]]) exposure_momentYK[i, 2*n] <- max(wader$TimeSinceHighTide[wader["plot"]==i&wader["transect"]==colnames(exposure_momentYK)[2*n-1]]) } } exposure_momentYK$XF[4] <- NA; exposure_momentYK$XF[5] <- NA # obvious outliers excluded wader <- wader_NP for (n in 1:(ncol(exposure_momentNP)/2)) { #exposure time approximated by observation time for (i in 1:9) { exposure_momentNP[i, 2*n-1] <- min(wader$TimeSinceHighTide[wader["plot"]==i&wader["transect"]==colnames(exposure_momentNP)[2*n-1]]) exposure_momentNP[i, 2*n] <- max(wader$TimeSinceHighTide[wader["plot"]==i&wader["transect"]==colnames(exposure_momentNP)[2*n-1]]) } } #1.2 linear regression to calculate the relationship between time and tide location (distance) #note the predicting factor here is time, and the response is distance (not the actual causal relationship) expo_coefYK <- data.frame(XF = c(1:6), SP = c(1:6), row.names = c("slope", "intercept", "p", "max_slope", "max_intercept", "max_p")) expo_coefNP <- data.frame(NPT1 = c(1:6), NPT2 = c(1:6), NPT3 = c(1:6), NPT4 = c(1:6), row.names = c("slope", "intercept", "p", "max_slope", "max_intercept", "max_p")) for (n in 1:ncol(expo_coefYK)) { distance <- (plot_distance[, colnames(expo_coefYK)[n]] -250) # max and min are both the upper edge of the plot # pay attention to the x and y values. Using time to predict distances lin.mod <- lm (distance ~ exposure_momentYK[, (2*n-1)]) expo_coefYK[1, n] <- coef(lin.mod)[2]; expo_coefYK[2, n] <- coef(lin.mod)[1] expo_coefYK[3, n] <- paste0("p=", round(summary(lin.mod)$coefficient[2, 4], digits = 4), ", rsqr=", round(summary(lin.mod)$r.squared, digits = 4)) lin.mod <- lm (distance ~ exposure_momentYK[, (2*n)]) expo_coefYK[4, n] <- coef(lin.mod)[2]; expo_coefYK[5, n] <- coef(lin.mod)[1] expo_coefYK[6, n] <- paste0("p=", round(summary(lin.mod)$coefficient[2, 4], digits = 4), ", rsqr=", round(summary(lin.mod)$r.squared, digits = 4)) } for (n in 1:ncol(expo_coefNP)) { distance <- (plot_distance[, colnames(expo_coefNP)[n]] -125) # max and min are both the upper edge of the plot # pay attention to the x and y values. Using time to predict distances lin.mod <- lm (distance ~ exposure_momentNP[, (2*n-1)]) expo_coefNP[1, n] <- coef(lin.mod)[2]; expo_coefNP[2, n] <- coef(lin.mod)[1] expo_coefNP[3, n] <- paste0("p=", round(summary(lin.mod)$coefficient[2, 4], digits = 4), ", rsqr=", round(summary(lin.mod)$r.squared, digits = 4)) lin.mod <- lm (distance ~ exposure_momentNP[, (2*n)]) expo_coefNP[4, n] <- coef(lin.mod)[2]; expo_coefNP[5, n] <- coef(lin.mod)[1] expo_coefNP[6, n] <- paste0("p=", round(summary(lin.mod)$coefficient[2, 4], digits = 4), ", rsqr=", round(summary(lin.mod)$r.squared, digits = 4)) } #1.3 calculate the proportion of mudflat expose using the time of observation #and linear movement of tide (from the relationship calculated above) for (n in 1:nrow(wader_YK)) { trans_name <- as.character(wader_YK$transect[n]) if (wader_YK$TimeSinceHighTide[n] < (exposure_momentYK[9, trans_name] + exposure_momentYK[9, paste0(trans_name, "_max")])/2) { wader_YK$prop[n] <- ((wader_YK$TimeSinceHighTide[n] * as.numeric(expo_coefYK[1, trans_name]) + as.numeric(expo_coefYK[2, trans_name])) / (plot_distance[9, trans_name] + 250)) } else { wader_YK$prop[n] <- ((wader_YK$TimeSinceHighTide[n] * as.numeric(expo_coefYK[4, trans_name]) + as.numeric(expo_coefYK[5, trans_name])) / (plot_distance[9, trans_name] + 250)) } if (wader_YK$prop[n] > 1) wader_YK$prop[n] <- 1 if (wader_YK$prop[n] <= 0) wader_YK$prop[n] <- 0.0001 } for (n in 1:nrow(wader_NP)) { trans_name <- as.character(wader_NP$transect[n]) if (wader_NP$TimeSinceHighTide[n] < (exposure_momentNP[9, trans_name] + exposure_momentNP[9, paste0(trans_name, "_max")])/2) { wader_NP$prop[n] <- ((wader_NP$TimeSinceHighTide[n] * as.numeric(expo_coefNP[1, trans_name]) + as.numeric(expo_coefNP[2, trans_name])) / (plot_distance[9, trans_name] + 125)) } else { wader_NP$prop[n] <- ((wader_NP$TimeSinceHighTide[n] * as.numeric(expo_coefNP[4, trans_name]) + as.numeric(expo_coefNP[5, trans_name])) / (plot_distance[9, trans_name] + 125)) } if (wader_NP$prop[n] > 1) wader_NP$prop[n] <- 1 if (wader_NP$prop[n] <= 0) wader_NP$prop[n] <- 0.0001 } #2 calculate the centroids for each time/exposure period (total 8/.125% each) for each transect #---- centroid_NP <- list("GreyP", "KentishP", "Dunlin", "Stint", "Sanderling", "GreatK", "CurlewS", "RedK", "Curlew", "Godwit", "Greenshanks", "Terek", "Turnstone") names(centroid_NP) <- c("GreyP", "KentishP", "Dunlin", "Stint", "Sanderling", "GreatK", "CurlewS", "RedK", "Curlew", "Godwit", "Greenshanks", "Terek", "Turnstone") centroid_YK <- list("GreyP", "KentishP", "SandP", "Dunlin", "Stint", "Knots", "Curlew", "Godwit", "Greenshanks", "Marsh", "Terek", "Turnstone", "Tattler") names(centroid_YK) <- c("GreyP", "KentishP", "SandP", "Dunlin", "Stint", "Knots", "Curlew", "Godwit", "Greenshanks", "Marsh", "Terek", "Turnstone", "Tattler") for (species in names(centroid_YK)) { time_segments <- 8 # the number of bins for the proportions to group to temp_centroid <- data.frame(YK = c(1:time_segments)) centroid_YK[[species]] <- list(centroid = temp_centroid) for (n in colnames(temp_centroid)) { plot_counts <- data.frame(per1 = c(1:9)) plot_counts.plot <- data.frame(per1 = c(1:9)) for (i in 1:time_segments) { for (p_num in 1:nrow(plot_counts)) { plot_counts[p_num, paste0("per", i)] <- mean(wader_YK[wader_YK["plot"]==p_num & wader_YK["prop"]<=(i/time_segments) & wader_YK["prop"]>((i-1)/time_segments), species]) plot_counts.plot[p_num, paste0("per", i)] <- plot_counts[p_num, paste0("per", i)] * plot_distance$Yangkou[p_num] } temp_centroid[i, n] <- sum(plot_counts.plot[ , paste0("per", i)], na.rm = TRUE) / sum(plot_counts[ , paste0("per", i)], na.rm = TRUE) } centroid_YK[[species]][n] <- list(plot_counts) } centroid_YK[[species]]$centroid <- temp_centroid } for (species in names(centroid_NP)) { time_segments <- 8 # the number of bins for the proportions to group to temp_centroid <- data.frame(NP = c(1:time_segments)) centroid_NP[[species]] <- list(centroid = temp_centroid) for (n in colnames(temp_centroid)) { plot_counts <- data.frame(per1 = c(1:9)) plot_counts.plot <- data.frame(per1 = c(1:9)) for (i in 1:time_segments) { for (p_num in 1:nrow(plot_counts)) { plot_counts[p_num, paste0("per", i)] <- mean(wader_NP[wader_NP["plot"]==p_num & wader_NP["prop"]<=(i/time_segments) & wader_NP["prop"]>((i-1)/time_segments), species]) plot_counts.plot[p_num, paste0("per", i)] <- plot_counts[p_num, paste0("per", i)] * plot_distance$Nanpu[p_num] } temp_centroid[i, n] <- sum(plot_counts.plot[ , paste0("per", i)], na.rm = TRUE) / sum(plot_counts[ , paste0("per", i)], na.rm = TRUE) } centroid_NP[[species]][n] <- list(plot_counts) } centroid_NP[[species]]$centroid <- temp_centroid } #3 calculate and test the slope between proportion of mudflat exposed (x) and centroid of birds (y) #---- centroid_test_YK <- data.frame() ## "slope", "slope_se", "r2", "AIC", "spec_slopes", "breakpoint" "spec_AIC" datasheet <- centroid_YK spp <- c("GreyP", "KentishP", "SandP", "Dunlin", "Stint", "Knots", "Curlew", "Godwit", "Greenshanks", "Marsh", "Terek", "Turnstone", "Tattler") for (species in spp) { # x <- rep(c(1:time_segments)/time_segments-1/time_segments/2, 2) #median time when the centroids were calculated # y <- c(datasheet[[species]]$centroid$XF, datasheet[[species]]$centroid$SP) / 9 # normalize to 0-1 x <- rep(c(1:time_segments)/time_segments - 1/time_segments/2) #median time when the centroids were calculated y <- c(datasheet[[species]]$centroid$YK) / (plot_distance$Yangkou[9] + 250) # normalize to 0-1 lin.mod <- lm(y ~ 0+x) temp <- summary(lin.mod) centroid_test_YK["slope", species] <- temp$coefficients[1] centroid_test_YK["slope_se", species] <- temp$coefficients[2] centroid_test_YK["r2, p", species] <- paste0("r2=", round(temp$r.squared, digits = 4), "; p=", round(temp$coefficients[4], digits = 4)) centroid_test_YK["AIC", species] <- AIC(lin.mod) centroid_test_YK["BIC", species] <- BIC(lin.mod) spec.mod <- 1 for (i in 1:20) try(spec.mod <- segmented(lin.mod, seg.Z = ~x), silent = TRUE) try(centroid_test_YK["spec_slopes", species] <- paste0("k1=", round(slope(spec.mod)$x[1, 1], digits = 4), "; se1=", round(slope(spec.mod)$x[1, 2], digits = 4), "; k2=", round(slope(spec.mod)$x[2, 1], digits = 4), "; se2=", round(slope(spec.mod)$x[2, 2], digits = 4))) try(centroid_test_YK["breakpoint", species] <- paste0("point=", spec.mod$psi[2], ", se=", spec.mod$psi[3])) try(centroid_test_YK["r2", species] <- round(summary(spec.mod)$r.squared, digits = 4)) try(centroid_test_YK["spec_AIC", species] <- AIC(spec.mod)) try(centroid_test_YK["spec_BIC", species] <- BIC(spec.mod)) pdf(paste0(species, "_YK_centroid.pdf")) plot(x, y, main = paste(species, "YK"), xlab = "proportion of mudflat exposed", xlim = c(0, 1), ylab = "normalized centroid position", ylim = c(0, 1)) abline(0, 1, col = "blue"); abline(0, 0.5, col = "yellow") try(segments(0, 0, spec.mod$psi[2], spec.mod$psi[2]*slope(spec.mod)$x[1, 1], col = "red")) try(segments(spec.mod$psi[2], spec.mod$psi[2]*slope(spec.mod)$x[1, 1], 1.2, spec.mod$psi[2]*slope(spec.mod)$x[1, 1] + (1.2-spec.mod$psi[2])* slope(spec.mod)$x[2, 1], col = "red")) try(segments(spec.mod$psi[2]-1.96*spec.mod$psi[3], spec.mod$psi[2]*slope(spec.mod)$x[1, 1], spec.mod$psi[2]+1.96*spec.mod$psi[3], spec.mod$psi[2]*slope(spec.mod)$x[1, 1], col = "grey")) abline(0, temp$coefficients[1], col = "black") abline(0, temp$coefficients[1]+1.96*temp$coefficients[2], col = "gray") abline(0, temp$coefficients[1]-1.96*temp$coefficients[2], col = "gray") dev.off() } write.csv(centroid_test_YK, file = "test_resultYK.csv") centroid_test_NP <- data.frame() ## "slope", "slope_se", "r2", "AIC", "spec_slopes", "breakpoint" "spec_AIC" datasheet <- centroid_NP spp <- c("GreyP", "KentishP", "Dunlin", "Stint", "Sanderling", "GreatK", "CurlewS", "RedK", "Curlew", "Godwit", "Greenshanks", "Terek", "Turnstone") for (species in spp) { # x <- rep(c(1:time_segments)/time_segments-1/time_segments/2, 4) #median time when the centroids were calculated # y <- c(datasheet[[species]]$centroid$NPT1, datasheet[[species]]$centroid$NPT2, # datasheet[[species]]$centroid$NPT3, datasheet[[species]]$centroid$NPT4) / 9 # normalize to 0-1 x <- rep(c(1:time_segments)/time_segments-1/time_segments/2) #median time when the centroids were calculated y <- c(datasheet[[species]]$centroid$NP) / (plot_distance$Nanpu[9] + 125) # normalize to 0-1 lin.mod <- lm(y ~ 0+x) temp <- summary(lin.mod) centroid_test_NP["slope", species] <- temp$coefficients[1] centroid_test_NP["slope_se", species] <- temp$coefficients[2] centroid_test_NP["r2, p", species] <- paste0("r2=", round(temp$r.squared, digits = 4), "; p=", round(temp$coefficients[4], digits = 4)) centroid_test_NP["AIC", species] <- AIC(lin.mod) centroid_test_NP["BIC", species] <- BIC(lin.mod) spec.mod <- 1 for (i in 1:20) try(spec.mod <- segmented(lin.mod, seg.Z = ~x), silent = TRUE) try(centroid_test_NP["spec_slopes", species] <- paste0("k1=", round(slope(spec.mod)$x[1, 1], digits = 4), "; se1=", round(slope(spec.mod)$x[1, 2], digits = 4), "; k2=", round(slope(spec.mod)$x[2, 1], digits = 4), "; se2=", round(slope(spec.mod)$x[2, 2], digits = 4))) try(centroid_test_NP["breakpoint", species] <- paste0("point=", spec.mod$psi[2], ", se=", spec.mod$psi[3])) try(centroid_test_NP["r2", species] <- round(summary(spec.mod)$r.squared, digits = 4)) try(centroid_test_NP["spec_AIC", species] <- AIC(spec.mod)) try(centroid_test_NP["spec_BIC", species] <- BIC(spec.mod)) pdf(paste0(species, "_NP_centroid.pdf")) plot(x, y, main = paste(species, "NP"), xlab = "proportion of mudflat exposed", xlim = c(0, 1), ylab = "normalized centroid position", ylim = c(0, 1)) abline(0, 1, col = "blue"); abline(0, 0.5, col = "yellow") try(segments(0, 0, spec.mod$psi[2], spec.mod$psi[2]*slope(spec.mod)$x[1, 1], col = "red")) try(segments(spec.mod$psi[2], spec.mod$psi[2]*slope(spec.mod)$x[1, 1], 1.2, spec.mod$psi[2]*slope(spec.mod)$x[1, 1] + (1.2-spec.mod$psi[2])* slope(spec.mod)$x[2, 1], col = "red")) try(segments(spec.mod$psi[2]-1.96*spec.mod$psi[3], spec.mod$psi[2]*slope(spec.mod)$x[1, 1], spec.mod$psi[2]+1.96*spec.mod$psi[3], spec.mod$psi[2]*slope(spec.mod)$x[1, 1], col = "grey")) abline(0, temp$coefficients[1], col = "black") abline(0, temp$coefficients[1]+1.96*temp$coefficients[2], col = "gray") abline(0, temp$coefficients[1]-1.96*temp$coefficients[2], col = "gray") dev.off() } write.csv(centroid_test_NP, file = "test_resultNP.csv") #4 make the distribution plot for each species #---- wader <- wader_NP; plottitle <- "NP" sppNP <- c("GreyP", "KentishP", "Dunlin", "Stint", "Sanderling", "GreatK", "CurlewS", "RedK", "Curlew", "Godwit", "Greenshanks", "Terek", "Turnstone") spp_nameNP <- data.frame(GreyP = "Grey Plover", KentishP = "Kentish Plover", Dunlin = "Dunlin", Stint = "Red-necked Stint", Sanderling = "Sanderling", GreatK = "Great Knot", CurlewS = "Curlew Sandpiper", RedK = "Red Knot", Curlew = "Curlews", Godwit = "Bar-tailed Godwit", Greenshanks = "Common Greenshank", Terek = "Terek Sandpiper", Turnstone = "Ruddy Turnstone") spp_nameNP <- rbind(spp_nameNP, setNames(data.frame("S", "S", "T", "N", "T", "T", "S", "G", "N", "S", "S", "S", "G"), sppNP)) for (species in sppNP) { legend <- factor(as.integer(wader[, species] > 0)) if (spp_nameNP[2, species] == "G") {forage_type <- "generalist"} if (spp_nameNP[2, species] == "S") {forage_type <- "zone specialist"} if (spp_nameNP[2, species] == "T") {forage_type <- "tide follower"} if (spp_nameNP[2, species] == "N") {forage_type <- "not assigned"} p1 <- ggplot(wader, aes(Distance, prop)) + labs(title = paste0(spp_nameNP[1, species], " at Nanpu: ", forage_type), x = "distance from sea wall", y = "proportion of tidal flat exposed") + # scale_color_grey(start = 0.75, end = 0, name = "") scale_color_discrete(name = "") p2 <- p1 + geom_point(aes(size = wader[, species], color = legend, shape = legend)) + scale_shape_manual(values = c("circle open", "circle")) pdf(file = paste0(species, "_", plottitle, "_distribution.pdf"), width = 6, height = 5) plot(p2) dev.off() } wader <- wader_YK; plottitle <- "YK" sppYK <- c("GreyP", "KentishP", "SandP", "Dunlin", "Stint", "Knots", "Curlew", "Godwit", "Greenshanks", "Marsh", "Terek", "Turnstone", "Tattler") spp_nameYK <- data.frame(GreyP = "Grey Plover", KentishP = "Kentish Plover", SandP = "Sand Plovers", Dunlin = "Dunlin", Stint = "Red-necked Stint", Knots = "Knots", Curlew = "Curlews", Godwit = "Bar-tailed Godwit", Greenshanks = "Common Greenshank", Marsh = "Marsh Sandpiper", Terek = "Terek Sandpiper", Turnstone = "Ruddy Turnstone", Tattler = "Grey-tailed Tattler") spp_nameYK <- rbind(spp_nameYK, setNames(data.frame("G", "G", "G", "S", "G", "N", "G", "S", "G", "S", "S", "S", "N"), sppYK)) for (species in sppYK) { legend <- factor(as.integer(wader[, species] > 0)) if (spp_nameYK[2, species] == "G") {forage_type <- "generalist"} if (spp_nameYK[2, species] == "S") {forage_type <- "zone specialist"} if (spp_nameYK[2, species] == "T") {forage_type <- "tide follower"} if (spp_nameYK[2, species] == "N") {forage_type <- "not assigned"} p1 <- ggplot(wader, aes(Distance, prop)) + labs(title = paste0(spp_nameYK[1, species], " at Rudong: ", forage_type), x = "distance from sea wall", y = "proportion of tidal flat exposed") + # scale_color_grey(start = 0.75, end = 0, name = "") scale_color_discrete(name = "") p2 <- p1 + geom_point(aes(size = wader[, species], color = legend, shape = legend)) + scale_shape_manual(values = c("circle open", "circle")) pdf(file = paste0(species, "_", plottitle, "_distribution.pdf"), width = 6, height = 5) plot(p2) dev.off() } #5 calculate the cumulative feeding time for each species in each transect #---- wader <- wader_YK; plottitle <- "YK" feeding_time <- data.frame(XF = c(1:9), SP = c(1:9)) spp <- c("GreyP", "KentishP", "SandP", "Dunlin", "Stint", "Knots", "Curlew", "Godwit", "Greenshanks", "Marsh", "Terek", "Turnstone", "Tattler") accumulated_feedingYK <- data.frame() for (species in spp) { for (n in 1:(length(feeding_time))) { for (i in 1:9) { feeding_time[i,n] <- auc(wader[wader["plot"]==i&wader["transect"]==colnames(feeding_time)[n], "TimeSinceHighTide"], wader[wader["plot"]==i&wader["transect"]==colnames(feeding_time)[n], species]) } #store it in the summary form accumulated_feedingYK[((n-1)*9+1):(9*n), species] <- feeding_time[,n] } accumulated_feedingYK[19:27, species] <- accumulated_feedingYK[1:9, species] + accumulated_feedingYK[10:18, species] } wader <- wader_NP; plottitle <- "NP" feeding_time <- data.frame(NPT1 = c(1:9), NPT2 = c(1:9), NPT3 = c(1:9), NPT4 = c(1:9)) spp <- c("GreyP", "KentishP", "Dunlin", "Stint", "Sanderling", "GreatK", "CurlewS", "RedK", "Curlew", "Godwit", "Greenshanks", "Terek", "Turnstone") accumulated_feedingNP <- data.frame() for (species in spp) { for (n in 1:(length(feeding_time))) { for (i in 1:9) { feeding_time[i,n] <- auc(wader[wader["plot"]==i&wader["transect"]==colnames(feeding_time)[n], "TimeSinceHighTide"], wader[wader["plot"]==i&wader["transect"]==colnames(feeding_time)[n], species]) } #store it in the summary form accumulated_feedingNP[((n-1)*9+1):(9*n), species] <- feeding_time[,n] } accumulated_feedingNP[37:45, species] <- accumulated_feedingNP[1:9, species] + accumulated_feedingNP[10:18, species] + accumulated_feedingNP[19:27, species] + accumulated_feedingNP[28:36, species] } #6 plot the proportion accumulated feeding time and the loss of total proportion to habitat loss #---- for (species in sppYK) { sum_feed <- sum(accumulated_feedingYK[19:27, species]) accumulated_feedingYK[28:36, species] <- accumulated_feedingYK[19:27, species]/sum_feed for(i in 1:9) { accumulated_feedingYK[(36+i), species] <- sum(accumulated_feedingYK[(27+i):36, species]) } tempdata <- data.frame(x = plot_distance$Yangkou, each = accumulated_feedingYK[28:36, species], total = accumulated_feedingYK[37:45, species]) tempdata_melted <- reshape::melt(tempdata, id = "x") if (spp_nameYK[2, species] == "G") {forage_type <- "generalist"; colorcode <- "#E69F00"} if (spp_nameYK[2, species] == "S") {forage_type <- "zone specialist"; colorcode <- "#009E73"} if (spp_nameYK[2, species] == "T") {forage_type <- "tide follower"; colorcode <- "#56B4E9"} if (spp_nameYK[2, species] == "N") {forage_type <- "not assigned"; colorcode <- "#999999"} pdf(paste0(species, "_YK_accumulated.pdf"), width = 6, height = 5) # plot(plot_distance$Yangkou, accumulated_feedingYK[37:45, species], main = paste(species, "YK"), # xlab = "distance from sea wall", ylab = "proportion of cumulative feeding time", # col = "gray") # lines(plot_distance$Yangkou, accumulated_feedingYK[37:45, species], col = "gray") # lines(plot_distance$Yangkou, accumulated_feedingYK[28:36, species]) p1 <- ggplot(tempdata_melted, aes(x = x, y = value, colour = variable, linetype = variable)) + labs(title = paste0(spp_nameYK[1, species], " at Rudong: ", forage_type), x = "distance from sea wall", y = "proportion of cumulative foraging time") + theme_classic() + geom_line() + scale_color_manual(values = c(colorcode, colorcode)) + scale_linetype_manual(values = c("solid", "dashed")) plot(p1) dev.off() } write.csv(accumulated_feedingYK, file = "accumulated_feedingYK.csv") for (species in sppNP) { sum_feed <- sum(accumulated_feedingNP[37:45, species]) accumulated_feedingNP[46:54, species] <- accumulated_feedingNP[37:45, species]/sum_feed for(i in 1:9) { accumulated_feedingNP[(54+i), species] <- sum(accumulated_feedingNP[(45+i):54, species]) } tempdata <- data.frame(x = plot_distance$Nanpu, each = accumulated_feedingNP[46:54, species], total = accumulated_feedingNP[55:63, species]) tempdata_melted <- reshape::melt(tempdata, id = "x") if (spp_nameNP[2, species] == "G") {forage_type <- "generalist"; colorcode <- "#E69F00"} if (spp_nameNP[2, species] == "S") {forage_type <- "zone specialist"; colorcode <- "#009E73"} if (spp_nameNP[2, species] == "T") {forage_type <- "tide follower"; colorcode <- "#56B4E9"} if (spp_nameNP[2, species] == "N") {forage_type <- "not assigned"; colorcode <- "#999999"} pdf(paste0(species, "_NP_accumulated.pdf"), width = 6, height = 5) # plot(plot_distance$Nanpu, accumulated_feedingNP[55:63, species], main = paste(species, "NP"), # xlab = "distance from sea wall", ylab = "proportion of cumulative feeding time", # col = "gray") # lines(plot_distance$Nanpu, accumulated_feedingNP[55:63, species], col = "gray") # lines(plot_distance$Nanpu, accumulated_feedingNP[46:54, species]) p1 <- ggplot(tempdata_melted, aes(x = x, y = value, colour = variable, linetype = variable)) + labs(title = paste0(spp_nameNP[1, species], " at Nanpu: ", forage_type), x = "distance from sea wall", y = "proportion of cumulative foraging time") + theme_classic() + geom_line() + scale_color_manual(values = c(colorcode, colorcode)) + scale_linetype_manual(values = c("solid", "dashed")) plot(p1) dev.off() } write.csv(accumulated_feedingNP, file = "accumulated_feedingNP.csv")