I wouldn't use Tweedie for this given what is shown. Species richness should be an integer value and the output is showing that you are effectively fitting a Poisson GAM (the $p = 1.01). I've rarely had good diagnostic plots when modelling species richness data with a Poisson model. Instead try the
nb() family for the negative binomial; that should help some with the extra-Poisson dispersion evident in the QQ plot.
I assume the response is not species richness of those birds that you computed the PCA on?
In general, I'm not a fan of using PC axes like this in models, mainly because you're going to run into problems explaining what the smooth of the two axes means in terms of the original data (bird abundances). I would also suggest that this smooth should be a tensor product, not a 2d thin plate spline, as the latter assumes the same wiggliness in all directions, which need not be the case with PCs.
As for the $k$ issue; the heuristic basis of the test starts to run into problems when your data are in the form of a time series or spatial layout. Even small amounts of autocorrelation left over in the residuals will cause the test to reject the null hypothesis (in my experience). Simon Wood, in his GAM book (2nd ed.) suggests an alternative:
take the deviance residuals from your model,
add them as the response in a GAM with the same covariate terms, the following family
family = quasi(link = "identity", variance = "constant")
and double the
k on the smooth terms,
look at the model summary and see if the EDF for the smooth terms in greater than 1 (the minimum),
Go back to your original model and double
k for any term that had a greater than 1 EDF in the model for the deviance residuals,
rinse and repeat if necessary
I probably wouldn't use the the shrinkage smooth for the spatial effect as my starting point, get the main features figured out first, then switch on extra penalties if you need them.
"mrf" basis will also work for this kind of spatial grid layout, but you'll either run into computational problems (if you run the full MRF you'll need a 89,000 x 89,000 penalty matrix) or you'll need to set
k lower to get the low-rank version. But give it a go if you want to try it. As the MRF has a different definition of wiggliness, it might work better in the spatial setting, but it will only do so by being very wiggly. Whether a smoother fit or a wigglier fit is desirable isn't clear from your question.