I try to model beak length in tropical birds by sampling year and age of the birds to see if beaks changed in size over the years using GAMMs. The data was collected for a number of birds (usually 20 birds per sampling site and sampling year) at 6 different sampling sites over 25 years. Unfortunately, the sampling was not done consistently so that I have two locations which have only been sampled once in the 25 years. I have no reason to believe that beak length is severely impacted by sampling site but in my model I included sampling site as random effect as in: s(sampling_site, bs="re"). However, I am not sure if this really makes sense considering that two locations have only been sampled once. Any advise would be much appreciated.
From what you describe, this should be fine. In fact, a main advantage of random effects is that they induce shrinkage (a nice illustration an be found here), which has also been called partial pooling: Information from the other sites is used to estimate the random intercept for the site with only one observation, because you are modelling them as coming from the same normal distribution (actually, it's not even clear from your description whether there is only one observation for a certain site, or if there are mutliple, but all from the same sampling date).
Due to this, you will probably have no model identifiability or convergence issues. But be aware that estimates of random intercepts for sites with one or few observations are pulled more strongly towards the overall mean than those with many. In case you are treating site as a mere nuisance variable, though, as typical in ecological studies, this should be irrelevant, because you are not going to interpret site effects.
Note that, as in any mixed model, you should ask yourself whether it is plausible that the site-specific effect is not correlated with your predictor of interest, i.e., bird age (the "random effects assumption"). This becomes especially relevant under stronger shrinkage.