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.rutishau...@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@data$ID <- seq(1,length(crown),1) > B <- as.data.frame(table(A$poly.ID)) > 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) > } > # 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@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo