How to interpret the plotted difference of in-factor smooths I used the gamm4 R package to analyze the location effect on heart rate over time (Inc_age), including some linear co-variates and a random effect. My model looks like this:
model <- gamm4(Heart rate ~ s(Inc_age, by = Location) + Location + Time_dec +
               Mass * Egg_temp, data = hr1, random = ~ (1 |Egg_ID),
               REML = TRUE)

I checked for colinearity of linear predictors, concurvity, autocorrelation, adequate k, and also the residuals look fine. Now I am struggling a bit with the interpretation of the results.
1) The summary table tells me that both smoothers per location are significant. The parametric coefficients table tells me that location has a significant effect on heart rate. The linear fixed effects and the interaction term are not significant. R-squared is 0.73. So far, so good. Now, I want to know in which exact time window the locations differ.

2) When plotting the two smoothers (one per location) the confidence levels are overlapping in most parts of the curves, only one part seems to be just not overlapping. Is this non-overlapping time window the only one in which the locations differ? This is the code I used:
plot_smooth(model$gam, view="Inc_age", cond=list(Location="North"), 
            rm.ranef=TRUE, rug=FALSE, col="red", ylim=c(0,400))
plot_smooth(model$gam, view="Inc_age", cond=list(Location="South"), 
            rm.ranef=TRUE, rug=FALSE, col="cyan", add=TRUE)


3) I also try to plot the difference between the smooths with the following code:
model$gam1 <- getViz(model$gam)

plotDiff(s1 = sm(model$gam1, 1), s2 = sm(model$gam1, 2)) + l_ciPoly() + 
         l_fitLine() + geom_hline(yintercept = 0, linetype = 2)

The y-axis is s(North)-s(South) = SD = the difference between the fitted smooth. X-axis is time. It shows a single curve with a confidence interval. Again, there is one time window in which the confidence bands do not cross the 0 line. This is not the same time window that I described in (2).

My question is whether this is the correct way to find the exact differences between locations. And if so, how do I interpret the plots? Is there also a way to calculate the exact time window instead of visualizing it?
Follow-up
I am following the steps you presented on link. 
I adjusted it slightly because I only have two groups (not more) that I want to compare. When I try to calculate the difference between the smoothers ( dif <- X %*% coef(model) ), I get an error saying "non-conformable arguments". I assume this is because I run this on the model that also includes my linear predictors. Because it does work for the same model, excluding those. If understand correctly, that is also what you describe in your comment above. 
Then, when I calculate dif, se, upr, and lwr of the reduced model (excluding linear covariates) and I plot this against time (inc_age), I get indeed the same plot as the second one in my original question. So, I guess this confirms that plotDiff computes the same as in your example. 
Then, !(upr_ci > 0 & lwr_ci < 0) tells me that the difference lies in inc_age 6 and 7. However, when I plot the raw data the lines are overlapping at that stage and the difference seems to be from inc_age 9 onwards. I am still struggling to understand this. Or did I make a mistake somewhere?

 A: The first plot showing the two smooths also includes the model constant term (the intercept), whereas the second plot is showing the difference between the two smooths themselves, excluding the constant term. As such you have to be careful to distinguish what you mean by the difference between the two smooths; the first figure conflates differences in group means and differences of smooth effects of the covariate between the groups.
Typically, however, one is interested in the difference in the estimated smooth effect of the covariate. As such the second plot is the one I would typically use to look at differences between smooths.
That said, the first plot is useful because it shows more clearly the differences between groups (i.e. groups means plus by group smooth effects).
There are ways to get the differences directly; I'm not sure how one does this with the software you are using (unstated) but it is reasonably trivial (if tedious) to compute them by hand. I have an example of how to do this here: https://fromthebottomoftheheap.net/2017/10/10/difference-splines-i/ and I'm reasonably sure that plotDiff is following the same approach to compute the differences between smooths. As such you could also look at the code behind plotDiff to see what it is doing.
Once you have the difference and the confidence interval, it only requires you to identify which points where the CI excludes 0. Find the points where 0 is inside the interval and then invert that selection: !(upr_ci > 0 & lwr_ci < 0) for example.
