# How to align GeoGraphics “GeoImage” with GeoPositions?

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}};

Point[GeoPosition /@ Rest@dat]}, GeoBackground -> None]


I create the density plot.

sdh = SmoothDensityHistogram[Reverse[dat, 2], PlotRange -> Full,


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

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


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

• Are you looking for something like this? To get it, I extracted the polygon and used it as the RegionFunction for the histogram plot – Jason B. Jul 12 '16 at 18:56
• @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. – Edmund Jul 12 '16 at 19:01
• @Edmund Probably you can use the method shown in this WC post: "Mass shootings and availability of gun dealerships." – Alexey Popkov Jul 12 '16 at 19:59
• I think your code is basically correct. Just change PlotRange->Full to PlotRange -> Reverse[GeoBounds[ad]] in the definition of sdh. – jose Jul 13 '16 at 12:12
• @jose - Nice catch, I just stole that and put it in my answer below – Jason B. Jul 13 '16 at 13:37

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],
GeoGraphics[{{Gray, GeoStyling[{"GeoImage", sdh}],
ad["Polygon"]}, {Red, PointSize[.01], Point[GeoPosition /@ dat]}},
GeoBackground -> None]


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[
SmoothDensityHistogram[
Reverse[dat, 2],
PlotRange -> Full,
Frame -> False,
RegionFunction -> Function[{x, y},
{x, y} ∈ utah],
AspectRatio -> Automatic,
Epilog -> {Red, PointSize[.01], Point@*Reverse /@ dat}]


## 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[

• 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) – Edmund Jul 12 '16 at 19:26