Hello,
I am playing around with the function below, which (using the dataset appended further below) produces the attached plot.
The functions requires the
To get the plot use:
Now, as you can se from the chart, 1000 cumulative curves are plotted (in gray), creating a sort of envelope for the black cumulative curve based on the observed data (distance of each point to its nearest neighbor). Plotting 1000 curves and stacking them on top of each other is not so computationallly sensible. I was wondering how I can 'peel out' the envelope in order to just plot the inner 95%, i.e. dropping the top and the bottom 5% of the gray curves.
Data:
I am playing around with the function below, which (using the dataset appended further below) produces the attached plot.
The functions requires the
rgeos
and spatstat
packages. I put some comments to the code to explain what it is actually going on.To get the plot use:
advNNa(springs)
Now, as you can se from the chart, 1000 cumulative curves are plotted (in gray), creating a sort of envelope for the black cumulative curve based on the observed data (distance of each point to its nearest neighbor). Plotting 1000 curves and stacking them on top of each other is not so computationallly sensible. I was wondering how I can 'peel out' the envelope in order to just plot the inner 95%, i.e. dropping the top and the bottom 5% of the gray curves.
Clike:
advNNa < function (feature, studyplot=NULL, buffer=0, B=1000, order=1) {
if(is.null(studyplot)==TRUE){
ch < rgeos::gConvexHull(feature)
region < rgeos::gBuffer(ch, width=buffer)
} else {
region < studyplot
}
dst < spatstat::nndist(coordinates(feature), k=order) #for each point in the input feature dataset, calculate the distance to its nearest neighbor
dst.ecdf < ecdf(dst) #calculate the ECDF of the observed NN distances
dist.rnd.mtrx < matrix(nrow=length(feature), ncol=B) #create a matrix to store the distance of each random point to its nearest neighbor; each column correspond to a random set of points
pb < txtProgressBar(min = 0, max = B, style = 3) # set the progress bar to be used later on within the loop
for (i in 1:B){
rnd < sp::spsample(region, n=length(feature), type='random') #draw a random sample of points within the study region
dist.rnd.mtrx[,i] < spatstat::nndist(coordinates(rnd), k=order) #calculate the NN distances of the random points and store them in the matrix (columnwise)
setTxtProgressBar(pb, i)
}
rnd.ecdf < ecdf(dist.rnd.mtrx[,1]) #calculate the ECDF for the first random dataset
plot(rnd.ecdf, verticals=TRUE, do.points=FALSE, #plot the ECDF of the first random dataset
col="#A9A9A988", xlab="Nearest Neighbor distance (d)",
ylab="G (d)",
main="Refined Nearest Neighbor analysis (G function)",
cex.main=0.95,
xlim=c(min(min(dist.rnd.mtrx), min(dst)), max(max(dist.rnd.mtrx), max(dst))))
for (i in 2:B){
rnd.ecdf < ecdf(dist.rnd.mtrx[,i]) #calculate the ECDF of the remaining random sets
plot(rnd.ecdf, #plot the above and add the plot to the preceding one
verticals=TRUE,
do.points=FALSE,
add=TRUE, col="#A9A9A988",
xlim=c(min(min(dist.rnd.mtrx), min(dst)), max(max(dist.rnd.mtrx), max(dst))))
}
plot(dst.ecdf, #plot the ECDF of the input dataset and add it to the preceding plots
verticals=TRUE,
do.points=FALSE,
add=TRUE,
xlim=c(min(min(dist.rnd.mtrx), min(dst)), max(max(dist.rnd.mtrx), max(dst))))
}
Clike:
springs < new("SpatialPointsDataFrame"
, data = structure(list(OBJECTID_1 = 1:49, OBJECTID = c(5L, 6L, 9L, 11L,
12L, 13L, 14L, 15L, 17L, 18L, 19L, 20L, 22L, 23L, 24L, 25L, 27L,
29L, 30L, 32L, 33L, 35L, 36L, 37L, 41L, 49L, 50L, 52L, 55L, 56L,
57L, 59L, 61L, 62L, 63L, 65L, 66L, 68L, 70L, 71L, 72L, 73L, 75L,
79L, 82L, 83L, 84L, 85L, 87L), location = c(NA, NA, "Marsa",
"Ras ilKnejjes", "Ghallis / Bahar icCaghaq", "Wied irRum",
"Gnien isSultan", "Mtahleb", "Korradino", "Marsa", "Mellieha",
"Ghajn Sfurija / Rdum ilPellegrin", "Mtahleb", "Gebel Ciantar",
"Selmun", "Mdina / Mtarfa", "Mizieb irRih, south of Mellieha",
"TasSanti", "Kalkara, near Mistra Village", "Mtahleb", "Girgenti",
"San Gorg talGhadir, B'buga", "West of Rabat", "West of Mdina",
"Marsaxlokk", "Gebel Ciantar", "West of Rabat", "Zebbiegh / Torri Falka",
"Bingemma", "Marsa", "TalVecca, St. Paul's Bay", "Near Burmarrad",
"Wardija (?)", "Fiddien", "Burmarrad", "Dellimara / Marnisi",
"Mgarr", "West of Pwales", "Mgarr", "Mtahleb / Wied irRum",
"Gnien isSultan, Rabat", "Ghajn Tuffieha", "Mellieha (TalFkieren)",
"Cliff end of Ghajn Zejtuna", "Ghajn Zejtuna", "East of Mellieha",
"West of Mellieha", "Mgarr", "Marsaskala"), comment = c("no correspondence in Grima's Appendix 3",
"no correspondence in Grima's Appendix 3", NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA), toponym = c("Ghajn Kaijet",
"Ghajn Clieb", "Ghajn, tal", "Ghajn Bierda", "Ghajn ta' Bir ilBahar",
"Ghajn Btejtes", "Ghajn Cirani", "Ghajn Corr", "Ghajn Dwieli",
"Ghajn Filep/Qortin/Ghajn\nFormag", "Ghajn talFkieren", "Ghajn Gifra",
"Ghajn Ghorab", "Ghajn Ghulem Alla", "Ghajn Hadid", "Ghajn Hamiem",
"Ghajn Hommed", "Ghajn talKalkara", "Ghajn talKalkara", "Ghajn ilKbira",
"Ghajn ilKbira", "Ghajn Kittien", "Ghajn ilKlieb", "Ghajn Kullija (Qollija)",
"Ghajn ta' Marsaxlokk, tal", "Ghajn Qadi", "Ghajn Qajjied",
"Ghajn ilQasab", "Ghajn Qattus, ta\u0092", "Ghajn Rabib", "Ghajn Razun (Rasul)",
"Ghajn Rihana", "Ghajn Saliba, ta\u0092", "Ghajn ta\u0092 San Pawl",
"Ghajn Selmet", "Ghajn Sender", "Ghajn Sfurija", "Ghajn Stas",
"Ghajn Targa", "Ghajn Tejtes", "Ghajn Tewzien", "Ghajn Tuffieha",
"Ghajn Tuta", "Ghajn Xajxa", "Ghajn Xorok, ta\u0092", "Ghajn Zejtuna",
"Ghajn Znuber, ta\u0092", "TalGhajn", "Wied ilGhajn"), Island = c("Malta",
"Malta", "Malta", "Malta", "Malta", "Malta", "Malta", "Malta",
"Malta", "Malta", "Malta", "Malta", "Malta", "Malta", "Malta",
"Malta", "Malta", "Malta", "Malta", "Malta", "Malta", "Malta",
"Malta", "Malta", "Malta", "Malta", "Malta", "Malta", "Malta",
"Malta", "Malta", "Malta", "Malta", "Malta", "Malta", "Malta",
"Malta", "Malta", "Malta", "Malta", "Malta", "Malta", "Malta",
"Malta", "Malta", "Malta", "Malta", "Malta", "Malta")), .Names = c("OBJECTID_1",
"OBJECTID", "location", "comment", "toponym", "Island"), row.names = c(NA,
49L), class = "data.frame")
, coords.nrs = numeric(0)
, coords = structure(c(444679.6213, 444489.1209, 453000, 440500, 450000,
442000, 445000, 441500, 455898.8655, 454000, 444000, 441000,
442000, 447000, 444689.442, 446000, 443000, 442000, 444500, 441500,
446600, 458000, 444500, 446000, 459000, 447500, 444800, 445000,
444569.3845, 453500, 445000, 447578.0522, 444000, 444000, 446500,
458000, 442500, 443932.5311, 442000, 441800, 444500, 441481.4156,
444000, 443500, 443000, 443131.7628, 440725.362, 443142.0498,
460307.5152, 3971824.5466, 3971673.0724, 3970000, 3973500, 3977500,
3970500, 3972000, 3971000, 3970643.9336, 3971000, 3979000, 3975000,
3970000, 3966500, 3980527.5167, 3972000, 3979000, 3974000, 3979100,
3970500, 3968200, 3966000, 3971600, 3971800, 3967000, 3967000,
3971700, 3975000, 3974490.602, 3971000, 3978200, 3976186.4681,
3976500, 3972000, 3977500, 3966500, 3975500, 3977485.4479, 3975500,
3970800, 3972000, 3976353.484, 3979000, 3980500, 3980000, 3980715.9004,
3978540.4813, 3975342.6171, 3969478.1718), .Dim = c(49L, 2L), .Dimnames = list(
NULL, c("coords.x1", "coords.x2")))
, bbox = structure(c(440500, 3966000, 460307.5152, 3980715.9004), .Dim = c(2L,
2L), .Dimnames = list(c("coords.x1", "coords.x2"), c("min", "max"
)))
, proj4string = new("CRS"
, projargs = "+proj=utm +zone=33 +ellps=intl +units=m +no_defs"
)
)
Attachments

48.2 KB Views: 12
Last edited: