pls mapping

plsmap.r

#PLS Dot Density Map

library(maps)
library(mapproj)
library(maptools)
library(geosphere)
library(splancs)
library(foreign)

# Map just the points
map("usa", col="gray", bg=NA, lwd=0.4)
points(pls12cai14$LONGITUDE, pls12cai14$LATITUDE, pch=21, cex=0.1, col="#00000020")


# Projected points
par(mar=c(0,0,0,0))
m.proj <- map("usa", proj="albers", param=c(39, 45), col="#000000", fill=FALSE, bg=NA, lwd=0.4)
loc.proj <- mapproject(pls12cai14$LONGITUDE, pls12cai14$LATITUDE)
pts <- as.points(list(x=loc.proj$x, y=loc.proj$y))
points(pts, pch=21, cex=0.1, col="#00000020")

# Color palette
blackWhite <- colorRampPalette(c("white", "black"))

# Hack to use US border as polygon
m.proj <- map("usa", proj="albers", param=c(39, 45), col=NA, fill=FALSE, bg=NA, lwd=0.4)
xlim <- match(NA, m.proj[["x"]]) - 1
ylim <- match(NA, m.proj[["y"]]) - 1
uspoly.proj <- as.matrix(cbind(m.proj[["x"]][1:xlim], m.proj[["y"]][1:ylim]))

# Density
smoo.proj <- kernel2d(pts, uspoly.proj, h0=0.02, nx=5000, ny=5000)

# Draw map
image(smoo.proj, uspoly.proj, add=TRUE, col=blackWhite(100))
map("state", proj="albers", param=c(39, 45), col="#999999", fill=FALSE, bg=NA, lwd=0.2, add=TRUE, resolution=1)