I’ve worked with similar kinds of data and have used tensor product splines to estimate the space-time components. You can use the t2() function in the formula to set up these kinds of terms. For example a smooth anisotropic spatial smooth interacting with a time component would be set up as
t2(x, y, time, k = c(10,10,15))
where x and y are the spatial coordinates, time is something like the observation year, and k is the basis size for each marginal smooth. An isotropic version would be
t2(x, y, time, d = c(2,1), k = c(50,15))
where the d sets the x and y coordinates to be from a 2d thin-plate spline, which has a single maximum wiggliness (here set to say 50).
The t2() in brms uses mgcv to set up the smoothers so as long as it goes into a t2() you can use any of a wide range of smoother types in the brms model.
You mention polygons (as if the data are areal observations?); if so, then this kind of spatial data can be accommodated in the spline framework via a Markov random field using the bs = 'mrf' option in mgcv and hence brms. You need to pass t2() a factor indexing the region/polygon for each observation rather than x,y pairs, but you can do this via a neighbourhood object generated from a polygon shapefile object read into R. See ?mgcv:::mrf for the detail.
These are not the kinds of models you mentioned but they are somewhat similar.