Add a vignette, and fix the row reordering it exposed

.Rbuildignore

@@ -3,3 +3,5 @@
 ^LICENSE\.md$
 ^\.github$
 ^README\.Rmd$
+^vignettes/.*\.html$
+^vignettes/.*_files$

.gitignore

@@ -8,3 +8,6 @@
 
 # Reference papers kept locally, not distributed with the package
 *.pdf
+
+# knitr scratch output from building the vignette outside the package machinery
+figure/

DESCRIPTION

@@ -38,9 +38,12 @@ Imports:
 Suggests:
     marmap,
     ncdf4,
+    knitr,
+    rmarkdown,
     testthat (>= 3.0.0),
     withr
 Config/testthat/edition: 3
+VignetteBuilder: knitr
 URL: https://github.com/chross22/datamatch
 BugReports: https://github.com/chross22/datamatch/issues
 Config/roxygen2/version: 8.1.0

NEWS.md

@@ -54,6 +54,12 @@
   daily dataset was treated as monthly and only day 1 of each month was
   downloaded. The frequency token is now matched as such.
 
+* **matchData() returned rows grouped by period, not in the order they went
+  in.** Rows are processed a period at a time, so a table whose periods were
+  interleaved came back reordered. Anyone aligning the result against the input
+  by position - `cbind()`, or assigning a column straight across - would have
+  got silently mismatched rows. The order is now restored before returning.
+
 * **Sparse daily data could be matched as though it were monthly.**
   `detect_temporal_resolution()` infers daily data from more than one day within
   a month, so a set of survey dates — one per month — was indistinguishable from

R/matchData.R

@@ -122,6 +122,16 @@ matchData <- function(dat, source,
     source_vars[match(collisions, source_vars)] <- renamed
   }
 
+  # Rows are processed a period at a time, so they come back grouped by period
+  # rather than in the order they arrived. That is a quiet hazard: anyone
+  # aligning the result against the input by position - cbind(), or assigning a
+  # column straight across - would get silently mismatched rows. The original
+  # position is carried through and used to restore the order at the end.
+  #
+  # The name is deliberately awkward so it cannot collide with a real column.
+  order_key <- ".datamatch_row_order"
+  dat[[order_key]] <- seq_len(nrow(dat))
+
   periods <- unique(sf::st_drop_geometry(dat)[match_keys])
   matched <- vector("list", nrow(periods))
   unmatched_periods <- character()
@@ -164,6 +174,11 @@ matchData <- function(dat, source,
   matched <- lapply(matched, function(x) x[column_order])
   matched_data <- do.call(rbind, matched)
 
+  # Back into the order the rows arrived in, and drop the bookkeeping column.
+  matched_data <- matched_data[order(matched_data[[order_key]]), ]
+  matched_data[[order_key]] <- NULL
+  rownames(matched_data) <- NULL
+
   matched_data$LON <- sf::st_coordinates(matched_data)[, 1]
   matched_data$LAT <- sf::st_coordinates(matched_data)[, 2]
 

tests/testthat/test-matchData.R

@@ -306,3 +306,36 @@ test_that("neither side has to be observations or a covariate grid", {
   expect_true(all(c("SST", "CHL") %in% names(result)))
   expect_false(anyNA(result$CHL))
 })
+
+test_that("rows come back in the order they went in", {
+  # Rows are processed a period at a time, so without restoring the order they
+  # come back grouped by period. Anyone aligning the result against the input by
+  # position - cbind(), or assigning a column straight across - would get
+  # silently mismatched rows, which is the worst kind of wrong.
+  env <- sf::st_as_sf(
+    do.call(rbind, lapply(1:3, function(m) {
+      g <- expand.grid(x = c(-70, -69), y = c(42, 43))
+      g$SST <- m
+      g$YEAR <- 2015L
+      g$MONTH <- as.integer(m)
+      g$DAY <- 1L
+      g
+    })),
+    coords = c("x", "y"), crs = 4326)
+
+  # Deliberately out of period order, and with a period repeated.
+  observations <- sf::st_as_sf(
+    data.frame(lon = rep(-69.5, 4), lat = 42.5, YEAR = 2015L,
+               MONTH = c(3L, 1L, 3L, 2L), id = 1:4),
+    coords = c("lon", "lat"), crs = 4326)
+
+  result <- matchData(observations, env)
+
+  expect_equal(result$id, observations$id)
+  # And the covariate follows its own row, not merely the row count.
+  expect_equal(result$SST, observations$MONTH)
+
+  # The bookkeeping column used to restore the order must not leak out.
+  expect_false(any(grepl("datamatch_row", names(result))))
+  expect_equal(rownames(result), as.character(seq_len(nrow(result))))
+})

vignettes/datamatch.Rmd

@@ -0,0 +1,218 @@
+---
+title: "Getting started with datamatch"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Getting started with datamatch}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.2)
+library(datamatch)
+set.seed(1)
+```
+
+This walks a set of observations from raw positions to a table a model can be
+fitted to, and shows what each step actually returns.
+
+Every number and figure below is computed when the vignette is built. Nothing is
+downloaded: the covariates are synthetic, generated to look like a Gulf of Maine
+grid with a seasonal cycle and satellite cloud gaps. The calls are the ones you
+would make against real data, so a fetch is the only thing to substitute.
+
+```{r fetch-for-real, eval = FALSE}
+bb <- list(xmin = -70, xmax = -66, ymin = 41, ymax = 44)
+
+env <- accessEnvDat(vars = c("SST", "MLD"), years = 2015, months = 1:12,
+                    bounding_box = bb)
+```
+
+```{r synthetic, echo = FALSE}
+grid <- expand.grid(x = seq(-70, -66, by = 0.25), y = seq(41, 44, by = 0.25))
+
+make_month <- function(m) {
+  d <- grid
+  # Warm in late summer, cool in late winter; colder to the north.
+  d$SST <- 9 + 5.5 * sin(2 * pi * (m - 4.5) / 12) - 1.4 * (d$y - 41) +
+    stats::rnorm(nrow(d), 0, 0.25)
+  # Deep mixed layer in winter, shallow in summer.
+  d$MLD <- 45 - 32 * sin(2 * pi * (m - 4.5) / 12) + stats::rnorm(nrow(d), 0, 3)
+  d$YEAR <- 2015L
+  d$MONTH <- as.integer(m)
+  d$DAY <- 1L
+  d
+}
+
+env <- sf::st_as_sf(do.call(rbind, lapply(1:12, make_month)),
+                    coords = c("x", "y"), crs = 4326)
+attr(env, "datamatch_step") <- "month"
+```
+
+`accessEnvDat()` returns one row per grid cell per time step, as an `sf` point
+object:
+
+```{r shape}
+env
+```
+
+## Look before you model
+
+Three plots answer the questions worth asking of a download, and each returns the
+data it drew rather than only drawing it.
+
+`plot_env()` maps one variable for one time step. This is the fastest way to
+catch a bounding box that landed somewhere unintended, or a variable that is
+entirely `NA`:
+
+```{r plot-env}
+plot_env(env, "SST", time = c(MONTH = 8))
+```
+
+`plot_series()` reduces each step to one number over the study area, so the
+seasonal cycle and its spatial spread are both visible:
+
+```{r plot-series}
+series <- plot_series(env, "SST")
+head(series, 4)
+```
+
+## Gaps are not spread evenly
+
+Satellite variables are missing wherever cloud blocked the view, and those gaps
+cluster in particular seasons. Adding a chlorophyll field with realistic winter
+gaps:
+
+```{r gappy}
+env$CHL <- 1.2 + 2.4 * exp(-((env$MONTH - 5)^2) / 5) + stats::rnorm(nrow(env), 0, 0.1)
+
+winter <- env$MONTH %in% c(11, 12, 1, 2)
+env$CHL[winter][sample(sum(winter), round(sum(winter) * 0.8))] <- NA
+env$CHL[!winter][sample(sum(!winter), round(sum(!winter) * 0.1))] <- NA
+```
+
+`plot_coverage()` is the one to run before trusting a monthly mean:
+
+```{r plot-coverage}
+coverage <- plot_coverage(env, c("SST", "CHL"))
+```
+
+```{r coverage-numbers}
+round(coverage$coverage[coverage$variable == "CHL"], 2)
+```
+
+A quarter of the grid in winter and near-complete in summer is exactly the shape
+that decides whether a winter value means anything.
+
+## Matching
+
+`matchData()` joins each observation to the nearest cell within the same time
+period. Some observations here are deliberately placed in a month the covariates
+do not cover, and some outside the grid:
+
+```{r observations}
+observations <- sf::st_as_sf(
+  data.frame(
+    lon   = c(-69.5, -68.2, -67.1, -66.5, -69.0, -64.0),
+    lat   = c(41.5,  42.3,  43.1,  41.9,  43.6,  42.0),
+    YEAR  = 2015L,
+    MONTH = c(3L, 6L, 8L, 6L, 11L, 6L),
+    count = c(12, 45, 7, 88, 3, 21)
+  ),
+  coords = c("lon", "lat"), crs = 4326)
+
+matched <- matchData(observations, env)
+sf::st_drop_geometry(matched)[c("MONTH", "count", "SST", "MLD", "CHL")]
+```
+
+Every observation comes back, in the same order, with the covariates attached.
+The last one sits well outside the grid and still matched — `matchData()` uses
+the *nearest* cell, and nearest has no maximum distance. That is worth knowing:
+an observation far outside the study area is joined to the closest edge cell
+rather than dropped, so check your bounding box covers your stations.
+
+The November observation has `NA` chlorophyll, because that cell was under cloud.
+That is a real gap rather than a failure, and the source of it is worth checking
+before dropping the row:
+
+```{r missing}
+colSums(is.na(sf::st_drop_geometry(matched)))
+```
+
+## Choosing a resolution
+
+Real products do not share a grid. Physics is 0.083°, biogeochemistry 0.25°, and
+satellite ocean colour 4 km, so combining them means deciding which grid to keep.
+
+Aggregating is the safe direction, because every value in the result summarises
+values that were really measured:
+
+```{r upscale}
+coarse <- upscale_grid(env, to = 0.5, vars = "SST", min_coverage = 0)
+
+c(cells_before = nrow(env) / 12, cells_after = nrow(coarse) / 12)
+```
+
+Interpolating is the direction to be careful in. It adds cells, not information:
+
+```{r downscale}
+fine <- downscale_grid(coarse, to = 0.25, vars = "SST")
+
+# `nearest`, the default, invents no values: every one was already in the source.
+all(stats::na.omit(fine$SST) %in% coarse$SST)
+```
+
+`bilinear` would return a smooth field that looks like a finely-resolved
+measurement and is not, which is why the blunt method is the default.
+
+### Partial coverage comes back NA
+
+Aggregating chlorophyll over the winter, when most cells are empty, is the case
+`min_coverage` exists for:
+
+```{r min-coverage}
+january <- env[env$MONTH == 1, ]
+
+strict <- upscale_grid(january, to = 0.5, vars = "CHL", min_coverage = 0.5)
+loose  <- upscale_grid(january, to = 0.5, vars = "CHL", min_coverage = 0)
+
+c(reported_at_half_coverage = sum(!is.na(strict$CHL)),
+  reported_with_no_guard    = sum(!is.na(loose$CHL)))
+```
+
+Both numbers are available; only one of them says how much went into it.
+`keep_counts = TRUE` returns the fraction behind each value.
+
+## Covariates that are not gridded
+
+Two other kinds attach differently. Seafloor terrain is static, so the same
+value goes to every time step at a location:
+
+```{r bathymetry, eval = FALSE}
+bathy <- fetch_bathymetry(bounding_box = bb)
+matched <- attach_bathymetry(matched, bathy, c("DEPTH", "SLOPE", "TPI"))
+```
+
+Climate indices have no spatial dimension at all — one value per month describes
+the whole basin, so every observation in a month receives the same number:
+
+```{r indices, eval = FALSE}
+matched <- attach_climate_index(matched, c("NAO", "AMOC"))
+```
+
+That makes them a different kind of covariate. They carry information about
+*when*, and none about *where* within a region conditions are better. A model
+given only indices cannot produce a map.
+
+```{r index-table}
+as.data.frame(index_dictionary())[c("name", "units", "source")]
+```
+
+## Where to go next
+
+- `vignette("datamatch")` is this document; the
+  [README](https://github.com/chross22/datamatch) covers the same ground in more
+  depth, including forecasts, daily data, and gap filling.
+- `variable_dictionary()` and `index_dictionary()` list what can be fetched.
+- The README's References section lists the DOI for every data source, since the
+  obligation to cite travels with the data.