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.
C-like:
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 (column-wise)
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))))
}
C-like:
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 il-Knejjes", "Ghallis / Bahar ic-Caghaq", "Wied ir-Rum",
"Gnien is-Sultan", "Mtahleb", "Korradino", "Marsa", "Mellieha",
"Ghajn Sfurija / Rdum il-Pellegrin", "Mtahleb", "Gebel Ciantar",
"Selmun", "Mdina / Mtarfa", "Mizieb ir-Rih, south of Mellieha",
"Tas-Santi", "Kalkara, near Mistra Village", "Mtahleb", "Girgenti",
"San Gorg tal-Ghadir, B'buga", "West of Rabat", "West of Mdina",
"Marsaxlokk", "Gebel Ciantar", "West of Rabat", "Zebbiegh / Torri Falka",
"Bingemma", "Marsa", "Tal-Vecca, St. Paul's Bay", "Near Burmarrad",
"Wardija (?)", "Fiddien", "Burmarrad", "Dellimara / Marnisi",
"Mgarr", "West of Pwales", "Mgarr", "Mtahleb / Wied ir-Rum",
"Gnien is-Sultan, Rabat", "Ghajn Tuffieha", "Mellieha (Tal-Fkieren)",
"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 il-Bahar",
"Ghajn Btejtes", "Ghajn Cirani", "Ghajn Corr", "Ghajn Dwieli",
"Ghajn Filep/Qortin/Ghajn\nFormag", "Ghajn tal-Fkieren", "Ghajn Gifra",
"Ghajn Ghorab", "Ghajn Ghulem Alla", "Ghajn Hadid", "Ghajn Hamiem",
"Ghajn Hommed", "Ghajn tal-Kalkara", "Ghajn tal-Kalkara", "Ghajn il-Kbira",
"Ghajn il-Kbira", "Ghajn Kittien", "Ghajn il-Klieb", "Ghajn Kullija (Qollija)",
"Ghajn ta' Marsaxlokk, tal-", "Ghajn Qadi", "Ghajn Qajjied",
"Ghajn il-Qasab", "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", "Tal-Ghajn", "Wied il-Ghajn"), 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: