[R] ggmap + geom_raster

[Ashim Kapoor]

Dear all,

I want to :

1. Estimate a weighted 2D kernel.
2. Paint a heatmap on a ggmap.

Here is some reproducible data / code (I got it from the internet) :

s_rit <- structure(list(score = c(45, 60, 38, 98, 98, 53, 90, 42, 96,
45, 89, 18, 66, 2, 45, 98, 6, 83, 63, 86, 63, 81, 70, 8, 78,
15, 7, 86, 15, 63, 55, 13, 83, 76, 78, 70, 64, 88, 61, 78, 4,
7, 1, 70, 88, 58, 70, 58, 11, 45, 28, 42, 45, 73, 85, 86, 25,
17, 53, 95, 49, 80, 70, 35, 94, 61, 39, 76, 28, 1, 18, 93, 73,
67, 56, 38, 45, 66, 18, 76, 91, 76, 52, 60, 2, 38, 73, 95, 1,
76, 6, 25, 76, 81, 35, 49, 85, 55, 66, 90), lat = c(28.040961,
27.321633, 27.457342, 26.541129, 27.889476, 26.365284, 28.555024,
26.541129, 26.272728, 28.279994, 27.889476, 28.279994, 26.6674,
26.272728, 25.776045, 26.541129, 30.247658, 26.365284, 25.450123,
27.889476, 26.541129, 27.264513, 26.718652, 28.044369, 28.251435,
27.264513, 26.272728, 26.272728, 28.040961, 30.312239, 27.889476,
26.541129, 26.6674, 27.321633, 26.365284, 28.279994, 26.718652,
30.23286, 28.040961, 30.193704, 30.312239, 28.044369, 27.457342,
25.450123, 30.23286, 30.312239, 30.193704, 28.279994, 30.247658,
26.541129, 26.365284, 28.279994, 27.321633, 25.776045, 26.272728,
30.23286, 30.312239, 26.718652, 26.541129, 25.450123, 28.251435,
28.185751, 25.450123, 28.040961, 27.321633, 28.279994, 27.321633,
27.321633, 27.321633, 28.279994, 26.718652, 28.362308, 27.264513,
26.365284, 28.279994, 30.23286, 25.450123, 28.362308, 25.450123,
25.776045, 30.193704, 28.251435, 27.457342, 27.321633, 28.185751,
27.457342, 27.889476, 26.541129, 26.541129, 30.23286, 30.312239,
26.718652, 25.450123, 26.139258, 28.040961, 30.23286, 26.718652,
28.185751, 28.044369, 28.555024), lon = c(-82.5498, -80.376729,
-82.525985, -81.843986, -82.317701, -81.796389, -81.276464, -81.843986,
-80.207508, -81.331178, -82.317701, -81.331178, -80.072089, -80.207508,
-80.199437, -81.843986, -81.808664, -81.796389, -80.433557, -82.317701,
-81.843986, -80.432125, -80.091078, -82.394639, -81.490407, -80.432125,
-80.207508, -80.207508, -82.5498, -81.575916, -82.317701, -81.843986,
-80.072089, -80.376729, -81.796389, -81.331178, -80.091078, -81.585975,
-82.5498, -81.579846, -81.575916, -82.394639, -82.525985, -80.433557,
-81.585975, -81.575916, -81.579846, -81.331178, -81.808664, -81.843986,
-81.796389, -81.331178, -80.376729, -80.199437, -80.207508, -81.585975,
-81.575916, -80.091078, -81.843986, -80.433557, -81.490407, -81.289394,
-80.433557, -82.5498, -80.376729, -81.331178, -80.376729, -80.376729,
-80.376729, -81.331178, -80.091078, -81.428494, -80.432125, -81.796389,
-81.331178, -81.585975, -80.433557, -81.428494, -80.433557, -80.199437,
-81.579846, -81.490407, -82.525985, -80.376729, -81.289394, -82.525985,
-82.317701, -81.843986, -81.843986, -81.585975, -81.575916, -80.091078,
-80.433557, -80.238901, -82.5498, -81.585975, -80.091078, -81.289394,
-82.394639, -81.276464)), .Names = c("score", "lat", "lon"), row.names =
c(3205L,
8275L, 4645L, 8962L, 9199L, 340L, 5381L, 8998L, 5476L, 4956L,
9256L, 4940L, 6681L, 5586L, 1046L, 9017L, 1878L, 323L, 4175L,
9236L, 8968L, 6885L, 5874L, 9412L, 6434L, 7168L, 5420L, 5680L,
3202L, 1486L, 9255L, 9009L, 6833L, 8271L, 261L, 5024L, 8028L,
1774L, 3329L, 1824L, 1464L, 9468L, 4643L, 4389L, 1506L, 1441L,
1826L, 4968L, 1952L, 8803L, 339L, 4868L, 8266L, 1334L, 5483L,
1727L, 1389L, 7944L, 8943L, 4416L, 6440L, 526L, 4478L, 3117L,
8308L, 4891L, 8290L, 8299L, 8233L, 4848L, 7922L, 5795L, 6971L,
179L, 4990L, 1776L, 4431L, 5718L, 4268L, 1157L, 1854L, 6433L,
4637L, 8194L, 560L, 4694L, 9274L, 8903L, 8877L, 1586L, 1398L,
5865L, 4209L, 6075L, 3307L, 1634L, 8108L, 514L, 9453L, 5210L), class =
"data.frame")


library(ggmap)
library(RColorBrewer)
MyMap <- get_map(location= "Orlando, FL",
                 source="google",
                 maptype="roadmap", crop=FALSE, zoom=7)
YlOrBr <- c("#FFFFD4", "#FED98E", "#FE9929", "#D95F0E", "#993404")
ggmap(MyMap) + stat_density_2d(data=s_rit, aes(x=lon, y=lat,
fill=..level.., alpha=..level..),
                              geom="polygon", size=0.01, bins=16) +
  scale_fill_gradient(low="red", high="green") +
  scale_alpha(range = c(0,0.3), guide=FALSE)

The above computes 2d density but it does not take the score (the weight
into account). To do a weighted KDE I do :

library(ks)

mydensity <- kde(x = s_rit[,c(2,3)],w = s_rit$score)

This computes a weighted KDE,but the data can't be passed to geom_raster
because the x,y's are 151 each while the estimated density is a matrix 151
x 151. However this can be accepted by image as shown below.


library(reshape2)
library(ggplot2)

x <- mydensity$eval.points[[1]]
y <- mydensity$eval.points[[2]]
z <- mydensity$estimate

image(x,y,z)

How can I convert this data so that I can do ggmap + geom_raster?

Also can I define eval.points in kde in a way that is relevant to the
underlying ggmap? How can I tweak the eval.points in the context of a map ?

Best Regards,
Ashim

	[[alternative HTML version deleted]]