IDW basemap.R

# IDW surface with basemap ------------------------------------------------

# carrying on from section 4.5.2 of the book
# https://gdsl-ul.github.io/san/points.html#a-surface-of-housing-prices

# version in the book using geom_raster()
ggplot(idw.hp, aes(x = X, y = Y, fill = var1.pred)) +
  geom_raster() +
  scale_fill_viridis_b() +
  theme_void() +
  geom_sf(alpha=0)

# the same thing but with geom_tile() (identical to the above)
ggplot(idw.hp, aes(x = X, y = Y, fill = var1.pred)) +
  geom_tile() +
  scale_fill_viridis_b() +
  theme_void() +
  geom_sf(alpha=0)

# can further interpolate with geom_raster 
# (although not very well - just blurs the edges really)
ggplot(idw.hp, aes(x = X, y = Y, fill = var1.pred)) +
  geom_raster(interpolate = T) +
  scale_fill_viridis_b() +
  theme_void() +
  geom_sf(alpha=0)

# using stat_summary_2d()
ggplot(idw.hp, aes(x = X, y = Y, z = var1.pred)) +
  stat_summary_2d(bins = 100) + # original: 100 (ie sd.grid resolution)
  scale_fill_viridis_b(alpha = 1) + 
  theme_void() +
  geom_sf(alpha=0)

# strip idw.hp of anything but the predicted house price, lat, and lon values
# otherwse qmplot throws an error
raw.idw.hp <- idw.hp %>% as.data.frame %>% select(var1.pred, X, Y)

# stat_summary_2d() with basemap
qmplot(X, Y, data = raw.idw.hp, geom = "blank") + 
  stat_summary_2d(aes(x = X, y = Y, z = var1.pred), 
                  bins = 100) + # original: 100 (ie sd.grid resolution)
  scale_fill_viridis_b(alpha = 0.3) + # uniform alpha value across the surface 
  theme_void()

# geom_tile() with alpha set by predicted house price
qmplot(X, Y, data = raw.idw.hp, geom = "blank") +
  geom_tile(data = raw.idw.hp,
            aes(x = X, y = Y, fill = var1.pred, alpha = var1.pred)) +
  scale_fill_viridis_b() + # direction = -1 to reverse
  theme_void()

## zooming in to see the main variation close up

# crop IDW surface to desired geographic bounds
crop.idw <- raw.idw.hp[raw.idw.hp$X > -117.28 
                       & raw.idw.hp$X < -117.2 
                       & raw.idw.hp$Y > 32.8 
                       & raw.idw.hp$Y < 32.87 , ]

# as above but with the cropped data
qmplot(X, Y, data = crop.idw, geom = "blank") +
  geom_tile(data = crop.idw,
            aes(x = X, y = Y, fill = var1.pred, alpha = var1.pred)) +
  scale_fill_viridis_b() + # direction = -1 to reverse
  theme_void()

# version where the cropping is done ad hoc by qmplot
# (the map zoom ends up at the wrong level)
qmplot(X, Y, data = raw.idw.hp, geom = "blank",
       xlim = c(-117.28, -117.2), ylim = c(32.8, 32.87)) +
  geom_tile(data = raw.idw.hp,
            aes(x = X, y = Y, fill = var1.pred),
            alpha = 0.3) +
  scale_fill_viridis_b() + # direction = -1 to reverse
  theme_void()