This paper describes the specification of spatially varying coefficient (SVC) regression models using Generalized Additive Models (GAMs). The GAMs include Gaussian Process (GP) splines (smooths) for each covariate parameterised with observation location and generate SVC estimates that capture spatially varying relationships. The ability of the proposed GAM approach to estimate true spatially varying coefficients was compared with that of Multiscale Geographically Weighted Regression (MGWR) using simulated data with complex spatial heterogeneities. The geographical GP GAM (GGP-GAM) was found to out-perform MGWR across a range of fit metrics and resulted in more accurate coefficient estimates and lower residual errors. The model for one of simulated datasets was investigated in detail to illustrate GAM diagnostics, model checks, spline / smooth convergence and basis evaluations, and tuning via the number knots. A larger simulated case study was investigated to explore the trade-offs between GGP-GAM complexity, performance and computation. Finally the GGP-GAM and MGWR approaches were applied to an empirical case study. The resulting models had very similar accuracies and fits, and generated subtly different spatially varying coefficient estimates. A number of areas of further work are identified.