Show density plot on a geographic map

Question

I have a bunch of data points in the form of {latitude, longitude}. How can I plot their density on a US map? Single points on the map do not clearly show their distribution, so I'd like it to be a density map or something clearer. Please suggest.

---

Note from Vitaliy Kaurov: this subject is also applicable to any geographic map (not necessarily US) and also to contour plots and general heat maps (besides density plots).

Answer 1

Update (15 FEB 2017)

In the light of new geo-functionality developments I would like to add another approach via GeoImage. Consider the magnitude of the geomagnetic field across the entire planet. Get the data:

data = QuantityMagnitude[
   GeomagneticModelData[{{-90, -180.}, {90, 180.}}, "Magnitude", 
    GeoZoomLevel -> -2]];

Build high-res contour plot:

plot = ListContourPlot[data, PlotRangePadding -> 0, 
  ColorFunction -> "TemperatureMap", Frame -> False, 
  InterpolationOrder -> 2, ImageSize -> 800, Contours -> 20, 
  ContourStyle -> Opacity[.2]]

Specify how you would like to place this rectangle on the map (we use world range to show the application of geo-projection too):

GeoGraphics[{GeoStyling[{"GeoImage", plot}, 
   GeoRange -> {{-90, 90}, {-180, 180}}], Opacity[0.6], 
  GeoBoundsRegion["World"]}, GeoRange -> "World", 
 GeoProjection -> "Robinson"]

Colored background interferes with colors of the heat map. Use just coastlines:

GeoGraphics[{
  GeoStyling[{"GeoImage", plot}, 
   GeoRange -> {{-90, 90}, {-180, 180}}], Opacity[.8], 
  GeoBoundsRegion["World"]}, GeoRange -> "World", 
 GeoProjection -> "Robinson", GeoBackground -> 
{"Coastlines", "Land" -> White, "Ocean" -> White, "Border" -> Black}]

Older versions

Well I do not have your data, so I'll just show things I have. I'll go with India, it has nice cities layout. You can upgrade for your data yourself. If you know only locations of cities, but not their population, then

data = DeleteCases[
   Reverse[CityData[#, "Coordinates"]] & /@ CityData[{All, "India"}], 
Missing["NotAvailable"]];

Show[SmoothDensityHistogram[data, .3, ColorFunction -> "SunsetColors",
   PlotPoints -> 150, Mesh -> 20, MeshStyle -> Opacity[.1], PlotRange -> All], 
 Graphics[{FaceForm[], EdgeForm[Directive[White, Opacity[.7]]], 
   CountryData["India", "Polygon"]}]]

If you know the population, then maybe something like this:

Graphics[{Opacity[0.1, Blue], 
  Cases[Disk[Reverse[CityData[#, "Coordinates"]], 
      Log[10.^12, CityData[#, "Population"]]] & /@ 
    CityData[{All, "India"}], Disk[{__Real}, _Real]], {FaceForm[], 
   EdgeForm[Directive[Red, Thick, Opacity[.9]]], 
   CountryData["India", "Polygon"]}}, Background -> Black]

For US use something like:

data = DeleteCases[Reverse[CityData[#, "Coordinates"]] & /@ 
    CityData[{All, "UnitedStates"}], Missing["NotAvailable"]];
Show[SmoothDensityHistogram[data, .5, ColorFunction -> "SunsetColors",
   PlotPoints -> 150, Mesh -> 20, MeshStyle -> Opacity[.1], PlotRange -> All], 
 Graphics[{FaceForm[], EdgeForm[Directive[White, Opacity[.7]]], 
   CountryData["UnitedStates", "FullPolygon"]}], 
 AspectRatio -> Automatic, PlotRange -> {{-180, -65}, {20, 70}}, 
 ImageSize -> 800, Frame -> False]

Answer 2

Due to the limitations on the accuracy of SmoothKernelDistribution (mentioned by Murta) upon which SmoothDensityHistogram is based, I prefer to work with the more exact KernelMixtureDistribution.

I will use the same data as Vitaliy here.

data = DeleteCases[
   Reverse[CityData[#, "Coordinates"]] & /@ CityData[{All, "India"}], 
   Missing["NotAvailable"]];

Note the use of MaxMixtureKernels -> All. This ensures that each data point will have a kernel placed at it. If we don't do this KernelMixtureDistribution may use a binned estimator.

dens = KernelMixtureDistribution[data, MaxMixtureKernels -> All];

For some added value credit goes to @rm-rf for pointing out the following point-in polygon algorithm. Note that it does not work for polygons with multiple components so we will need to table over those.

inPolyQ[poly_, pt_] := Graphics`Mesh`InPolygonQ[poly, pt];

Lets obtain the polygon data for India...

poly = CountryData["India", "Polygon"][[1]];

And now to display the density (I prefer the use of ContourPlot but we could just as easily have used DensityPlot). I use Table here so that we account for the small islands to the south east of the mainland.

p1 = ContourPlot[Evaluate@PDF[dens, {x, y}], {x, 65, 100}, {y, 0, 40},
       PlotRange -> All, ColorFunction -> "TemperatureMap", 
       PlotPoints -> 100, RegionFunction -> (inPolyQ[poly[[1]], {#1, #2}] &), 
       PlotLegends -> Automatic];

Show[p1, Table[
       ContourPlot[Evaluate@PDF[dens, {x, y}], {x, 65, 100}, {y, 0, 40}, 
         PlotRange -> All, ColorFunction -> "TemperatureMap", 
         PlotPoints -> 100, RegionFunction ->(inPolyQ[i, {#1, #2}] &)], {i,Rest@poly}], 
       Graphics[{FaceForm[], EdgeForm[Black], CountryData["India", "Polygon"]}]]

Answer 3

Here's just a different suggestion. Playing around with the parameters can yield very different results, so it's encouraged.

I'm using the same data as Vitaliy:

data = DeleteCases[Reverse[CityData[#, "Coordinates"]] & /@ CityData[{All, "India"}], Missing["NotAvailable"]];
outline = ColorNegate@Graphics[CountryData["India", "Polygon"]];

My solution plots each point as a disk and then blurs them together, hence creating a heat map.

density[data_, r_] := Blur[Graphics[Disk[#, 0.1] & /@ data], r];
map[outline_, density_, col_] := ImageCrop@ImageMultiply[ ImageApply[(ColorData[col]@Mean[#] &@#) /. RGBColor[r_, g_, b_] :> {r, g, b} &, density], outline];

Examples:

map[outline, density[data, 20], "SunsetColors"]

map[outline, density[data, 5], "AvocadoColors"]

If you want a different background you can replace ImageMultiply by SetAlphaChannel and then use for example Show[myMap,Background->Blue]. It is possible to add a border using EdgeDetect, which can be put on top of the image. Dilation can be used to adapt the width of the border.

EDIT

Creating these charts seem to require some manual changing of the plot range for certain polygons (depending on how CountryData specifies them, I suppose). I honestly thought the above would always work. I tried making a general approach by setting PlotRange and PlotRangePadding explicitly, and this worked but the solution is not very short anymore (though this partly depends on a few new options, like borders, that I added).

Options[densityMap] = {diskSize -> 1, blurRadius -> 10, padding -> 10,
    background -> Black, showBorder -> False, borderMagnify -> 0, 
   borderColor -> Black};
densityMap[country_, data_, col_, OptionsPattern[]] := 
 Module[{plotrange, disks, density, outline, outlineObj, border, 
   rendered},
  outlineObj = 
   Graphics[{White, CountryData[country, "Polygon"]}, 
    Background -> OptionValue[background], 
    PlotRangePadding -> OptionValue[padding]];
  outline = Rasterize[outlineObj];
  plotrange = PlotRange /. AbsoluteOptions[outlineObj];
  disks = Graphics[{
     Disk[#, OptionValue[diskSize]] & /@ data
     }, PlotRange -> plotrange, 
    PlotRangePadding -> OptionValue[padding]];
  density = SetAlphaChannel[
    SetAlphaChannel[
     ImageApply[(ColorData[col]@Mean[#] &@#) /. 
        RGBColor[r_, g_, b_] :> {r, g, b} &, 
      ColorNegate@Blur[disks, OptionValue[blurRadius]]], 
     ColorNegate@Blur[disks, OptionValue[blurRadius]]],
    Graphics[{White, CountryData[country, "Polygon"]}, 
     Background -> Black, PlotRangePadding -> OptionValue[padding]]
    ];

  If[OptionValue[showBorder],
   border = 
    Dilation[
     EdgeDetect[
      Graphics[{White, CountryData["United States", "Polygon"]}, 
       Background -> Black, 
       PlotRangePadding -> OptionValue[padding]]], 
     OptionValue[borderMagnify]];
   rendered = 
    SetAlphaChannel[
     ImageAdd[ColorConvert[ColorNegate@border, "RGB"], 
      OptionValue[borderColor] /. RGBColor[r_, g_, b_] :> {r, g, b}], 
     border];
   Show[outline, density, rendered],
   Show[outline, density]
   ]

  ]

Using the OP suggested:

data = Reverse /@ {CityData[{"New York City", "USA"}, "Coordinates"], 
    CityData[{"Boston", "Massachusetts", "USA"}, "Coordinates"]} ;

densityMap["United States", data, "SunsetColors", blurRadius -> 5, 
 diskSize -> 1, background -> White]

Unfortunately this method introduces white artifacts, which I believe are due to a bug in ImageApply. For that reason this function works best with either a white background or one can use the border option which will cover them up. See the beginning of the code for the list of options.

Answer 4

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, I illustrated the "spherical smooth histogram" function from this paper; 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.