Correct R syntax for changing draw() plot - gratia package

Question

I need help translating "modelling the long-term trend that varies spatially in my data" and for changing the plot below into the correct R syntax.

I'd like to have Longitude on the x-axis, latitude on the y-axis, and a panel for each calendar year ("CYR"). Changing the order of the variables in my smooth function gives me this error though:

     ...ti(Latitude, CYR, Longitude, d = c(2,1), bs = c('ds','tp'), k = c(25, 14)), ..
     ...ti(CYR, Latitude, Longitude, d = c(2,1), bs = c('ds','tp'), k = c(25, 14)), ..

Error in check.term(termi, rec) : 
  bam can not discretize with this nesting structure

It's also somehow connected to my other covariates as well (removing them stops the error, but I can't get the plot to "draw" correctly).

library(mgcv)
library(gratia)

m <- bam(occur ~ s(temp) + 
           sal + 
           s(DO) + 
           sed_depth + 
           water_depth + 

           s(fCYR, bs = "re") + # Long-term trend
           
           # Spatial variation
           s(Longitude, Latitude, k = 100, bs = 'tp') + # Try bs = tp (default)
           
           s(fSite, bs = "re") + # Repeated measures design
         
         # long-term trend varies spatially (code I need to change, I think)
         ti(Latitude, CYR, Longitude, d = c(2,1), bs = c('ds','tp'), k = c(25, 14)),
         
         data = toad, 
         method = 'fREML', 
         # knots = knots,
         nthreads = 4, 
         discrete = TRUE,
         family = binomial(link = "logit"), 
         # select = TRUE,
         gamma = 1.5)

draw(m, select = 6)

Answer 1

As of version 0.8.1.25 (the development version of the package at the time of writing), {gratia} should now be able to handle this automatically.

library("gratia")
library("mgcv")
library("ggplot2")
library("dplyr")

df <- data_sim("eg1", n = 1000,  dist = "normal", scale = 2, seed = 1)
m <- gam(y ~ te(x0, x1, x2, k = c(25, 10), d = c(1,2), bs = c("cr", "ds")),
         data = df, method = "REML")
draw(m, n = 25)

produces

which, as can be seen, has plotted the 2D Duchon spline in the panels and x0 on the facets even though this is not how the terms were ordered when the te() term was defined in the formula.

Contray to my comment on the OP's question, this (draw.gam() behaviour) was only added in 0.8.1.25; .24 added support for smooth_estimates() and its draw() method only.

I posted this as a separate answer because this new behaviour only works for 2D marginal smooths; if we had te(x0, x1, x2, x3, bs = c("cr", "ds"), d = c(1, 3)), the code currently doesn't identify that we could do better by generating data for x1 and x2 to go on the x and y axes of the panels, and then facet by x3 and x1. In that case you will still need the ideas in the original answer I posted to generate the desired output manually.