2

I am attempting to create a disc for each point in a pattern; each disc will have the same radius. Then for each disc, I want to count the number of points falling within the disc. Each pattern has 100-400 points. I have written code to do this, but it is quite slow. The code is below. I cannot provide the shapefile and points as that would be very difficult, but I could create some dummy data if need be. 


  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


}
James Fischer
  • 31
  • 3

1 Answers1

1

If I understand you correctly you want to loop through each point and check how many points fall within a disc of radius R centered in that point. This is done very efficiently in spatstat with the function closepaircounts:

closepaircounts(SP2, r = meanlength)

This simply returns a vector with the number of points contained in the disc of radius r for each point in SP2.

I have just tried this for 100,000 points where each point on average had almost 3000 other points in the disc around it, and it took 8 seconds on my laptop. If you have many more points or in particular if the disc radius is so big that each disc contains many more points it may become very slow to calculate this.

Ege Rubak
  • 4,347
  • 1
  • 10
  • 18
  • I knew there would be a much more simple solution! Thank you! Do you happen to know if the closepaircounts function takes into account the window the points are contained in? I want to calculate density of the points in each circle, but if the circles aren't clipped by the window, I cannot do that accurately. – James Fischer Mar 10 '20 at 00:49
  • I'm not sure I understand you completely and it may be fiddly to answer in a comment. Maybe you can state the problem clearly in a new question? For each point you want to draw a disc of radius R and clip it to the overall window (only relevant for points close to the border) and then measure the area of the potentially clipped disc? – Ege Rubak Mar 10 '20 at 15:21
  • Exactly! I can use discs() then in a loop clip each disc to the window, and measure the area of the clipped edge discs, but this takes a while. – James Fischer Mar 12 '20 at 17:42
  • If you want full details you should really open a new question. Here is a quick version that fits in a comment. Points where disc overlaps border: `i <- which(bdist.points(SP2) – Ege Rubak Mar 12 '20 at 23:56
  • Works perfectly! Thank you for the help! – James Fischer Mar 13 '20 at 03:14