2
$\begingroup$

I have latitude and longitude locations that I need density (or contour) plot to overlay on a map. For example:

dat = {{40.5`, -111.35`}, {37.84`, -112.87`}, {38.79`, -112.05`}, 
{40.28`, -109.23`}, {38.79`, -112.04`}, {41.42`, -112.87`}, {37.09`, 
-112.9`}, {37.82`, -113.1`}, {42, -112.6`}, {39.58`, -111.69`}, 
{40.69`, -109.31`}, {38.074`, -112.192`}, {38.71`, -112.43`}, 
{38.73`, -112.44`}, {38.73`, -112.45`}, {38.72`, -112.46`}, {38.71`, 
-112.44`}, {38.73`, -112.44`}, {41.42`, -112.87`}, {38.88`, 
-111.98`}, {38.86`, -111.97`}, {38.89`, -111.98`}, {39.75`, 
-110.84`}, {40.61`, -109.41`}};
ad = AdministrativeDivisionData[{"Utah", "UnitedStates"}];

GeoGraphics[{Gray, ad["Polygon"], Red, PointSize[.01], 
  Point[GeoPosition /@ Rest@dat]}, GeoBackground -> None]

enter image description here

I create the density plot.

sdh = SmoothDensityHistogram[Reverse[dat, 2], PlotRange -> Full, 
  PlotRangePadding -> None, Frame -> False, ImagePadding -> None]

enter image description here

However, when I overlay it on the geographic region the density plot is not aligned with the geographic locations.

GeoGraphics[{
  {Gray, GeoStyling[{"GeoImage", sdh}], ad["Polygon"]},
  {Red, PointSize[.01], Point[GeoPosition /@ Rest@dat]}},
 GeoBackground -> None]

enter image description here

What am I doing wrong? Should I be using something other than "GeoImage"?

$\endgroup$
  • $\begingroup$ Are you looking for something like this? To get it, I extracted the polygon and used it as the RegionFunction for the histogram plot $\endgroup$ – Jason B. Jul 12 '16 at 18:56
  • $\begingroup$ @JasonB Yes. That is what I am looking for. Although it is not a GeoGraphics it may be a good workaround if nothing else materializes. Please post as an answer. $\endgroup$ – Edmund Jul 12 '16 at 19:01
  • $\begingroup$ @Edmund Probably you can use the method shown in this WC post: "Mass shootings and availability of gun dealerships." $\endgroup$ – Alexey Popkov Jul 12 '16 at 19:59
  • 2
    $\begingroup$ I think your code is basically correct. Just change PlotRange->Full to PlotRange -> Reverse[GeoBounds[ad]] in the definition of sdh. $\endgroup$ – jose Jul 13 '16 at 12:12
  • 1
    $\begingroup$ @jose - Nice catch, I just stole that and put it in my answer below $\endgroup$ – Jason B. Jul 13 '16 at 13:37
4
$\begingroup$

jose gave the best answer this question in comments, so I will add it here. As plotted, the PlotRange of the density plot is {{-113.032, -109.298}, {36.3135, 42.7765}}, but the GeoBounds for the map is {{-114.051, -109.045}, {36.9991, 42.0016}}. Apparently this mismatch is the problem.

sdh = SmoothDensityHistogram[Reverse[dat, 2], 
   PlotRangePadding -> None, Frame -> False, ImagePadding -> None, 
   PlotRange -> Reverse[GeoBounds[ad]]];
GeoGraphics[{{Gray, GeoStyling[{"GeoImage", sdh}], 
   ad["Polygon"]}, {Red, PointSize[.01], Point[GeoPosition /@ dat]}}, 
 GeoBackground -> None]

Mathematica graphics

I give a couple of alternate methods below as well.

Workaround 1 - Graphics using the geo-polygon

I remembered doing something similar here,

utah = DiscretizeGraphics[
   AdministrativeDivisionData[{"Utah", "UnitedStates"}, "Polygon"]];
SmoothDensityHistogram[
 Reverse[dat, 2],
 PlotRange -> Full,
 PlotRangePadding -> None,
 Frame -> False,
 ImagePadding -> None,
 RegionFunction -> Function[{x, y},
   {x, y} ∈ utah],
 AspectRatio -> Automatic,
 Epilog -> {Red, PointSize[.01], Point@*Reverse /@ dat}]

Mathematica graphics

Workaround 2 - Extract the Graphics primitives and feed them to GeoGraphics

This is actually my favorite method. Here we take the trick from the of forming the plot using a RegionFunction, and combine it with Mr. Wizard's method for inserting Graphics into GeoGraphics:

reg1 = DiscretizeGraphics[
   AdministrativeDivisionData[{"Utah", "UnitedStates"}, "Polygon"]];
sdh = SmoothDensityHistogram[Reverse[dat, 2], PlotRange -> Full, 
   PlotRangePadding -> None, Frame -> False, ImagePadding -> None, 
   RegionFunction -> Function[{x, y}, {x, y} \[Element] reg1]];
prim = First@Cases[sdh, Graphics[a_, ___] :> a, {0, -1}, 1];
GeoGraphics[{prim, Red, PointSize[.01], 
  Point[GeoPosition /@ Rest@dat]}, GeoBackground -> None]

Mathematica graphics

$\endgroup$
  • $\begingroup$ The regions lose their proportions when this method is used. They don't look like their GeoGraphics counterparts. This one is ok because it is essentially a box but others do not far as well. Also does not handle complicated borders very well. Thanks for the effort, though. (+1) $\endgroup$ – Edmund Jul 12 '16 at 19:26
  • $\begingroup$ Nice wrap-up of the workaround and the answer. $\endgroup$ – Edmund Jul 14 '16 at 13:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.