Silently drop wl where sensdata is 0 if required

NEWS.md

@@ -14,6 +14,10 @@
 
 ## MINOR FEATURES AND BUG FIXES
 
+* `vismodel()` now works with spectral data that doesn't span the whole 
+300-700nm range as long as the spectral sensitivity is zero for the missing 
+range. This is useful if you have data covering the 390-700nm range and work 
+the CIE model for example.
 * `getimg()` now imports image files with uppercase extensions (e.g., JPG or 
 PNG), such as those produced by some camera brands or processing software.
 * Maximum quantum catches computation (`data.maxqcatches` attribute) now works

R/sensdata.R

@@ -2,6 +2,7 @@
 #'
 #' Retrieve (as an rspec object) or plot pavo's in-built spectral sensitivity data.
 #'
+#' @inheritParams sensmodel
 #' @param visual visual systems. Options are:
 #' - `"none"`: no visual sensitivity data.
 #' - `"all"`: all visual sensitivity data.
@@ -77,14 +78,17 @@ sensdata <- function(visual = c(
                      trans = c("none", "all", "bluetit", "blackbird"),
                      bkg = c("none", "all", "green"),
                      plot = FALSE,
+                     range = c(300, 700),
                      ...) {
   visual2 <- match.arg(visual, several.ok = TRUE)
   achro2 <- match.arg(achromatic, several.ok = TRUE)
   illum2 <- match.arg(illum, several.ok = TRUE)
   bkg2 <- match.arg(bkg, several.ok = TRUE)
   trans2 <- match.arg(trans, several.ok = TRUE)
 
-  dat <- data.frame("wl" = 300:700)
+  wl <- seq(range[1], range[2])
+
+  dat <- data.frame("wl" = wl)
 
   # Visual system
   if (!isTRUE("none" %in% visual2)) {
@@ -95,7 +99,7 @@ sensdata <- function(visual = c(
         "ctenophorus"
       )
     }
-    S <- vissyst[, grepl(paste(visual2, collapse = "|"), names(vissyst)), drop = FALSE]
+    S <- vissyst[vissyst$wl %in% wl, grepl(paste(visual2, collapse = "|"), names(vissyst)), drop = FALSE]
     dat <- cbind(dat, S)
   }
 
@@ -104,7 +108,7 @@ sensdata <- function(visual = c(
     if (isTRUE("all" %in% achro2)) {
       achro2 <- c("bt.dc", "ch.dc", "st.dc", "md.r1", "ra.dc", "cf.r")
     }
-    achro <- vissyst[, grepl(paste(achro2, collapse = "|"), names(vissyst)), drop = FALSE]
+    achro <- vissyst[vissyst$wl %in% wl, grepl(paste(achro2, collapse = "|"), names(vissyst)), drop = FALSE]
     dat <- cbind(dat, achro)
   }
 
@@ -113,7 +117,7 @@ sensdata <- function(visual = c(
     if (isTRUE("all" %in% illum2)) {
       illum2 <- c("bluesky", "D65", "forestshade")
     }
-    illum <- bgandilum[, grepl(paste(illum2, collapse = "|"), names(bgandilum)), drop = FALSE]
+    illum <- bgandilum[bgandilum$wl %in% wl, grepl(paste(illum2, collapse = "|"), names(bgandilum)), drop = FALSE]
     dat <- cbind(dat, illum)
   }
 
@@ -122,7 +126,7 @@ sensdata <- function(visual = c(
     if (isTRUE("all" %in% bkg2)) {
       bkg2 <- "green"
     }
-    bkg <- bgandilum[, grepl(paste(bkg2, collapse = "|"), names(bgandilum)), drop = FALSE]
+    bkg <- bgandilum[bgandilum$wl %in% wl, grepl(paste(bkg2, collapse = "|"), names(bgandilum)), drop = FALSE]
     dat <- cbind(dat, bkg)
   }
 
@@ -131,7 +135,7 @@ sensdata <- function(visual = c(
     if (isTRUE("all" %in% trans2)) {
       trans2 <- c("bluetit", "blackbird")
     }
-    trans <- transmissiondata[, grepl(paste(trans2, collapse = "|"), names(transmissiondata)), drop = FALSE]
+    trans <- transmissiondata[transmissiondata$wl %in% wl, grepl(paste(trans2, collapse = "|"), names(transmissiondata)), drop = FALSE]
     dat <- cbind(dat, trans)
   }
 

R/vismodel.R

@@ -267,10 +267,6 @@ vismodel <- function(rspecdata,
     )
   }
 
-  if (!isTRUE(all.equal(wl, sens_wl, check.attributes = FALSE))) {
-    stop("wavelength range in spectra and visual system data do not match")
-  }
-
   # DEFINING ILLUMINANT & BACKGROUND
 
   bgil <- bgandilum
@@ -279,14 +275,14 @@ vismodel <- function(rspecdata,
     illum <- bgil[, grep(illum2, names(bgil))]
   }
   if (illum2 == "ideal") {
-    illum <- rep(1, dim(rspecdata)[1])
+    illum <- rep_len(1, 401)
   }
 
   if (bg2 != "user-defined") {
     bkg <- bgil[, grep(bg2, names(bgil))]
   }
   if (bg2 == "ideal") {
-    bkg <- rep(1, dim(rspecdata)[1])
+    bkg <- rep_len(1, 401)
   }
 
   # Defining ocular  <- mission
@@ -296,7 +292,7 @@ vismodel <- function(rspecdata,
     trans <- trdat[, grep(tr2, names(trdat))]
   }
   if (tr2 == "ideal") {
-    trans <- rep(1, dim(rspecdata)[1])
+    trans <- rep_len(1, 401)
   }
 
   if (tr2 != "ideal" & visual2 == "user-defined") {
@@ -338,6 +334,18 @@ vismodel <- function(rspecdata,
   illum <- prepare_userdefined(illum)
   achromatic <- prepare_userdefined(achromatic)
 
+  if (!isTRUE(all.equal(wl, sens_wl, check.attributes = FALSE))) {
+    if (length(sens_wl) > length(wl) && all(S[!sens_wl %in% wl, ] == 0)) {
+      S <- S[sens_wl %in% wl, ]
+      trans <- trans[sens_wl %in% wl]
+      bkg <- bkg[sens_wl %in% wl]
+      illum <- illum[sens_wl %in% wl]
+      achromatic <- achromatic[sens_wl %in% wl]
+    } else {
+      stop("wavelength range in spectra and visual system data do not match")
+    }
+  }
+
   # Transform from percentages to proportions (Vorobyev 2003)
   if (max(y) > 1) {
     y <- y / 100

data-raw/sysdata.R

@@ -3,31 +3,41 @@ library(pavo)
 # FIXME: at some point, this line should be removed because all data will be
 # generated from this file and we won't have to rely on existing data anymore.
 load("R/sysdata.rda")
+bgandilum <- bgandilum[, c("wl", "bluesky", "forestshade", "green")]
+vissyst <- vissyst[, !grepl("cie", colnames(vissyst))]
 
 # bgandilum
 
-d65 <- readxl::read_xls("data-raw/ciedata.xls", skip = 5, col_names = c("wl", "d65"), sheet = "D65")
-d65 <- as.rspec(d65, lim = c(300, 700))
+d65 <- readxl::read_xls("data-raw/ciedata.xls", skip = 5, col_names = c("wl", "D65"), sheet = "D65")
+d65 <- as.rspec(d65)
 d65 <- procspec(d65, opt = "maximum")
 
-bgandilum$D65 <- d65$d65
-bgandilum$ideal <- 1
+bgandilum <- merge(bgandilum, d65, all = TRUE)
+bgandilum[is.na(bgandilum)] <- 0
 
 # transmissiondata
 
 # vissyst
 
-cie2 <- readxl::read_xls("data-raw/ciedata.xls", range = "1931 col observer!A6:D86", col_names = c("wl", "x", "y", "z"))
-cie2 <- as.rspec(cie2, lim = c(300, 700), exceed.range = FALSE)
-cie2[is.na(cie2)] <- 0
+cie2 <- readxl::read_xls(
+  "data-raw/ciedata.xls",
+  range = "1931 col observer!A6:D86",
+  col_names = c("wl", paste0("cie2_", c("X", "Y", "Z")))
+)
+cie2 <- as.rspec(cie2)
+
+vissyst <- merge(vissyst, cie2, all = TRUE)
 
-vissyst[, paste0("cie2_", c("X", "Y", "Z"))] <- cie2[, c("x", "y", "z")]
+cie10 <- readxl::read_xls(
+  "data-raw/ciedata.xls",
+  range = "1964 col observer!A6:D86",
+  col_names = c("wl", paste0("cie10_", c("X", "Y", "Z")))
+)
+cie10 <- as.rspec(cie10)
 
-cie10 <- readxl::read_xls("data-raw/ciedata.xls", range = "1964 col observer!A6:D86", col_names = c("wl", "x", "y", "z"))
-cie10 <- as.rspec(cie10, lim = c(300, 700), exceed.range = FALSE)
-cie10[is.na(cie10)] <- 0
+vissyst <- merge(vissyst, cie10, all = TRUE)
 
-vissyst[, paste0("cie10_", c("X", "Y", "Z")]) <- cie10[, c("x", "y", "z")]
+vissyst[is.na(vissyst)] <- 0
 
 usethis::use_data(
   bgandilum,

tests/testthat/test-vismodel.R

@@ -35,6 +35,11 @@ test_that("Warnings", {
 
   expect_silent(vismodel(flowers_NIR, sensmodel(c(350, 450, 550, 650), range = c(300, 1200))))
 
+  expect_identical(
+    vismodel(flowers, visual = "cie10", illum = "D65", vonkries = TRUE, relative = FALSE),
+    vismodel(flowers[91:401,], visual = "cie10", illum = "D65", vonkries = TRUE, relative = FALSE)
+  )
+
 })
 
 test_that("Sensmodel", {
@@ -63,8 +68,13 @@ test_that("sensdata()", {
   # No negative values, no NA
   expect_false(any(vis_all < 0))
   expect_false(anyNA(vis_all))
+  expect_identical(dim(vis_all), c(401L, 57L))
+  expect_s3_class(vis_all, "rspec")
 
-  colspace(vismodel(sensdata(illum = "D65"), visual = "cie10"))
+  vis_part <- sensdata(visual = "cie2", range = c(400, 700))
+
+  expect_identical(dim(vis_part), c(301L, 4L))
 
+  colspace(vismodel(sensdata(illum = "D65"), visual = "cie10"))
 
 })