Smoothing in Time and Space, with Random Effects Setup
My city has recently run a pilot program to mandate that only public transportation use particular sections of the downtown core. 
The city has released data about the program.  They have sensors at strategic locations across the downtown core to monitor the travel time between segments.  The data include:


*

*The latitude and longitude of the sensor

*A code for the segment (e.g. Yo_Ds_Yo_Ds).

*The time the measurements were made.  

*The is the average number of seconds required to travel through a segment.  This is a measure of congestion.  More time means higher congestion.


The data span 3 months before the pilot began and nearly 8 months after the pilot started.  Data is recorded for all hours of the day.
I'm interested in determining how congestion has changed in surrounding areas due to the program coming into effect. I anticipate this will require a smooth in space and time.  
The Data and Model
The dataset has 5 columns:


*

*TravelTime -- Avg time in seconds to traverse the segment

*Lat/Lon -- coordinates of the sensor

*Time -- Time the observation was made

*IsProgram -- Binary indicator for if the observation is from the program.


As of now, I use the following to model the data
bam(TravelTime ~ s(Lat,Lon, bs = 'gp', by = IsProgram) + s(Time,bs = 'cc'))
here is my rationale:


*

*I use a bam because I have something like 1,000,000+ observations

*IsProgram is a factor indicating if the observation is from when the program was in place.  Passing this into the by argument will allow me to estimate congestion for pre and post program.

*I use a s(Time) to control for the effect of time on congestion.  The cc is because the effect of time is periodic. 
Question
Is this a reasonable approach to using GAMS to model this problem?  Should I be using a random effects model since I have repeated observations from sensors?
If you need data, I can link to a GitHub repo.
 A: In general, you don't need to include a sensor random effect if you think that the reason sensor $s_i$ records higher transit times on average than sensor $s_j$ is due to a smoothly varying spatial process.
All the random effect is doing is saying that the mean of the response at each sensor is a random variable with mean equal to the model intercept plus a mean zero random variable with some variance to be estimated.
All the spatial smooth is saying it that the mean of the response at each sensor is a random variable with mean equal to the model intercept plus a contribution from a smoothly varying spatial process.
That spatial process need not be smooth but it may have more spatial structure than implied by the random effect. For example, we might consider the sensors to be yielding data for an entire block of space and then we might want to have a spatial process that uses areal information (think data aggregated at a state or county level, where there is spatial structure but it is not in the form necessarily of continuous spatial coordinates because the data are aggregated areally at the state or county level). Then we might model the spatial process as an Markov Random Field, which can also be represented as a penalised spline basis. This would be in the GAM framework. There are other ways to approach this kind of spatial modelling.
I would suggest not fitting your model as you indicated in your question. You are making a very strong assumption that the length scale of the correlation functions that comprise the Gaussian process basis is the largest distance between any two locations. For the equivalence of the Gaussian process and the penalized spline basis / quadratic penalty models estimated by mgcv, the GP basis cannot depend on additional parameters and yet the correlation functions do depend on the length scale parameter. Hence an heuristic is use by default if the user doens't supply this via the m argument. That default is to use the maximum distance between any two points as the length scale.
To estimate this model without fixing the length scale, you need to write a loop over values for the length scale and plug that into the model, estimate model, and then extract the REML score from the fit. Then you can plot the REML score against the length scale values you fitted at and see if there is a minimum in the curve produced, which would suggest an optimum length scale for this model.
More simply, just replace the gp baiss with the default tp basis (thinplate regression spline) or even the ds basis for a Duchon spline and you'll invariably do just as well without having to profile over the length scale parameter.
