I have gone through your posted code and have a couple of comments/ideas here: 

 1. The boundary object that you're supplying does not have the variable `f` in each boundary list. This variable tells the model what to do right at the boundary. You will want to set `f` to `0` considering you want to model occurrence in the lake. Please the [vignette][1]/[blog post][2] on when in the code and how to create `f`. 
 2. I have created the following [GitHub repo][3] with my edits to the code. I still cannot get [{gratia}][4] to draw the model, as it produces the same error, but you should be able to predict densities across the lake by following the blog/the GitHub repo uses the same methodology
 3. To check fixed and random effects and smoothers. I suggest looking at the main effects of the model first prior to running `summary()`. This can be done using `anova.gam()` or `anova()`. You will notice for the soap-film model that `origin` is not statistically significant which for the thin-plate model it is.
 4. I did not change this in my code but I do agree with @robert-lennox that your `k` value is likely to high resulting in the model overfitting. 

  


 


  [1]: https://github.com/dill/soapcheckr
  [2]: https://blog.benjaminhlina.com/posts/post-with-code/soapcheckr/
  [3]: https://github.com/benjaminhlina/soap-film-help
  [4]: https://gavinsimpson.github.io/gratia/