We propose a novel model selection algorithm based on a penalized maximum likelihood estimator (PMLE) for functional hidden dynamic geostatistical models (f-HDGM). These models employ a classic mixed-effect regression structure with embedded spatiotemporal dynamics to model georeferenced data observed in a functional domain. Thus, the parameters of interest are functions across this domain. The algorithm simultaneously selects the relevant spline basis functions and regressors that are used to model the fixed-effects relationship between the response variable and the covariates. In this way, it automatically shrinks to zero irrelevant parts of the functional coefficients or the entire effect of irrelevant regressors. The algorithm is based on iterative optimisation and uses an adaptive least absolute shrinkage and selector operator (LASSO) penalty function, wherein the weights are obtained by the unpenalised f-HDGM maximum-likelihood estimators. The computational burden of maximisation is drastically reduced by a local quadratic approximation of the likelihood. Through a Monte Carlo simulation study, we analysed the performance of the algorithm under different scenarios, including strong correlations among the regressors. We showed that the penalised estimator outperformed the unpenalised estimator in all the cases we considered. We applied the algorithm to a real case study in which the recording of the hourly nitrogen dioxide concentrations in the Lombardy region in Italy was modelled as a functional process with several weather and land cover covariates.