Here is another possibility. The trouble with the use of either `SmoothKernelDistribution[]` or `KernelMixtureDistribution[]` is that neither of the two distributions take into account that the data are, well, points on a sphere! (Fine, more or less a sphere... :P)
In this [previous answer](https://mathematica.stackexchange.com/a/14020), I illustrated the "spherical smooth histogram" function from [this paper](http://dx.doi.org/10.1016/0098-3004(85)90015-9); among other things, the function will not suffer from boundary problems if, say, one wants to study the islands clustered around the International Date Line.
Only a few modifications are needed so that the density function is expressed in terms of longitude/latitude:
(* U.S. outline *)
uspoly = Reverse /@ First[CountryData["UnitedStates", "Coordinates"]];
(* some cities *)
uscities = Select[Reverse[CityData[#, "Coordinates"]] & /@
CityData[{Large, "UnitedStates"}], FreeQ[#, _Missing] &];
(* unit vectors *)
cityVecs = {Cos[#1] Cos[#2], Sin[#1] Cos[#2], Sin[#2]} & @@@ (uscities °);
(* "robust" estimate of smoothing parameter *)
c1 = (Length[cityVecs]^(1/3) (Length[cityVecs] - 1))/
Total[ArrayPad[Drop[Sort[1 - cityVecs.Normalize[Total[cityVecs]]], -2],
{0, 2}, "Fixed"]];
(* smoothing parameter, automatically determined by maximizing "pseudo-log likelihood" *)
(* as a conservative estimate, the larger of the "pseudo-log likelihood" estimate
and "robust" estimate is chosen *)
c = Max[c1,
With[{n = Length[cityVecs]}, First @ FindArgMax[
Sum[Log[Total[\[FormalC] Csch[\[FormalC]] Exp[\[FormalC]
Delete[cityVecs, k].Extract[cityVecs, k]]/
(4 π (n - 1))]], {k, n}], {\[FormalC], c1},
Method -> "PrincipalAxis"]]];
(* spherical smooth density histogram *)
With[{n = Length[cityVecs]},
DensityPlot[Total[c Csch[c] Exp[c cityVecs.
{Cos[λ °] Cos[ϕ °], Sin[λ °] Cos[ϕ °], Sin[ϕ °]}]/(4 π n)],
{λ, -130, -60}, {ϕ, 20, 50},
AspectRatio -> Automatic, ColorFunction -> "ThermometerColors",
Mesh -> True, MeshFunctions -> {#3 &}, PerformanceGoal -> "Quality",
PlotPoints -> 95, PlotRange -> All,
RegionFunction -> (Graphics`PolygonUtils`InPolygonQ[uspoly,
{#1, #2}] &)]]
---
As another demonstration, here is a smooth histogram of the capital cities of the world:
capitals = Select[Reverse[CityData[CountryData[#, "CapitalCity"], "Coordinates"]] & /@
CountryData[], FreeQ[#, _Missing] &];
cVecs = {Cos[#1] Cos[#2], Sin[#1] Cos[#2], Sin[#2]} & @@@ (capitals °);
(* smoothing parameter *)
c1 = (Length[cVecs]^(1/3) (Length[cVecs] - 1))/
Total[ArrayPad[Drop[Sort[1 - cVecs.Normalize[Total[cVecs]]], -2],
{0, 2}, "Fixed"]];
c = Max[c1,
With[{n = Length[cVecs]}, First @ FindArgMax[
Sum[Log[Total[\[FormalC] Csch[\[FormalC]] Exp[\[FormalC] Delete[
cVecs, k].Extract[cVecs, k]]/(4 π (n - 1))]], {k, n}],
{\[FormalC], c1}, Method -> "PrincipalAxis"]]];
hist = With[{n = Length[cVecs]}, Image[DensityPlot[
Total[c Csch[c] Exp[c cVecs.
{Cos[λ °] Cos[ϕ °], Sin[λ °] Cos[ϕ °], Sin[ϕ °]}]/(4 π n)],
{λ, -180, 180}, {ϕ, -90, 90}, AspectRatio -> Automatic,
ColorFunction -> "ThermometerColors", Frame -> False, Mesh -> True,
MeshFunctions -> {#3 &}, ImagePadding -> None, PerformanceGoal -> "Quality",
PlotPoints -> 95, PlotRange -> All, PlotRangePadding -> None],
ImageResolution -> 300]];
ParametricPlot3D[{Cos[λ °] Cos[ϕ °], Sin[λ °] Cos[ϕ °], Sin[ϕ °]},
{λ, -180, 180}, {ϕ, -90, 90}, Axes -> None, Boxed -> False,
Lighting -> "Neutral", Mesh -> None, PlotStyle -> Texture[hist],
TextureCoordinateFunction -> ({#4, #5} &)]
---
For some reason, the performance of this method is rather slow for Vitaliy's India example; I'll do more research on how to make this method more efficient.