Merge pull request #51 from emlab-ucsb/dev-gfw

Dev gfw and geomorphology

DESCRIPTION

@@ -19,10 +19,11 @@ Imports:
     stats,
     spatialgridr (>= 0.0.2.0),
     rlang
-Remotes: 
+Remotes:  
     github::emlab-ucsb/spatialgridr,
     github::bio-oracle/biooracler,
-    github::lifewatch/mregions2
+    github::lifewatch/mregions2,
+    github::GlobalFishingWatch/gfwr
 RoxygenNote: 7.3.2
 Depends: 
     R (>= 3.5.0)
@@ -34,7 +35,8 @@ Suggests:
     parallel,
     remotes,
     biooracler (>= 0.0.0.9000),
-    mregions2 (>= 1.0.0)
+    mregions2 (>= 1.0.0),
+    gfwr (>= 2.0.0)
 Config/testthat/edition: 3
 VignetteBuilder: knitr
 URL: https://emlab-ucsb.github.io/oceandatr/

NAMESPACE

@@ -10,6 +10,7 @@ export(get_ecoregion)
 export(get_enviro_regions)
 export(get_features)
 export(get_geomorphology)
+export(get_gfw)
 export(get_grid)
 export(get_knolls)
 export(get_seamounts)

R/get_enviro_regions.R

@@ -201,9 +201,17 @@ get_enviro_data <- function(spatial_grid = NULL){
 
   grid_bbox <- sf::st_bbox(polygon4326)
 
+  #queries to the ERDDAP server where Bio-Oracle is hosted only allow coordinates -89.975 to 89.975, and -179.975 to 179.975, e.g: https://erddap.bio-oracle.org/erddap/info/tas_baseline_2000_2020_depthsurf/index.html
+  x_min <- if(grid_bbox["xmin"][[1]] < -179.975) -179.975 else grid_bbox["xmin"][[1]]
+  x_max <- if(grid_bbox["xmax"][[1]] > 179.975) 179.975 else grid_bbox["xmax"][[1]]
+
+  y_min <- if(grid_bbox["ymin"][[1]] < -89.975) -89.975 else grid_bbox["ymin"][[1]]
+  y_max <- if(grid_bbox["ymax"][[1]] > 89.975) 89.975 else grid_bbox["ymax"][[1]]
+
+
   constraints <- list(time = c('2010-01-01T00:00:00Z', '2010-01-01T00:00:00Z'),
-                      latitude = c(as.numeric(grid_bbox["ymin"]), as.numeric(grid_bbox["ymax"])),
-                      longitude = c(as.numeric(grid_bbox["xmin"]), as.numeric(grid_bbox["xmax"])))
+                      latitude = c(y_min, y_max),
+                      longitude = c(x_min, x_max))
 
   biooracle_data <- list()
 

R/get_geomorphology.R

@@ -2,9 +2,11 @@
 #' 
 #' @description Get geomorphological data for a spatial grid or polygon
 #' 
-#' @details Geomorphological features are from the [Harris et al. 2014](https://doi.org/10.1016/j.margeo.2014.01.011) dataset, available at [https://www.bluehabitats.org](https://www.bluehabitats.org), excluding depth related features which can be created using `get_bathymetry()`:
+#' @details Geomorphological features are from the [Harris et al. 2014](https://doi.org/10.1016/j.margeo.2014.01.011) dataset, available at [https://www.bluehabitats.org](https://www.bluehabitats.org). Data is included in this package, except depth classification features which can be created using `get_bathymetry()` and seamounts which can be retrieved from a more recent dataset using `get_seamounts()`. List of features:
 #' 
 #' \itemize{
+#' \item Abyssal hills
+#' \item Abyssal plains
 #' \item Basins: 
 #'    \itemize{
 #'    \item large basins of seas and oceans
@@ -63,12 +65,14 @@ get_geomorphology <- function(spatial_grid = NULL, raw = FALSE, antimeridian = N
   meth <- if(check_raster(spatial_grid)) 'near' else 'mode'
 
   sf::sf_use_s2(FALSE)
-  geomorph_data <- system.file("extdata/geomorphology", package = "oceandatr") %>% 
-    list.files() %>% 
-    system.file("extdata/geomorphology", ., package = "oceandatr") %>% 
-    lapply(readRDS) %>% 
-    do.call(rbind, .) %>%
-    get_data_in_grid(spatial_grid = spatial_grid, dat = ., raw = raw, meth = meth, feature_names = "geomorph_type", antimeridian = antimeridian)
+  suppressWarnings(
+    geomorph_data <- system.file("extdata/geomorphology", package = "oceandatr") %>% 
+      list.files() %>% 
+      system.file("extdata/geomorphology", ., package = "oceandatr") %>% 
+      lapply(readRDS) %>% 
+      do.call(rbind, .) %>%
+      get_data_in_grid(spatial_grid = spatial_grid, dat = ., raw = raw, meth = meth, feature_names = "geomorph_type", antimeridian = antimeridian)
+  )
 
   sf::sf_use_s2(TRUE)
 

R/get_gfw.R

@@ -0,0 +1,123 @@
+#' Get data from Global Fishing Watch
+#'
+#' @description Global Fishing Watch (GFW) provides data on apparent fishing effort (hours) at 0.01 degree spatial resolution, based on automatic identification system (AIS) broadcasts from vessels. Data is principally for larger vessel (> 24m in length); less than 1% of vessels <12 m length are represented in the data (see [GFW website](https://globalfishingwatch.org/dataset-and-code-fishing-effort/) for detailed information). This function is primarily a wrapper for the [`gfwr` package](https://github.com/GlobalFishingWatch/gfwr) function `get_raster()`, but allows the user to return multiple years of data in a summarized and gridded format. An API key is required to retrieve GFW data; see the package website for instructions on how to get and save one (free).
+#' 
+#' @param spatial_grid `sf` or `terra::rast()` grid, e.g. created using `get_grid()`. Alternatively, if raw data is required, an `sf` polygon can be provided, e.g. created using `get_boundary()`, and set `raw = TRUE`.
+#' @param raw `logical` if TRUE, `spatial_grid` can be an `sf` polygon, and the raw GFW data, in tibble format, is returned for a bounding box covering the polygon +/-1 degree. The bounding box is expanded to ensure that the entire polygon has data coverage once the point data is rasterized. This data will not be summarised, i.e. `summarise` is ignored.
+#' @param resolution `string` either `"HIGH"` = 0.01 degree spatial resolution, or `"LOW"`  = 0.1.
+#' @param start_year `numeric` must be 2012 or more recent year. Note that GFW added data from Spire and Orbcomm AIS providers in 2017, so data from 2017 is likely to have greater spatial and temporal coverage ([Welch et al. 2022](https://doi.org/10.1126/sciadv.abq2109)). 
+#' @param end_year `numeric` any year between 2012 and the current year
+#' @param group_by `string` can be `"geartype"`, `"flag"`, or `"location"`. 
+#' @param summarise `string` can be `"total_annual_effort"`; the sum of all fishing effort for all years, from `start_year` to `end_year` at each location is returned, grouped by `group_by`; or `"mean_total_annual_effort"`; the mean of the annual sums of fishing effort for all years at each location is returned, grouped by `group_by`.
+#' @param key `string` Authorization token. Can be obtained with gfw_auth() function. See `gfwr` [website](https://github.com/GlobalFishingWatch/gfwr?tab=readme-ov-file#authorization) for details on how to request a token.
+#'
+#' @return For gridded data, a `terra::rast()` or `sf` object, depending on the `spatial_grid` format. If `raw = TRUE`, non-summarised data in `tibble` format is returned for the polygon area direct from the GFW query `gfwr::get_raster()`.
+#' @export
+#'
+#' @examplesIf nchar(Sys.getenv("GFW_TOKEN"))>0
+#' #get mean total annual fishing effort for Bermuda for the years 2022-2023
+#' #first get a grid for Bermuda
+#' bermuda_grid <- get_boundary(name = "Bermuda") %>% get_grid(resolution = 0.1, crs = 4326)
+#' 
+#' bermuda_gfw_effort <- get_gfw(spatial_grid = bermuda_grid, start_year = 2022)
+#' 
+#' #plot the data
+#' terra::plot(bermuda_gfw_effort)
+#' 
+#' #get total fishing effort for each gear type in Fiji's EEZ for 2022
+#' fiji_grid <- get_boundary(name = "Fiji") %>% get_grid(resolution = 1e4, crs = "+proj=tcea +lon_0=178 +datum=WGS84 +units=m +no_defs", output = "sf_square")
+#' 
+#' fiji_gfw_effort <- get_gfw(spatial_grid = fiji_grid, start_year = 2022, end_year = 2022, group_by = "geartype", summarise = "total_annual_effort")
+#' 
+#' plot(fiji_gfw_effort, border = FALSE)
+#' 
+#' #quantile is better for viewing the fishing effort distribution due to the long tail of values
+#' plot(fiji_gfw_effort[1], border= FALSE, breaks = "quantile")
+get_gfw <- function(spatial_grid = NULL, raw = FALSE, resolution = "LOW", start_year = 2018, end_year = 2023, group_by = "location", summarise = "mean_total_annual_effort", key = gfwr::gfw_auth()){
+
+  current_year <- as.numeric(format(Sys.Date(), "%Y"))
+
+  if(end_year > current_year) stop('"end_year" must be ', current_year, " or before.")
+
+  if(start_year < 2012) stop("The first available year of Global Fishing Watch data is 2012")
+
+  number_years <- end_year - start_year + 1
+
+  if(group_by == "location") gfw_group <- "GEARTYPE" else gfw_group <- toupper(group_by)
+
+  gfw_shape <- spatial_grid %>%
+    polygon_in_4326() %>%
+    sf::st_buffer(dist = 0.5) %>% 
+    sf::st_wrap_dateline()
+
+  #create unique filename for caching - assumes that bounding box can be used to distinguish different spatial queries
+  bounding_box <- sf::st_bbox(spatial_grid)
+  file_name <- paste("gfwdata_x1", bounding_box[[1]], "y1", bounding_box[[2]], "x2", bounding_box[[3]], "y2", bounding_box[[4]], gfw_group, start_year, end_year, resolution, ".rds", sep = "_")
+
+  if(file_name %in% list.files(path = tempdir())) {
+    message("GFW data already downloaded, using cached version", 
+            sep = "")
+    fishing_effort <- readRDS(file.path(tempdir(), file_name))
+  } else{
+    fishing_effort <- lapply(start_year:end_year, function(yr){
+      gfwr::get_raster(spatial_resolution = resolution,
+                       temporal_resolution = 'YEARLY',
+                       group_by = gfw_group,
+                       start_date = paste0(yr, "-01-01"),
+                       end_date = paste0(yr, "-12-31"),
+                       region = gfw_shape,
+                       region_source = 'USER_SHAPEFILE',
+                       key = key)
+    }) 
+    fishing_effort <- do.call(dplyr::bind_rows, fishing_effort)
+
+    saveRDS(fishing_effort, file = file.path(tempdir(), file_name))
+  }
+
+  if(raw){
+    return(fishing_effort)
+  } 
+
+    grouping_vars <- c("Lon", "Lat", "Time Range") %>% 
+      {if(group_by == "location") . else c(., tolower(group_by))} 
+
+    annual_effort <- fishing_effort %>% 
+      dplyr::group_by(dplyr::across(dplyr::all_of(grouping_vars))) %>% 
+      dplyr::summarise(total_annual_effort = sum(`Apparent Fishing Hours`, na.rm = TRUE)) %>% 
+      dplyr::ungroup()  
+
+    if(summarise == "total_annual_effort"){
+      if(group_by == "location"){
+        final_effort <- annual_effort %>% 
+          tidyr::pivot_wider(names_from = "Time Range",
+                             values_from = total_annual_effort) 
+
+      }else{
+        final_effort <- annual_effort %>% 
+          tidyr::pivot_wider(names_from = dplyr::all_of(c(group_by, "Time Range")),
+                             values_from = total_annual_effort) 
+      }
+    } else{
+      mean_total_effort <- annual_effort %>% 
+        dplyr::group_by(., dplyr::across(-c("Time Range", total_annual_effort))) %>%
+        dplyr::summarise(mean_total_annual_effort = sum(total_annual_effort, na.rm = TRUE)/number_years) %>% #calculate mean manually to ensure that years that have NA catch in a cell are still included in the denominator
+        dplyr::ungroup()
+
+      if(group_by == "location"){
+        final_effort <- mean_total_effort
+
+      }else{
+        final_effort <- mean_total_effort %>% 
+          tidyr::pivot_wider(names_from = dplyr::all_of(group_by),
+                             values_from = mean_total_annual_effort)
+      }
+    }
+
+    final_effort %>% 
+      terra::rast(type = "xyz", crs = "epsg:4326") %>% 
+        # terra::subst(NA, 0.01) %>% #too many zero cost value cells leads to prioritization of too much area because they are 'free'
+          get_data_in_grid(spatial_grid = spatial_grid,
+                           dat = ., 
+                           meth = if(check_sf(spatial_grid)) "mean" else "average")
+
+}

R/utils.R

@@ -134,6 +134,13 @@ classify_layers <- function(dat, dat_breaks = NULL, classification_names = NULL)
       dplyr::select(3:ncol(.), 2) #put classification before geometry and drop original values
   }
 }
+#' Get an sf polygons in lonlat (EPSG 4326) from terra or sf input object
+#'
+#' @param spatial_grid 
+#'
+#' @return `sf` polygons
+#'
+#' @noRd
 polygon_in_4326 <-
   function(spatial_grid) {
     crs_is_4326 <- check_matching_crs(spatial_grid, 4326)
@@ -147,15 +154,17 @@ polygon_in_4326 <-
               terra::project(., "epsg:4326")
           } %>%
           sf::st_as_sf() %>% 
+          sf::st_union() %>%
+          sf::st_sf() %>% 
+          sf::st_wrap_dateline() %>% 
           {if(sf::st_is_valid(.)) . else sf::st_make_valid(.)}
       } else{
         spatial_grid %>%
-          {
-            if (crs_is_4326)
-              .
-            else
-              sf::st_transform(., 4326)
-          }
+          {if (crs_is_4326) . else sf::st_transform(., 4326)} %>% 
+          sf::st_union() %>% 
+          sf::st_sf() %>% 
+          sf::st_wrap_dateline() %>% 
+          {if(sf::st_is_valid(.)) . else sf::st_make_valid(.)}
       }
   }
 

README.Rmd

@@ -127,10 +127,11 @@ terra::plot(depth_zones, col = c("grey60", "navyblue"), axes = FALSE, fun = func
 
 ## Get geomorphological data
 
-The seafloor has its own mountains, plains and other geomorphological features just as on land. These data come from [Harris et al. 2014, Geomorphology of the Oceans](https://doi.org/10.1016/j.margeo.2014.01.011) and are available for download from https://www.bluehabitats.org. The features that are suggested as major habitats for inclusion in no-take MPAs by [Ceccarelli et al. 2021](https://doi.org/10.3389/fmars.2021.634574) are included in this package, so it is not necessary to download them.
+The seafloor has its own valleys, plains and other geomorphological features just as on land. These data come from [Harris et al. 2014, Geomorphology of the Oceans](https://doi.org/10.1016/j.margeo.2014.01.011) and are available for download from https://www.bluehabitats.org, but all but depth classifications (which can be created using `get_bathymetry()`) and seamounts (which can be retrieved from a more recent dataset using `get_seamounts()`) are included in this package. 
 
 ```{r geomorphology, warning=FALSE, message=FALSE}
-geomorphology <- get_geomorphology(spatial_grid = bermuda_grid)
+geomorphology <- get_geomorphology(spatial_grid = bermuda_grid) %>% 
+  remove_empty_layers() #can remove any empty layers so we don't have so many layers to plot
 
 #brown colour indicates that geomorphological feature is present
 terra::plot(geomorphology, col = data.frame(c(0,1), c("grey60", "sienna")), axes = FALSE, legend = FALSE, fun = function(){terra::lines(terra::vect(bermuda_eez_projected))})
@@ -206,18 +207,3 @@ enviro_regions <- get_enviro_regions(spatial_grid = bermuda_grid, show_plots = T
 #value of 1 indicates that environmental region is present
 terra::plot(enviro_regions, col = c("grey60", "forestgreen"), axes = FALSE, fun = function(){terra::lines(terra::vect(bermuda_eez_projected))})
 ```
-
-
-```{r pressure, echo = FALSE, include=FALSE}
-## REMINDER FOR ME, WILL DELETE:
-
-# You'll still need to render `README.Rmd` regularly, to keep `README.md` up-to-date. `devtools::build_readme()` is handy for this. You could also use GitHub Actions to re-render `README.Rmd` every time you push. An example workflow can be found here: <https://github.com/r-lib/actions/tree/v1/examples>.
-# 
-# You can also embed plots, for example:
-
-#plot(pressure)
-
-#In that case, don't forget to commit and push the resulting figure files, so they display on GitHub and CRAN.
-```
-
-

README.md

@@ -124,9 +124,6 @@ package.
 ``` r
 bathymetry <- get_bathymetry(spatial_grid = bermuda_grid, classify_bathymetry = FALSE)
 #> This may take seconds to minutes, depending on grid size
-```
-
-``` r
 
 terra::plot(bathymetry, col = hcl.colors(n=255, "Blues"), axes = FALSE) 
 plot(bermuda_eez_projected, add=TRUE)
@@ -150,9 +147,6 @@ We can get the depth zones for Bermuda simply by setting the
 ``` r
 depth_zones <- get_bathymetry(spatial_grid = bermuda_grid, classify_bathymetry = TRUE)
 #> This may take seconds to minutes, depending on grid size
-```
-
-``` r
 
 #value of 1 indicates that depth zone is present
 terra::plot(depth_zones, col = c("grey60", "navyblue"), axes = FALSE, fun = function(){terra::lines(terra::vect(bermuda_eez_projected))})
@@ -162,17 +156,18 @@ terra::plot(depth_zones, col = c("grey60", "navyblue"), axes = FALSE, fun = func
 
 ## Get geomorphological data
 
-The seafloor has its own mountains, plains and other geomorphological
+The seafloor has its own valleys, plains and other geomorphological
 features just as on land. These data come from [Harris et al. 2014,
 Geomorphology of the
 Oceans](https://doi.org/10.1016/j.margeo.2014.01.011) and are available
-for download from <https://www.bluehabitats.org>. The features that are
-suggested as major habitats for inclusion in no-take MPAs by [Ceccarelli
-et al. 2021](https://doi.org/10.3389/fmars.2021.634574) are included in
-this package, so it is not necessary to download them.
+for download from <https://www.bluehabitats.org>, but all but depth
+classifications (which can be created using `get_bathymetry()`) and
+seamounts (which can be retrieved from a more recent dataset using
+`get_seamounts()`) are included in this package.
 
 ``` r
-geomorphology <- get_geomorphology(spatial_grid = bermuda_grid)
+geomorphology <- get_geomorphology(spatial_grid = bermuda_grid) %>% 
+  remove_empty_layers() #can remove any empty layers so we don't have so many layers to plot
 
 #brown colour indicates that geomorphological feature is present
 terra::plot(geomorphology, col = data.frame(c(0,1), c("grey60", "sienna")), axes = FALSE, legend = FALSE, fun = function(){terra::lines(terra::vect(bermuda_eez_projected))})

data-raw/geomorphology_data_prep.R

@@ -10,12 +10,11 @@ library(dplyr)
 data_file_path <- "temp_raw/global-seafloor-geomorphology"
 
 geomorph_files <- list.files(path = data_file_path, full.names = TRUE, pattern = '\\.shp$') %>% 
-  #Remove features not suggested for prioritization in Ceccarelli et al. 2021 Table 5. Want to include Shelf_valleys, hence removal of specific shelf strings
-  {.[-grep("Abyss|Hadal|Seamounts|Shelf_Classification|Shelf\\.|Slope", .)]}
+  #Remove features that are just classification of depths - will do this using latest GEBCO data. Want to include Shelf_valleys, hence removal of specific shelf strings. Removal classification of abyssal depths ("Abyss") but want to keep "Abyssal_Classification", hence the specific string.
+  {.[-grep("Abyss\\.|Hadal|Seamounts|Shelf_Classification|Shelf\\.|Slope", .)]}
 
 sf_use_s2(FALSE)
 
-
 for (file_name in geomorph_files) {
   feature_name <- gsub(pattern =  ".shp",replacement =  "", basename(file_name))
 
@@ -40,6 +39,15 @@ for (file_name in geomorph_files) {
         # 
         # do.call("use_data", list(as.name(paste0(feature_name, "_", geomorph_type)), overwrite = TRUE))
     } 
+  } else if(feature_name == "Abyssal_Classification"){
+    for(abyssal_class in c("Hills", "Plains")){ #only want Hill and Plains, not seamounts since these will come from more recent data
+      geomorph_sf_object %>% 
+        filter(class == abyssal_class) %>% 
+        st_union() %>% 
+        st_sf %>% 
+        dplyr::mutate(geomorph_type = paste0("Abyssal_", abyssal_class), .before = 1) %>% 
+        saveRDS(file = file.path("inst/extdata/geomorphology", paste0("Abssyal_", abyssal_class, ".rds")))
+    }
   } else{
     geomorph_sf_object %>%  
       st_union() %>%