A quick demo of calcuating distances between points and then working with that data. Makes it through thinking of distance pairs as a network. Feel free to fork and improve.

# libraries
library(curl)
library(ggplot2)
library(igraph)
library(sp)
library(tidyr)

# load amphitheater data
library(cawd)

#generate distances of all amphitheaters to all other amphitheaters
ramphs.dists <- spDists(ramphs.sp,ramphs.sp,longlat = TRUE)

# this is in the form of a matrix
#print the distance between the Dura and Arles Amphitheaters
ramphs.dists[1,2]
## [1] 3241.864
# but let's use names
colnames(ramphs.dists) <- ramphs$id
rownames(ramphs.dists) <- ramphs$id

# now print distane between Arles Amphitheater and Lyon Amphitheater
ramphs.dists["arlesAmphitheater", "lyonAmphitheater"]
## [1] 232.9044
# we also have useless info that amphitheaters are 0Km from themselves
ramphs.dists["arlesAmphitheater", "arlesAmphitheater"]
## [1] 0
# Would be nice to have a d.f in the form from,to,km
as.data.frame(ramphs.dists) -> ramphs.dists.df

#add column of names, in anticipation of using gather()
cbind(from = ramphs$id,ramphs.dists.df) -> ramphs.dists.df

gather(ramphs.dists.df, to, km, 2:dim(ramphs.dists.df)[2]) -> ramphs.dists.df

# data now in nice columns but...
# we still have the 0Km rows.
head(ramphs.dists.df)
##                         from                      to       km
## 1    duraEuroposAmphitheater duraEuroposAmphitheater    0.000
## 2          arlesAmphitheater duraEuroposAmphitheater 3241.864
## 3           lyonAmphitheater duraEuroposAmphitheater 3255.734
## 4           ludusMagnusArena duraEuroposAmphitheater 2579.209
## 5    romeFlavianAmphitheater duraEuroposAmphitheater 2579.429
## 6 romeAmphitheatrumCastrense duraEuroposAmphitheater 2577.511
#set phasers to kill...
ramphs.dists.df[ramphs.dists.df$km > 0,] -> ramphs.dists.df

# yay, they're gone
head(ramphs.dists.df)
##                         from                      to        km
## 2          arlesAmphitheater duraEuroposAmphitheater 3241.8637
## 3           lyonAmphitheater duraEuroposAmphitheater 3255.7344
## 4           ludusMagnusArena duraEuroposAmphitheater 2579.2087
## 5    romeFlavianAmphitheater duraEuroposAmphitheater 2579.4293
## 6 romeAmphitheatrumCastrense duraEuroposAmphitheater 2577.5114
## 7 eleutheropolisAmphitheater duraEuroposAmphitheater  645.7789
# printing pairs is wordier than using the matrix
ramphs.dists.df$km[ramphs.dists.df$from == 'arlesAmphitheater' & ramphs.dists.df$to == 'lyonAmphitheater']
## [1] 232.9044
# we record both directions, that's probably not a problem but do keep it in mind
ramphs.dists.df$km[ramphs.dists.df$from == 'lyonAmphitheater' & ramphs.dists.df$to == 'arlesAmphitheater']
## [1] 232.9044
#reasonably straigtforward to find closest neighbor (target for improvement)
ramphs.dists.df[ramphs.dists.df$from == "duraEuroposAmphitheater",] -> tmp.df
as.character(tmp.df$to[which.min(tmp.df$km)]) # avoid extra factor output w/as.character()
## [1] "palmyraAmphitheater"
# and note that we can think of this df as edges of a network
graph.data.frame(ramphs.dists.df, directed = TRUE) -> ramphs.dists.nw

# let's see groups of amphitheaters close to each other
subgraph.edges(ramphs.dists.nw, eids = E(ramphs.dists.nw)[km < 99]) -> close.amphitheaters.nw
close.amphitheaters.nw
## IGRAPH DN-- 177 1206 -- 
## + attr: name (v/c), km (e/n)
par(mai = c(0,0,0,0))
plot(close.amphitheaters.nw,
     edge.arrow.size = 0,
     vertex.size = 2,
     vertex.label.cex = .3)

# now let's look for "close" amphitheaters with high network degree.
# those are deeply integrated into dense groupings.

degree(close.amphitheaters.nw) -> tmp.nw

set.vertex.attribute(close.amphitheaters.nw,
                     name = "degree",
                     index = V(close.amphitheaters.nw),
                     value = tmp.nw) -> close.amphitheaters.nw

summary(V(close.amphitheaters.nw)$degree)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00    2.00    6.00   13.63   22.00   46.00
# top quartile
summary(V(close.amphitheaters.nw)$degree)[5]
## 3rd Qu. 
##      22
# print the top quartile
V(close.amphitheaters.nw)[degree >= summary(V(close.amphitheaters.nw)$degree)[5]]
## Vertex sequence:
##  [1] "ludusMagnusArena"                 "romeFlavianAmphitheater"         
##  [3] "romeAmphitheatrumCastrense"       "pompeiiAmphitheater"             
##  [5] "cumaeAmphitheater"                "pozzuoliFlavianAmphitheater"     
##  [7] "pozzuoliEarlyAmphitheater"        "hadrianVillaAmphitheater"        
##  [9] "paestumAmphitheater"              "casinumAmphitheater"             
## [11] "avellaAmphitheater"               "oudnaAmphitheater"               
## [13] "carsulaeAmphitheater"             "ocriculumAmphitheater"           
## [15] "lucusFeroniaeAmphitheater"        "allifaeAmphitheater"             
## [17] "hispellumAmphitheatre"            "albaFucensAmphitheater"          
## [19] "sutriumAmphitheatre"              "tiburAmphitheater"               
## [21] "telesiaAmphitheatre"              "assisiAmphitheatre"              
## [23] "castraAlbanaAmphitheater"         "sanBenedettoDeiMarsiAmphitheater"
## [25] "monteleoneAmphitheater"           "bolsenaAmphitheater"             
## [27] "ascoliPicenoAmphitheater"         "beneventoAmphitheater"           
## [29] "calviRisortaAmphitheater"         "sanVittorinoAmphitheater"        
## [31] "capuaRepublicanAmphitheater"      "capuaImperialAmphitheater"       
## [33] "liternoAmphitheater"              "noleAmphitheater"                
## [35] "noceraSuperioreAmphitheater"      "aquinoAmphitheater"              
## [37] "seressiAmphitheater"              "genzanoPrivateAmphitheater"      
## [39] "hencherElKasbatAmphitheater"      "portusOvalStructure"             
## [41] "statiliusTaurusAmphitheater"      "faleriiNoviAmphitheater"         
## [43] "bevanaAmphitheater"               "minturnaeAmphitheater"           
## [45] "sinuessaAmphitheater"
# next up, exploring where those are (Italy, anyone?)