You might consider using a kernel especially suitable for the sphere, such as a von Mises-Fisher density
$$f(\mathbf{x};\kappa,\mu) \propto \exp(\kappa \mu^\prime \mathbf{x})$$
where $\mu$ and $\mathbf{x}$ are locations on the unit sphere expressed in 3D Cartesian coordinates.
The analog of the bandwidth is the parameter $\kappa$. The contribution to a location $x$ from an input point at location $\mu$ on the sphere, having weight $\omega(\mu)$, therefore is
$$\omega(\mu) f(\mathbf{x};\kappa,\mu).$$
For each $\mathbf{x}$, sum these contributions over all input points $\mu_i$.
To illustrate, here is R code to compute the von Mises-Fisher density, generate some random locations $\mu_i$ and weights $\omega(\mu_i)$ (12 of them in the code), and display a map of the resulting kernel density for a specified value of $\kappa$ (equal to $6$ in the code).
![[Figure]](https://i.sstatic.net/4Z2fg.png)
The points $\mu_i$ are shown as black dots sized to have areas proportional to their weights $\omega(\mu_i)$. The contribution of the large dot near $(100,60)$ is evident throughout the northern latitudes. The bright yellow-white patch around it would be approximately circular when shown in a suitable projection, such as an Orthographic (earth from space).
# # von Mises-Fisher density. # mu is the location and x the point of evaluation, *each in lon-lat* coordinates. # Optionally, x is a two-column array. # dvonMises <- function(x, mu, kappa, inDegrees=TRUE) { lambda <- ifelse(inDegrees, pi/180, 1) SphereToCartesian <- function(x) { x <- matrix(x, ncol=2) t(apply(x, 1, function(y) c(cos(y[2])*c(cos(y[1]), sin(y[1])), sin(y[2])))) } x <- SphereToCartesian(x * lambda) mu <- matrix(SphereToCartesian(mu * lambda), ncol=1) c.kappa <- kappa / (2*pi*(exp(kappa) - exp(-kappa))) c.kappa * exp(kappa * x %*% mu) } # # Define a grid on which to compute the kernel density estimate. # x.coord <- seq(-180, 180, by=2) y.coord <- seq(-90, 90, by=1) x <- as.matrix(expand.grid(lon=x.coord, lat=y.coord)) # # Give the locations. # n <- 12 set.seed(17) mu <- cbind(runif(n, -180, 180), asin(runif(n, -1, 1))*180/pi) # # Weight them. # weights <- rexp(n) # # Compute the kernel density. # kappa <- 6 z <- numeric(nrow(x)) for (i in 1:nrow(mu)) { z <- z + weights[i] * dvonMises(x, mu[i, ], kappa) } z <- matrix(z, nrow=length(x.coord)) # # Plot the result. # image(x.coord, y.coord, z, xlab="Longitude", ylab="Latitude") points(mu[, 1], mu[, 2], pch=16, cex=sqrt(weights))