Coral trend was calculated using condition data from 1975-2006 (Bruno and Selig 2007, Schutte et al. 2010).
Using a new dataset for the global distribution of coral reefs, so we are going to gapfill the trend of the new regions that have been added to the coral extent data.
## This file makes it easier to process data for the OHI global assessment
## by creating the following objects:
##
## * dir_M = identifies correct file path to Mazu (internal server) based on your operating system
## * mollCRS = the crs code for the mollweide coordinate reference system we use in the global assessment
## * regions_shape() = function to load global shapefile for land/eez/high seas/antarctica regions
## * ohi_rasters() = function to load two rasters: global eez regions and ocean region
## * region_data() = function to load 2 dataframes describing global regions
## * rgn_syns() = function to load dataframe of region synonyms (used to convert country names to OHI regions)
## * low_pop() = function to load dataframe of regions with low and no human population
## * UNgeorgn = function to load dataframe of UN geopolitical designations used to gapfill missing data
## loads 2 dataframes: rgns_all and rgns_eez
## rgns_all = includes eez/high seas/antarctica regions, IDs correspond with region shapefile and raster
## rgns_eez = includes only eez regions
Now match the old trend data to the extent data, and gapfill any missing trend data with georegion averages. This will also exclude any regions that have trend data but no extent data.
<- read_csv(file.path(dir_git_old, "data/habitat_trend_coral_updated.csv"))
trend
## read in georegions data
<- read.csv(file.path(here('globalprep/supplementary_information/v2018/rgn_georegions_wide_2013b.csv')))
regions
## read in extent data
<- read.csv(file.path(here("globalprep/hab_coral/v2021/data/habitat_extent_coral_updated.csv")))
all
<- read.csv(file.path(here("globalprep/hab_coral/v2012/data/habitat_extent_coral_updated.csv")))
all_old
<- all %>%
all filter(km2 > 0) %>%
filter(rgn_id < 255) %>%
::select(-year)
dplyr
setdiff(all$rgn_id, trend$rgn_id) # coral extent data but no trend score # 6 10 11 12 28 38 44 45 46 48 49 50 51 54 79 107 114 119 125 137 146 148 149 150 155 157 159 190 191 204 212 231 237 136
setdiff(trend$rgn_id, all$rgn_id) # have trend but no extent... two are small islands. One is new zealand. New Zealand does not have a coral reef.. so probably not a big deal.
# [1] 3 90 162
<- all %>% #eliminates the regions with health scores but no mangroves
tmp left_join(trend)
<- tmp %>%
tmp left_join(regions) %>%
group_by(r2) %>%
mutate(avg_trend_r2 = mean(trend, na.rm=TRUE)) %>%
ungroup() %>%
group_by(r1) %>%
mutate(avg_trend_r1 = mean(trend, na.rm=TRUE)) %>%
ungroup() %>%
mutate(trend2 = ifelse(is.na(trend), avg_trend_r2, trend)) %>%
mutate(trend3 = ifelse(is.na(trend2), avg_health_r1, trend2)) %>%
mutate(gap_fill = ifelse(is.na(trend), "georegion", NA)) %>%
mutate(gap_fill = ifelse(is.na(trend) & is.na(avg_trend_r2), "georegion", gap_fill)) %>%
mutate(gap_fill = ifelse(is.na(gap_fill), 0, gap_fill)) %>%
::select(rgn_id, habitat, gap_fill, trend=trend3)
dplyrsummary(tmp)
### save summary of gap-filling:
<- read.csv(file.path(here("globalprep/hab_coral/v2012/data/trend_coral_gf.csv")))
old_gaps
<- tmp %>%
trend_gaps left_join(old_gaps, by = c("rgn_id", "habitat")) %>%
mutate(gap_fill = ifelse(gap_fill.x != 0, gap_fill.x, gap_fill.y)) %>%
mutate(gap_fill = ifelse(is.na(gap_fill), 0, gap_fill)) %>%
full_join(rgns_eez) %>%
::select(rgn_id, habitat, variable, gap_fill) %>%
dplyrmutate(variable = "trend") %>%
mutate(habitat = "coral")
write.csv(trend_gaps, 'globalprep/hab_coral/v2021/data/trend_coral_gf.csv', row.names=FALSE)
### save health data:
<- tmp %>%
trend ::select(rgn_id, habitat, trend) %>%
dplyrmutate(year = 2021)
write.csv(trend, 'globalprep/hab_coral/v2021/data/habitat_trend_coral_updated.csv', row.names=FALSE)
Datacheck: compare to old health data
<- read_csv(file.path(here("globalprep/hab_coral/v2012/data/habitat_trend_coral_updated.csv")))
old_trend
<- trend %>%
compare_habitat_trend filter(rgn_id <= 250) %>%
left_join(old_trend, by = "rgn_id") %>%
mutate(difference = trend.x - trend.y) %>%
left_join(rgns_eez)
<- compare_habitat_trend %>%
test filter(is.na(trend.y)) ## we gained 34 regions with trend data
ggplot(compare_habitat_trend, aes(x = trend.y, y = trend.x)) +
geom_point() +
geom_abline() +
labs(title = "mangrove Habitat version old (v2015) vs new (v2021)", x = "old trend", y= "new trend") +
theme_bw()