[R-sig-Geo] Optimizing spatial query in R

[Tim Keitt]

I generally use postgis for this sort of thing.

THK

http://www.keittlab.org/

On Mon, Feb 20, 2017 at 2:45 AM, Ervan Rutishauser <er.rutishauser at gmail.com
> wrote:

> Dear All,
>
> I am desperately trying to fasten my algorithm to estimate the fraction of
> tree crown that overlap a given 10x10 subplot in a forest plot. I have
> combined a set of spatial functions (gDistance, extract)  & objects
> (SpatialGrid, SpatialPolygons) in a way that is probably not the most
> efficient, as it takes dozen of hours to run for a 50000 subplots (50 ha
> forest plot).
>
> My detailed problem and a reproducible example are posted on Stackoverflow
> <http://stackoverflow.com/questions/42303559/optimizing-spatial-query-in-r>
> and
> append below, if you wanna have a look. Apart of Amazon Web Server, is
> anyone aware of a sever where to execute (and save results back) R codes
> online?
> Thanks for any help.
> Best regards,
>
>
> # 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's 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 at data$ID <- seq(1,length(crown),1)
>         B <- as.data.frame(table(A$poly.ID))
>         B <- merge(B,Crown at 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)
>     }
> # C. Check the result
>     res  # sum of crown fraction overlapping each cell (works fine)
>
> --
> _____________________________
> Ervan Rutishauser, PhD
> ​ <http://carboforexpert.ch/>STRI post-doctoral fellow
> CarboForExpert​
> ​ <http://carboforexpert.ch/>
> ​Skype: ervan-CH
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

	[[alternative HTML version deleted]]