OHI Science | Citation policy
This script generates the extent of kelp for each OHI region for the latest year of data. We do this using polygon data. We extract the area of each polygon per each ohi region using our combined land and eez raster.
This is an entirely new layer for the 2021 assessment!
Reference: Jayathilake, D.R.M., Costello, M.J., 2020. A modelled global distribution of the kelp biome. Biological Conservation 252, 108815. https://doi.org/10.1016/j.biocon.2020.108815
Downloaded: 04/28/2021
Description:
This is a MaxEnt model map of the global distribution of the Laminarian kelp biome. This was modelled using 44,265 cleaned primary occurrence records from the Global Biodiversity Information Facility (GBIF), and the Ocean Biodiversity Information System (OBIS) and 13 environmental var
Time range: 1900 - 2020
Read in kelp data and explore
kelp_shp <- st_read(file.path(dir_kelp, "01_Data/Kelp UNEP.shp"))
mapview(head(kelp_shp))Match the data to OHI regions
kelp_moll <- st_transform(kelp_shp, crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
fasterize_py <- fasterize::fasterize(kelp_moll, raster = land_ocean,
field = NULL) # all polygons given value of 1
fasterize_py
plot(fasterize_py)
print(cellStats(fasterize_py, "sum", na.rm=TRUE)*0.934478877011218970**2) # 1,403,760 km2 - paper says 1,469,900 km2
sum(kelp_moll$Shape_Area) # 1.469908e+12 - we are losing about 60000 km2... this seems reasonable though.
writeRaster(fasterize_py, filename = file.path(dir_M, "/git-annex/globalprep/hab_kelp/v2021/int/py_kelp_rast.tif"), overwrite=TRUE)
# Calculate zonal stats with zones raster and new combined kelp. Convert to km^2 and save int/output files*
zonal_sums_combined <- raster::zonal(fasterize_py,
land_ocean,
fun = "sum",
na.rm = TRUE,
progress="text") #sum all kelp cells for each ohi zone
zonal_sums_combined_df <- data.frame(zonal_sums_combined)
zonal_sums_km2 <- zonal_sums_combined_df %>%
mutate(year = 2020, habitat = "kelp",
km2 = (0.934478877011218970**2*sum)) %>% #one cell is equal to ~0.93 km
rename("rgn_id" = "zone") %>%
dplyr::select(-sum) %>%
filter(rgn_id <= 250)
print(sum(zonal_sums_km2$km2)) # 1394773
write.csv(zonal_sums_km2, file.path(here("globalprep/hab_kelp/v2021/data/habitat_extent_kelp.csv")), row.names = FALSE)
## make gapfilling file
kelp_extent_gf <- read_csv(file.path(dir_git, "data/habitat_extent_kelp.csv")) %>%
mutate(variable = "extent") %>%
dplyr::select(-km2, -year) %>%
mutate(gap_fill = 0)
write.csv(kelp_extent_gf, file.path(dir_git, "data/extent_kelp_gf.csv"), row.names = FALSE)