Spatial query optimization in R-crown-distributed biomass in the forest

Question

I am trying to calculate the area of ​​cross trees overlapping with square grid cells in a forest. Further reproducible example:

# A. Define objects
    require(sp)
    require(raster)
    require(rgdal)
    require(rgeos)
    require(dismo)
    radius=25   # max search radius around 10 x 10 m cells
    res <- vector() # where to store results

    # Create a fake set of trees with x,y coordinates and trunk diameter (=dbh)
    set.seed(0) 
    survey <- data.frame(x=sample(99,1000,replace=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))  
    coordinates(survey) <- ~x+y

    # Define 10 x 10 subplots
    grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10))) 
    survey$subplot <- over(survey,grid10)   

# B. Now find fraction of tree crown overlapping each subplot
    for (i in 1:100) {
        # Extract centroïd of each the ith cell
        centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]
        corner <- data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))


        # Find trees in a max radius (define above)
        tem <- survey[which((centro$x-survey$x)^2+(centro$y-survey$y)^2<=radius^2),]


        # Define tree crown based on tree diameter
        tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in meter

        # Compute the distance from each tree to cell borders
        pDist <- vector()
        for (k in 1:nrow(tem))  {
            pDist[k] <- gDistance(tem[k,],SpatialPolygons(list(Polygons(list(Polygon(corner)),1))))
        }

        # Keeps only trees whose crown is lower than the above distance (=overlap)
        overlap.trees <- tem[which(pDist<=tem$crownr),]
        overlap.trees$crowna <-overlap.trees$crownr^2*pi  # compute crown area

        # Creat polygons from overlapping crowns
        c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr, lonlat=F, dissolve=F)
        crown <- polygons(c1)
        Crown <-    SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=overlap.trees$dbh,crown.area=overlap.trees$crowna))

        # Create a fine grid points to retrieve the fraction of overlapping crowns
        max.dist <- ceiling(sqrt(which.max((centro$x - overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance to narrow search

        finegrid <- as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))
        coordinates(finegrid) <- ~ x+y
        A <- extract(Crown,finegrid)
        Crown@data$ID <- seq(1,length(crown),1)
        B <- as.data.frame(table(A$poly.ID))
        if (nrow(B)>0) {
        B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)
        B$overlap <- B$Freq/B$crown.area
        B$overlap[B$overlap>1] <- 1
        res[i] <- sum(B$overlap) } else {
        res[i] <- 0 }
    }

# C. Check the result
    res  # sum of crown fraction overlapping each cell (works fine)

This algo takes about 3 minutes to work on 100 cells. I have a large dataset with 35,000 cells, so 150 * 7 = 1050 minutes or 17.5 hours. Any hint at fixing and / or optimizing this algorithm

Answer 1

profvis, , , . , , .

pDist <- vector()
for (k in 1:nrow(tem))  {
    pDist[k] <- gDistance(tem[k,],SpatialPolygons(list(Polygons(list(Polygon(corner)),1))))
}

pDist <- rep(NA, nrow(tem))
my.poly <- SpatialPolygons(list(Polygons(list(Polygon(corner)),1)))
for (k in 1:nrow(tem))  {
  pDist[k] <- gDistance(tem[k,], my.poly)
}

SpatialPolygons . , ( ).

, .

# load only necessary package for code until parSapplyLB
# LB is load-balancing, which means it will distribute task to cores
# which are idle. This is great if jobs take an uneven amount of time
# to run.

library(parallel)
library(sp)

system.time({

  # prepare the cluster, default is PSOCK on windows but can be FORK form *nix
  cl <- makeCluster(4)
  # worker is just a new instance of fresh vanilla R so you need to load the 
  # necessary libraries to all the workers
  clusterEvalQ(cl = cl, library(sp))
  clusterEvalQ(cl = cl, library(raster))
  clusterEvalQ(cl = cl, library(rgdal))
  clusterEvalQ(cl = cl, library(rgeos))
  clusterEvalQ(cl = cl, library(dismo))

  radius <- 25   # max search radius around 10 x 10 m cells
  # res <- rep(NA, 100) # where to store results

  # Create a fake set of trees with x,y coordinates and trunk diameter (=dbh)
  set.seed(0) 
  survey <- data.frame(x=sample(99,1000,replace=T),y=sample(99,1000,replace=T),dbh=sample(100,1000,replace=T))  
  coordinates(survey) <- ~x+y

  # Define 10 x 10 subplots
  grid10 <- SpatialGrid(GridTopology(c(5,5),c(10,10),c(10,10))) 
  survey$subplot <- over(survey,grid10)   

  # Export needed variables to workers
  clusterExport(cl = cl, varlist = c("survey", "radius"))

  # this function is your former for() loop, increase X = 1:100 to suit your needs

  res <- parSapplyLB(cl = cl, X = 1:100, FUN = function(i, survey) {
  # B. Now find fraction of tree crown overlapping each subplot
    # Extract centroïd of each the ith cell
    centro <- expand.grid(x=seq(5,95,10),y=seq(5,95,10))[i,]
    corner <- data.frame(x=c(centro$x-5,centro$x+5,centro$x+5,centro$x-5),y=c(centro$y-5,centro$y-5,centro$y+5,centro$y+5))

    # Find trees in a max radius (define above)
    tem <- survey[which((centro$x-survey$x)^2+(centro$y-survey$y)^2<=radius^2),]

    # Define tree crown based on tree diameter
    tem$crownr <- exp(-.438+.658*log(tem$dbh/10)) # crown radius in meter

    # Compute the distance from each tree to cell borders
    pDist <- vector()
    my.poly <- SpatialPolygons(list(Polygons(list(Polygon(corner)),1)))
    for (k in 1:nrow(tem))  {
      pDist[k] <- gDistance(tem[k,], my.poly)
    }

    # Keeps only trees whose crown is lower than the above distance (=overlap)
    overlap.trees <- tem[which(pDist<=tem$crownr),]
    overlap.trees$crowna <-overlap.trees$crownr^2*pi  # compute crown area

    # Creat polygons from overlapping crowns
    c1 <- circles(coordinates(overlap.trees),overlap.trees$crownr, lonlat=F, dissolve=F)
    crown <- polygons(c1)
    Crown <-    SpatialPolygonsDataFrame(polygons(c1),data=data.frame(dbh=overlap.trees$dbh,crown.area=overlap.trees$crowna))

    # Create a fine grid points to retrieve the fraction of overlapping crowns
    max.dist <- ceiling(sqrt(which.max((centro$x - overlap.trees$x)^2 + (centro$y - overlap.trees$y)^2))) # max distance to narrow search

    finegrid <- as.data.frame(expand.grid(x=seq(centro$x-max.dist,centro$x+max.dist,1),y=seq(centro$y-max.dist,centro$y+max.dist,1)))
    coordinates(finegrid) <- ~ x+y
    A <- extract(Crown,finegrid)
    Crown@data$ID <- seq(1,length(crown),1)
    B <- as.data.frame(table(A$poly.ID))
    if (nrow(B)>0) {
      B <- merge(B,Crown@data,by.x="Var1",by.y="ID",all.x=T)
      B$overlap <- B$Freq/B$crown.area
      B$overlap[B$overlap>1] <- 1
      res <- sum(B$overlap) } else {
        res <- 0 }
  }, survey = survey)
stopCluster(cl = cl)
})

Answer 2

, , , ( . -). " " 1x1m. 6 . !

# Create a fake 1-ha forest stand:
 trees <- data.frame(x=sample(99.5,1000,replace=T),y=sample(99.5,1000,replace=T),dbh=sample(100,1000,replace=T))  


# Create a 1x1m cell matrix where to  store the result
cdagb=matrix(0,nrow=100,ncol=100)


#Calculate the crownradius for every stem (fake proportion)
trees$crownradius = 2*trees$dbh^0.5

#Calculate the index of the 1x1 m quadrat in which the tree stem falls 
trees$quadx=ceiling(trees$x)
trees$quady=ceiling(trees$y)

# Run the algo stem-by-stem
for (i in 1:nrow(trees)) {
    # xdisp and ydisp are the integer cell position differences in x and y that should be checked to see if the crown of the focal tree overlaps
    xdisp=seq(ceiling(trees$quadx[i]-trees$crownradius[i]),floor((trees$quadx[i]+trees$crownradius[i])),1)
    xdisp[xdisp>=1000] <- 1000 +(1000 - xdisp[xdisp>=1000])  # mirror values on edges onto adjacent cells
    xdisp[xdisp<1] <-  -xdisp[xdisp<1] + 1 # avoid XY to be 0
    ydisp=seq(ceiling(trees$quady[i]-trees$crownradius[i]),floor((trees$quady[i]+trees$crownradius[i])),1)
    ydisp[ydisp>=500] <- 500 +(500 - ydisp[ydisp>=500])
    ydisp[ydisp<1] <-  -ydisp[ydisp<1] + 1

    # Calculate the square of the x and y distances from the focal tree to the center of each of these cells
    xdistsqr=(xdisp-trees$quadx[i])^2
    ydistsqr=(ydisp-trees$quady[i])^2
    nx=length(xdisp)
    ny=length(ydisp)

    # Calculate the distance from the center of each cell in the neighborhood to the focal tree               
    distmatrix=matrix(sqrt(rep(xdistsqr,each=ny)+rep(ydistsqr,nx)),nrow=nx,ncol=ny) 

    # includes only trees that overlap the grid cells
    incmatrix=ifelse(distmatrix<trees$crownradius[i],1,0)
    ncells=sum(incmatrix)
    agbpercell=trees$agb[i]/ncells   # divide the biomass by cell
    addagbmatrix=incmatrix*agbpercell  # relloacte biomass by cell

    # add the biomass divided in square meter to each grid point
    cdagb[xdisp,ydisp] = cdagb[xdisp,ydisp] + addagbmatrix
}