OHI Science | Citation policy
nutrient <- list.files(file.path(dir_M, "git-annex/impact_acceleration/stressors/land_based/int"),
pattern="fert", full=TRUE)
nutrient <- grep("raw", nutrient, value=TRUE)
nutrient <- grep("2002", nutrient, value=TRUE, invert=TRUE) # only want 10 years of data
stack_nutrient <- stack(c(zones_all, nutrient))nutrient_trend <- read.csv(file.path(dir_M, "git-annex/globalprep/hab_seagrass/v2019/sg_trend/seagrass_cell_nutrients.csv"))
21537/215951 # about 10% NA values...more than I would have expected. Might want to gapfill these.
sum(duplicated(nutrient_trend$ID)) # is a unique ID
sum(duplicated(nutrient_trend$cells)) # also is a unique ID
nutrient_trend <- nutrient_trend %>%
rename(rgn_id = layer) %>%
gather("year", "fert", starts_with("global"))
nutrient_trend <- nutrient_trend %>%
mutate(year = gsub("global_plumes_fert_", "", year)) %>%
mutate(year = gsub("_raw", "", year)) %>%
mutate(year = as.numeric(year))
nutrient_lm <- nutrient_trend %>%
filter(!is.na(fert)) %>%
group_by(ID, cells, rgn_id) %>%
do(trend_fert = lm(fert ~ year, data = ., na.action=na.exclude))
# get the coefficients by group in a tidy data_frame
nutrient_coef = tidy(nutrient_lm, trend_fert)
nutrient_coef2 <- nutrient_coef %>%
dplyr::filter(term == "year") %>%
mutate(inverse_estimate = estimate*(-1)) %>%
dplyr::filter(inverse_estimate>-20 & inverse_estimate<30) %>%
data.frame()
ggplot(dplyr::filter(nutrient_coef2, estimate <= 1 & estimate >-1), aes(x=estimate)) +
geom_histogram()From 1990-2000 (the most recent years of trend analysis) there were 115 sites with decreasing seagrass, 72 sites with increasing seagrass and 88 sites with no detectable change. I will assume the proportion of seagrass beds sampled in each category is representative of global seagrass beds.
I also assume that each ~1km2 raster cell is an independent seagrass bed. This is similar to the average size based on polygon areas in original data.
# proportion sites decreasing:
115/(115+72+88) # 0.42
# proportion sites increasing:
72/(115+72+88) # 0.26
# proportion sites remaining the same:
88/(115+72+88) # 0.32
# find values that correspond to these quantiles in the raster data
quants <- quantile(nutrient_coef2$inverse_estimate, c(0.42, 0.74))
nutrient_coef2$trend_cat <- cut(nutrient_coef2$inverse_estimate, c(-Inf, quants[1], quants[2], Inf),
labels=c("Decreasing", "Constant", "Increasing"))
# check I did this correctly:
table(nutrient_coef2$trend_cat)
81651/(81651+62210+50546)The constant values will have a trend of 0. For the Increasing and Decreasing categories this will scale linearly from 0 to the 99th predicted quantile based on SD.
library("scales")
decreasing_quant <- quantile(nutrient_coef2$inverse_estimate[nutrient_coef2$trend_cat=="Decreasing"], c(0.5))
nutrient_coef2$rescale_0to1 <- ifelse(nutrient_coef2$trend_cat=="Decreasing",
rescale(nutrient_coef2$inverse_estimate[nutrient_coef2$trend_cat=="Decreasing"]), NA)
sd_sample <- 1.92*sqrt(265)
hist(rnorm(1000, mean=-1.26, sd=sd_sample))
install.packages("fGarch")
library(fGarch)
test <- rsnorm(100000, mean=-1.26, sd = sd_sample, xi=-1.1)
hist(test)
mean(test)
median(test)