Splines free of Gibbs phenomenon for smoothing I'm trying to do 3-dimensional smoothing of some observations we made. To simplify things, here's a visualisation of one slice of the data we have:

Note that the observations are gridded and contain a lot of zeros. When I run a cubic smoothing spline (I use the package csaps in python) we get something like the plot below.

(My apologies for the changing colour scales).
As can be seen on the second plot, the splines oscillate and become even negative, which in our application does not make sense and is an issue. I researched a bit and found that this phenomenon is often referred to as Gibbs phenomenon that can affect splines (for example in this or that application).
I have a vague understanding of the math behind it and am more interested in the practical application. Does anybody know of an implementation of splines that are free of the Gibbs phenomenon in Python (or R or something else really) that we could use for our application?
Some more background on the application:
We are estimating the probability that two particles are in the same volume. The observations are occupation rates per volume (in some arbitrary unit, such as hours/year). We would like to fit the smooth to the data in an analogous way as you can fit a function to a 1D histogram (in our case it’s just 3D). Since the observed data is noisy (and we don't have access to the underlying generating process), we would like to use the fitted function to draw or to estimate some probabilities instead of using the observations directly. The probabilities (or occupancy rate, which is somewhat similar in our case) are then used to estimate the probability that two particles are in the same volume (the other particle is generated by a different process which we know how to model).
 A: Thanks to the hints in the comments, I managed to find a suitable solution to my problem. More precisely, I found two that overcome the issue with the splines. But neither of my solutions uses splines anymore.
I think it's worth highlighting that the suggestion of @Gavin is probably also possible, but I didn't try this since I'm quite happy with what I have so far.
Gaussian Filter
I applied a 3D gaussian filter to the data and the results look reasonable, as shown in the figure below.

From my point of view, the drawback of this method is that the choice of sigma (standard deviation of the Gaussian kernel) is somewhat arbitrary. Instead, I settled on the next method.
Kernel Density Estimation
Thanks to the hint from @Aksakal I looked into kernel density estimation (KDE) and it seems to me suitable for the problem at hand. The result is shown in the figure below.

What I like about this method is that the bandwidth of the kernel can be estimated from the data with cross validation (in my opinion this is an advantage over just simply going for the value that seems to be looking fine).
If anybody is interested in KDE in Python, an excellent comparison of the different implementations to be found here.
