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