How to properly smooth a 2d map? Lately I've been trying to model the probability of success for an event given it's location. Basically, the data I have consists of locations in x and y coordinates, and a corresponding value of 1 or 0 (indicating success or failure). When I bin the events into boxes and calculate their average success rate for that bin, I get:

My goal is to smooth this out, as there is obviously some noise in here. I've tried multiple methods (nonlinear regression, logit model, mgcv package) of which none seemed to give a proper smooth. For example, a full tensor product smooth from the binomial family (using the mgcv package), yielded this result:

This is kind of in the right direction, but it does seem to underfit the very high probabilities around the 'white' area. Upping the degrees of freedom (currently 15) in the model is limited due to computatoinal complexity.
I've been trying to get a more accurate model of the probabilities for some weeks without success. Any tips/ideas on how to improve this? (I'm using R at the moment.)
 A: I do wonder if there is a big problem here at all? If you look at the very white pixels in your upper figure, they are surrounded by points that are less intensely-coloured (?) white. What you have here are data and the fitted spline is balancing the highly successful point with the less successful ones in the same spatial location.
A basis dimensionality of 15 sounds far too few for a 2-d thin plate or tensor product spline. In this are you referring to the k in s(x, y, k = 15) or te(x, y, k = 15), or did you set k to something else? If the former, then you actually have 15*15 as the basis dimension, which would seem more reasonable as a starting point: the smoothness selection will reduce this via the wiggliness penalty.
As you have both known 0 and 1's, I would start, if using mgcv, using the actual coordinates of the shots (not aggregated to pixels) with:
m <- gam(scored ~ te(x, y), data = myDF, family = binomial, method = "REML")

and then check the resulting fit for adequacy, especially via gam.check(). This will tell you whether the basis dimension for the te() term was high enough or not (look for K close to 1 or a high p value in the printed output). You can also assess the model diagnostic plots but some of these may be less than useful as the response is binary.
If gam.check() suggests too low a basis dimension (the default is $2^5$) and you should increase k.
It might be that you need an adaptive smoother: you can do te(x, y, bs = "ad") but this will impact the computational efficiency a lot as this is effectively going to fit several different GAMs to section of the data.
If you are on Linux you could consider setting the nthreads option via gam.control to a number equal to the available CPU cores (or thereabouts) in the gam() call: control = gam.control(nthreads = 4). This will utilise multiple threads in some parts of the computation and can for some models greatly decrease the compute time.
You could also look at bam() as a drop-in replacement for gam(), which is designed to work on big data sets.
