How to align GeoGraphics "GeoImage" with GeoPositions?

Question

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]

I create the density plot.

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

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]

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

Answer 1

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]

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

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]