Функция R для создания дисков вокруг каждой точки в шаблоне с последующим подсчетом количества точек в каждом диске c [пространственный]

Question

Я пытаюсь создать дис c для каждой точки в паттерне; каждый диск c будет иметь одинаковый радиус. Затем для каждого дис c я хочу посчитать количество точек, попадающих в дис c. Каждый шаблон имеет 100-400 баллов. Я написал код для этого, но он довольно медленный. Код ниже. Я не могу предоставить шейп-файл и точки, так как это будет очень сложно, но я мог бы создать некоторые фиктивные данные, если это будет необходимо.

  W <- as.owin(shape) 
  #Converts created .shp file into a "window" 
  #in which everything is plotted and calculated
  SPDF <- SpatialPointsDataFrame(P[,1:2], P) 
  #Converts data frame to spatial points data frame
  SP <- as(SPDF, "SpatialPoints") #Converts SPDF to spatial points
  SP1  <- as.ppp(coordinates(SP), W)

  SP2 <- as.ppp(SP1)

  attr(SP1, "rejects")
  attr(SP2, "rejects")  



  aw <- area.owin(W) #Area, in pixels squared, of leaf window created earlier
  #awm <- aw * (meas)^2 * 100 #Area window in millimeters squared

  # Trichome_Density_Count-----------------------------------------------------------------------------------------------

  TC <- nrow(P) #Counts number of rows in XY data points file,
  #this is number of trichomes from ImageJ

  TD <- TC/awm #Trichome density, trichomes per mm^2




#SPDF2 <- as.SpatialPoints.ppp(SP2)


#kg <- knn.graph(SPDF2, k = 1) 
#Creates the lines connecting each NND pairwise connection
#dfkg <- data.frame(kg) #Converts lines into a data frame
#dfkgl <- dfkg$length

meanlength <- 78

discstest <- discs(SP2, radii = meanlength, 
                   separate = TRUE, mask = FALSE, trim = FALSE,
                   delta = NULL, npoly=NULL) 
#Function creates discs for each trichome
#Using nearest neighbor lengths as radii


#NEED TO ADD CLIPPING

ratiolist <- c()

for (i in 1:length(discstest)) {



  ow2sp <- owin2SP(discstest[[i]])

  leafsp <- owin2SP(W)

  tic("gIntersection")

  intersect <-  rgeos::gIntersection(ow2sp, leafsp)

  Sys.sleep(1)
  toc()


  tic("over")


  res <- as.data.frame(sp::over(SP, intersect, returnList = FALSE))

  Sys.sleep(1)
  toc()

  res[is.na(res)] <- 0

  newowin <- as.owin(intersect)

  circarea <- area.owin(newowin)

  trichactual <- sum(res)

  trichexpect <- (TC / aw) * circarea

  ratio <- trichactual / trichexpect


  ratiolist[[i]] <- ratio


}

Answer 1

Если я вас правильно понимаю, вы хотите провести l oop через каждую точку и проверить, сколько точек попадает в радиус dis c радиуса R с центром в этой точке. Это очень эффективно выполняется в spatstat с помощью функции closepaircounts:

closepaircounts(SP2, r = meanlength)

Это просто возвращает вектор с количеством точек, содержащихся в dis c радиуса r для каждой точки в SP2.

Я только что попробовал это для 100 000 точек, где каждая точка в среднем имела почти 3000 других точек в диске c, и на моем ноутбуке это заняло 8 секунд. Если у вас есть еще много точек или, в частности, если радиус дис c настолько велик, что каждый дис c содержит гораздо больше точек, его вычисление может стать очень медленным.