Merge pull request #100 from Rafnuss/v3.1.0

V3.1.0

.github/workflows/update-citation-cff.yaml

@@ -0,0 +1,58 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/master/examples
+# The action runs when:
+# - A new release is published
+# - The DESCRIPTION or inst/CITATION are modified
+# - Can be run manually
+# For customizing the triggers, visit https://docs.github.com/en/actions/learn-github-actions/events-that-trigger-workflows
+on:
+  release:
+    types: [published]
+  push:
+    branches: [master, main]
+    paths:
+      - DESCRIPTION
+      - inst/CITATION
+  workflow_dispatch:
+
+name: Update CITATION.cff
+
+jobs:
+  update-citation-cff:
+    runs-on: macos-latest
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+    steps:
+      - uses: actions/checkout@v3
+      - uses: r-lib/actions/setup-r@v2
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: |
+            any::cffr
+            any::V8
+
+      - name: Update CITATION.cff
+        run: |
+
+          library(cffr)
+
+          # Customize with your own code
+          # See https://docs.ropensci.org/cffr/articles/cffr.html
+
+          # Write your own keys
+          mykeys <- list()
+
+          # Create your CITATION.cff file
+          cff_write(keys = mykeys)
+
+        shell: Rscript {0}
+
+      - name: Commit results
+        run: |
+          git config --local user.name "github-actions[bot]"
+          git config --local user.email "41898282+github-actions[bot]@users.noreply.github.com"
+          git add CITATION.cff
+          git commit -m 'Update CITATION.cff' || echo "No changes to commit"
+          git push origin || echo "No changes to commit"
+
+
+

.lintr

@@ -2,5 +2,5 @@ linters: linters_with_defaults(
   line_length_linter(100),
   commented_code_linter = NULL,
   spaces_left_parentheses_linter = NULL,
-  indentation = NULL,
+  indentation_linter = NULL,
   cyclocomp_linter = NULL)

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: GeoPressureR
 Title: Global Positioning by Atmospheric Pressure
-Version: 3.0.0
+Version: 3.1.0
 Authors@R: c(
     person("Raphaël", "Nussbaumer", , "rafnuss@gmail.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-8185-1020")),
@@ -31,9 +31,9 @@ Imports:
     methods,
     ncdf4,
     pracma,
-    progressr,
-    terra,
+    readr,
     stats,
+    terra,
     tools,
     utils
 Suggests: 
@@ -45,9 +45,9 @@ Suggests:
     shiny,
     RColorBrewer,
     readxl,
-    readr,
     shinyWidgets,
-    covr
+    covr,
+    maps
 biocViews:
 Config/testthat/edition: 3
 Encoding: UTF-8

NEWS.md

@@ -1,3 +1,17 @@
+# GeoPressureR v3.1.0
+## Major
+* Update to GeoPressureAPI v2, using `thr_mask` in `geopressure_map_mismatch()` and splitting `keep_mse_mask`.
+* Adjust computation of ground speed to account for grid resolution. 
+
+## Minor
+* Use negative indexing for `known`
+* remove trailing `/` to default directories.
+* documentations and minor fixes.
+
+
+## Full Changelog
+[https://github.com/Rafnuss/GeoPressureR/compare/v3.0.0...v3.1.0](https://github.com/Rafnuss/GeoPressureR/compare/v3.0.0...v3.1.0)
+
 # GeoPressureR v3.0.0
 
 ## Guiding principles of v3

R/geolight_map.R

@@ -33,7 +33,7 @@
 #' @return a `tag` with the likelihood of light as `tag$map_light`
 #'
 #' @examples
-#' setwd(system.file("extdata/", package = "GeoPressureR"))
+#' setwd(system.file("extdata", package = "GeoPressureR"))
 #' # Read geolocator data and build twilight
 #' tag <- tag_create("18LX", quiet = TRUE) |>
 #'   tag_label(quiet = TRUE) |>
@@ -58,7 +58,7 @@
 #' @export
 geolight_map <- function(tag,
                          twl_calib_adjust = 1.4,
-                         twl_llp = \(n) 0.1,
+                         twl_llp = \(n) log(n) / n,
                          compute_known = FALSE,
                          quiet = FALSE) {
   # Check tag

R/geopressure_map.R

@@ -1,71 +1,150 @@
 #' Compute likelihood map from pressure data
 #'
 #' @description
-#' This function computes a likelihood map for each stationary period based on pressure
-#' measurements. It performs the following operations:
+#' The `geopressure_map()` function computes a likelihood map for each stationary period based on
+#' pressure measurements. It is a wrapper of the following two child functions:
 #'
-#' 1. Compute the mismatch map between the pressure sensor measurements and the ERA5 reanalysis
-#' database. See [`geopressure_map_mismatch()`] for details.
-#' 2. Convert the mismatch map into a likelihood map with [`geopressure_map_likelihood()`].
+#' 1. `geopressure_map_mismatch()` computes the mismatch maps between the pressure sensor
+#' measurements and the ERA5 reanalysis database.
+#' 2. `geopressure_map_likelihood()` converts the mismatch maps into a likelihood map.
+#'
+#' See below for details.For more background on the method behind these functions, please refer to
+#' Nussbaumer et al. ([2023a](https://doi.org/10.1111/2041-210X.14043)).
+#'
+#' @details
 #'
 #' A map will only be computed for the stationary periods included in `tag$stap$include` and without
-#' a known position `tag$stap$known_l**`. If the position is known, the function will
-#' create a likelihood map with a single 1 value at the grid cell closest to the known position.
+#' a known position `tag$stap$known_l**` as defined by `tag_set_map()`. At known stationary
+#' periods, the likelihood map is a map of `0`s with a single `1` at the pixel closest to the known
+#' position.
 #'
-#' For more background and details on the method behind these functions, please refer to the
-#' Nussbaumer et al. ([2023a](https://doi.org/10.1111/2041-210X.14043)).
+
+#' @section Mismatch map with GeoPressureAPI:
+#'
+#' `geopressure_map_mismatch()` computes the mismatch maps on Google Earth Engine via the map
+#' entry point of the [GeoPressure API](https://raphaelnussbaumer.com/GeoPressureAPI/#description).
+#' This consists of the following steps:
+#'
+#' 1. **Pre-process pressure**: the pressure measurements are first smoothed and downscaled to a
+#' 1-hour resolution in order to match ERA-5 resolution (see `geopressure_map_preprocess()`).
+#' 2. **Generate requests**: Send a single request to the GeoPressureAPI to generate the Google
+#' Earth Engine (GEE) URLs, one for each elevation/stationary period which can be used to compute
+#' the maps on the GEE server.
+#' 3. **Send the requests**: Call the URLs, which will start the computation on the GEE server.
+#' At this step, it does not wait for an answer, but send the requests in parallel to be faster.
+#' 4. **Compute and download the maps**: When all requests are sent, we wait for the GEE server to
+#' return a geotiff file (map) for each elevation/ stationary period.
+#' 5. **Post-process maps**: Read these geotiff maps as matrix, create the corresponding
+#' GeoPressureR `map` object with `create_map()`.
+#'
+#' The following two maps are returned for each stationary period:
+#'
+#' 1. **map_pressure_mse**: The Mean Square Error (MSE) between the data logger pressure timeseries
+#' and the reanalysis. The mean error is removed because we assume no specific altitude of the
+#' geolocator, thus allowing an altitudinal shift of the pressure timeseries.
+#' 2. **map_pressure_mask** (optionally): The mask of the proportion of pressure measurements
+#' corresponding to altitude values found within the min and max ground elevation at each location.
+#' The altitude value of the geolocator pressure timeseries is computed with the barometric formula
+#' accounting for the temporal variation of pressure (surface-pressure) and temperature
+#' (2m-temperature) based on ERA5 data. The min and max ground elevation of each pixel is computed
+#' from SRTM-90. This map is only returned if `keep_mask` is `TRUE`.
+#'
+#' For more details, read the [GeoPressure API documentation
+#' ](https://raphaelnussbaumer.com/GeoPressureAPI/).
+#'
+#' @section Elevation levels:
+#'
+#' It is possible to indicate different elevation levels when the bird was spending time at
+#' locations with different elevations within a general area (~10km), and thus within the same
+#' stationary period. This can be done by using `tag$label="elev_x"`for all measurements of the same
+#' elevation level *x*. See more information on the labeling of elevation levels in [the
+#' corresponding section in the GeoPressureManual](
+#' https://raphaelnussbaumer.com/GeoPressureManual/labelling-tracks.html#elevation-period).
+#'
+#' Behind the scen, each of these elevations levels produce a new requests on the GeoPressureAPI.
+#' The mismatch maps of all elevation levels belonging to the same stationay periods are combined
+#' (as a weighted average) to results in a single mismatch maps per stationary period.
+#'
+#' @section Convert mismatch map into likelihood map:
+#'
+#' We convert the map of the mean square error \eqn{MSE} and altitude mask \eqn{z_{mask}} computed
+#' by [`geopressure_map_mismatch()`] into a likelihood map with,
+#'
+#' \deqn{L = \exp \left(-w \frac{MSE}{\sigma} \right) \left[z_{mask}>T \right],}
+#'
+#' where \eqn{\sigma} is the standard deviation (`sd`) of pressure and \eqn{T} is the mask threshold
+#' (`thr_mask`).
+#'
+#' Because the auto-correlation of the timeseries is not accounted for in this equation, we use a
+#' log-linear pooling weight of \eqn{w=\log(n)/n} by default, where \eqn{n} is the number of samples
+#' in the timeserie (i.e., data points used to compute the MSE). See [GeoPressureManual |
+#' Probability aggregation
+#' ](https://raphaelnussbaumer.com/GeoPressureManual/probability-aggregation.html) for details.
+#'
+#' **Important Note**: Since GeoPressure v3.1.0, the threshold of the mask is happening directly
+#' on the GEE server (i.e., during `geopressure_map_mismatch()`). This allows to compute the MSE
+#' only for pixels which are within the threshold, thus reducing the computational cost
+#' significantly.
 #'
 #' @param tag a GeoPressureR `tag` object
-#' @param max_sample The computation of the maps is only performed on `max_sample` datapoints of
+#' @param max_sample the computation of the maps is only performed on `max_sample` datapoints of
 #' pressure to reduce computational time. The samples are randomly (uniformly) selected on the
 #' timeseries.
-#' @param margin The margin is used in the mask map to accept measurement errors, small-scale
+#' @param margin the margin is used in the mask map to accept measurement errors, small-scale
 #' topography, and vertical movements of the bird (unit in meters, 1hPa~10m).
-#' @param sd Standard deviation of the pressure error . numeric of lenght 1 or number of stationary
+#' @param sd standard deviation of the pressure error . numeric of lenght 1 or number of stationary
 #' periods.
-#' @param thr_mask Threshold of the percentage of data points outside the elevation range to be
+#' @param thr_mask threshold of the percentage of data points outside the elevation range to be
 #' considered not possible.
-#' @param log_linear_pooling_weight Weighting function of the log-linear pooling, taking the number
+#' @param log_linear_pooling_weight weighting function of the log-linear pooling, taking the number
 #' of samples of the stationary periods used and returning the weight of the aggregation. See
 #' [GeoPressureManual | Probability aggregation
 #' ](https://raphaelnussbaumer.com/GeoPressureManual/probability-aggregation.html) for more details.
-#' @param keep_mse_mask logical defining if the matrix of the MSE and mask are returned in the list.
+#' @param keep_mask logical defining if the mask map is returned in `tag`.
+#' @param keep_mse logical defining if the MSE map is returned in `tag`.
 #' @param timeout Duration before the code is interrupted both for the request on
-#' GeoPressureAPI and GEE (in seconds, see [`httr::timeout()`]).
+#' GeoPressureAPI and on GEE (in seconds, see [`httr::timeout()`]).
 #' @param workers Number of parallel requests on GEE. Integer between 1 and 99. `"auto"` adjust the
 #' number of workers to the number of `stap_elev` to query.
 #' @param compute_known logical defining if the map(s) for known stationary period should be
 #' estimated based on twilight or hard defined by the known location `stap$known_l**`
 #' @param quiet logical to hide messages about the progress
+#' @param debug logical to display additional information to debug a request
 #'
-#' @return Same list as parameter `tag` but with a `likelihood` field containing a list of maps
-#' for each stationary period stored as matrix. Note that stationary periods not marked as model in
-#' `tag$stap$include` will be included with a `NULL`value in `likelihood`.
-#' `tag$param` is a list of all the parameters used to compute the likelihood map.
-#' If `keep_mse_mask` is
-#' true, the `mse` and `mask` maps computed with `geopressure_map_mismatch` will also be kept, as
-#' well as `tag$stap$nb_sample`, indicating the number of datapoints used to compute the MSE.
+#' @return Returns the same GeoPressureR `tag` object including the GeoPressureR `map` object
+#' `tag$map_pressure` containing the likelihood map of each statationay period. See `map_create()`
+#' for details.
+#' If `keep_mask` and `keep_mse` are each true, `tag` also includes the `tag$map_pressure_mse` and
+#' `tag$map_pressure_mask` maps, as well as `tag$stap$nb_sample`, indicating the number of
+#' datapoints used to compute the MSE.
 #'
 #' @examples
-#' setwd(system.file("extdata/", package = "GeoPressureR"))
+#' setwd(system.file("extdata", package = "GeoPressureR"))
 #' tag <- tag_create("18LX", quiet = TRUE) |>
 #'   tag_label(quiet = TRUE) |>
 #'   tag_set_map(
 #'     extent = c(-16, 23, 0, 50),
 #'     scale = 4
 #'   )
 #'
-#' tag <- geopressure_map(tag,
+#' tag <- geopressure_map_mismatch(tag,
 #'   max_sample = 50,
-#'   sd = 0.7,
-#'   keep_mse_mask = TRUE,
-#'   quiet = TRUE
+#'   margin = 20,
+#'   thr_mask = 0.95,
+#'   keep_mask = TRUE,
+#'   quiet = TRUE,
+#'   workers = 1
 #' )
 #'
 #' plot(tag, type = "map_pressure_mse", plot_leaflet = FALSE)
 #'
 #' plot(tag, type = "map_pressure_mask", plot_leaflet = FALSE)
 #'
+#' tag <- geopressure_map_likelihood(tag,
+#'   sd = 1,
+#'   log_linear_pooling_weight = function(n) log(n) / n
+#' )
+#'
 #' plot(tag, type = "map_pressure")
 #'
 #' @references{ Nussbaumer, Raphaël, Mathieu Gravey, Martins Briedis, and Felix Liechti. 2023.
@@ -82,13 +161,16 @@ geopressure_map <- function(tag,
                             sd = 1,
                             thr_mask = 0.9,
                             log_linear_pooling_weight = \(n) log(n) / n,
-                            keep_mse_mask = FALSE,
+                            keep_mask = FALSE,
+                            keep_mse = FALSE,
                             compute_known = FALSE,
                             quiet = FALSE) {
   # Compute mean square error maps
   tag <- geopressure_map_mismatch(tag,
     max_sample = max_sample,
     margin = margin,
+    thr_mask = thr_mask,
+    keep_mask = keep_mask,
     timeout = timeout,
     workers = workers,
     compute_known = compute_known,
@@ -98,15 +180,9 @@ geopressure_map <- function(tag,
   # Compute likelihood maps from the MSE maps
   tag <- geopressure_map_likelihood(tag,
     sd = sd,
-    thr_mask = thr_mask,
-    log_linear_pooling_weight = log_linear_pooling_weight
+    log_linear_pooling_weight = log_linear_pooling_weight,
+    keep_mse = keep_mse
   )
 
-  # remove intermediate maps computed by geopressure_map_mismatch()
-  if (!keep_mse_mask) {
-    tag[names(tag) %in% c("map_pressure_mse", "map_pressure_mask")] <- NULL
-    tag$stap <- tag$stap[names(tag$stap) != "nb_sample"]
-  }
-
   return(tag)
 }