How to structure random slope with interactions using glmmTMB in R

I am trying to:

Build a model to investigate the effect of landcover (a continuous variable) on presence of a species.
Test if and how this effect changes with season (a categorical variable, spring and fall).

However, no animal in my data set has observations in both seasons.

I am using R package glmmTMB to estimate a binomial GLMM with random intercept and random slope for animal ID, but I don't know how to represent the effects of landcover (the variable of interest) and season in the model formula.

I have tried these two variations:

1)

glmmTMB(presence ~ landcover:season + (1 | ID) +  (0 + landcover:season | ID),  
        family = binomial(), data = dat.rsf, doFit = FALSE, weights = weight)

This one is more intuitive to me, but does not converge. Is that just an indication that I cannot estimate this model with my data? Or is this model not the best way to get at the relationships that I am trying to test?

2)

glmmTMB(presence ~ landcover:season + (1 | ID) +  (0 + landcover | ID), 
        family = binomial(), data = dat.rsf, doFit = FALSE, weights = weight)

This converges, but is not as intuitive. Does it make sense to not include the interaction in the random slope when I include it as a fixed effect? Or does that violate the principle of a random slope?

Solution

tl;dr I would suggest

glmmTMB(presence ~ landcover*season + (1 + landcover | ID), 
        family = binomial(), data = dat.rsf, 
         weights = weight,
         map = list(theta=factor(c(NA, 1, 2))),
         start = list(theta = c(log(1e3), rep(0, 2))
    )

(noting that landcover is a continuous predictor).

If you're using this for habitat selection analysis and need to fix the values of some parameters, see e.g. here and search for "one-step way" to see how to use the map and start parameters to achieve your goal ... what I've done here is to fix the intercept variation among individuals to a large value (SD=1000). The other two parameters in theta are the (among-individual variation in land cover effects, and the correlations between intercept and slope.
You can separate (1|ID) and (0+landcover|ID); this will estimate intercept and landcover-effect variation, but assume tha they are independent for each individual. That's not a ridiculous assumption, but the model listed here is more general. diag(landcover|ID) is equivalent to (1|ID) + (0+landcover|ID)
including season as a term as in season|ID doesn't work because it's not possible to estimate among-individual variation in season effects when no individuals are observed in multiple seasons
As for landcover:season: in general, it doesn't make sense to include an interaction term without including the main effects. landcover*season is equivalent to landcover + season + landcover:season, and would be preferred. (There are sometimes reasons to violate marginality by including interactions without their corresponding main effects, but if you're going to do this you should know exactly why you're doing it ...)