Two ways to make “heatmaps” of the distributioin of Roman Amphitheaters.

library(curl)
library(ggplot2)
library(raster)
library(rgdal)
library(RColorBrewer)
library(rgl)
library(sp)
library(spatstat)
# See http://
library(cawd)

#plot Roman empire to confirm load
par(mai=c(0,0,0,0))
plot(awmc.roman.empire.200.sp)

# but that's not a very nice projection when done in plot, so...
plot(spTransform(awmc.roman.empire.200.sp, CRS("+init=epsg:3857")))

# ggplot2 makes nice heatmaps and its coord_map() handles WGS84 by default

#first convert to form ggplot2 likes
fortify(awmc.roman.empire.200.sp) -> awmc.roman.empire.200.gg

# do a little imputation to create missing data
ramphs$imputed.capacities <- NA
ramphs$imputed.capacities[!(is.na(ramphs$capacity))] <- ramphs[!(is.na(ramphs$capacity)),]$capacity
ramphs$imputed.capacities[is.na(ramphs$capacity)] <- ramphs[is.na(ramphs$capacity),]$ext.major * 250 # 250 seats per exterior meter
ramphs$imputed.capacities[is.na(ramphs$imputed.capacities)] <- 7000

# the stat_density2d() call is what makes the heatmap
ggplot(ramphs, aes(x = longitude, y = latitude))  +
  coord_map() +
  geom_polygon(data=awmc.roman.empire.200.gg,aes(x=long, y=lat,group=group), colour="purple", fill = "purple", alpha = .5) +
  geom_density2d(aes(fill = ..level..),kernel = "cosine",) +
  geom_point( colour = 'black', size = 2) +
  geom_point(aes(size = capacity), colour = 'red') +
  scale_fill_gradientn(colours=rev(brewer.pal(7,"Spectral")))   +
  theme_bw() 

# we can also use spatstat to make the heat map

# yes, I could calculate these bounds
longitudeExtent <- c(-9, 44)
latitudeExtent <- c(31, 56)

# first make a spatstat ppp object with a rectangular owin (that's the limits of the point pattern)
tmp.ppp <- as.ppp(dplyr::select(ramphs, longitude,latitude),owin(longitudeExtent,latitudeExtent))

# expand owin for analysis by calculating the convex hull of the amphitheaters.
# but expand it a little bit so heat map doesn't cut off too much

tmp.hull <- spatstat::convexhull(tmp.ppp)
# while we're at it, plot the hull
par(mai=c(0,0,0,0)) # do I really have to repeat this?
plot(awmc.roman.empire.200.sp)
plot(tmp.hull, add=T)

tmp.owin <- expand.owin(tmp.hull, distance = .8)

# so now we can make ramphs.ppp with a decent owin
ramphs.ppp <- as.ppp(dplyr::select(ramphs, longitude,latitude, imputed.capacities),W = tmp.owin)

# make a grid of the density weighted by the imputed capacity
# but right now, the imputed capacity is way too fancy a term
# i'm just filling in 10K for when I don't have that info
density <- density(ramphs.ppp, weights = ramphs$imputed.capacities)
raster(density) -> r
plot(r)
plot(awmc.roman.empire.200.sp, add = T)

# and in 3d sort of
persp(r,phi = 30,
      theta = -30,
      expand = .25,
      border = NA,
      shade = .6) # yeah, needs a color scheme. i'll return to that...