Dear Chris, Thanks for your prompt reply and pointing out the "sf" package . I wasn't aware of it, but I am glad to see it has been released only 2 weeks ago ;) I'll have a look right now.
Tschüss Ervan On 20 February 2017 at 10:56, Chris Reudenbach <reudenb...@uni-marburg.de> wrote: > Ervan, > > Just on the run, especially if you deal with real data and some 100K > -1000K of vectors I think the fastest way to do so is to engage GDAL > GRASS7/SAGA and their ability to deal highly efficient with this kind of > topological and geometrical queries. > > A very good pure R alternative is to use the sf package that is ways > faster and more satble than sp. It also provides this type of GIS > functionality like st_overlaps() etc. > > cheers Chris > > > > On 20.02.2017 09:45, Ervan Rutishauser 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-spat >> ial-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,re >> place=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) >> >> > -- > Dr Christoph Reudenbach, Philipps-University of Marburg, Faculty of > Geography, GIS and Environmental Modeling, Deutschhausstr. 10, D-35032 > Marburg, fon: ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: > gis-ma.org, giswerk.org, moc.environmentalinformatics-marburg.de > > -- _____________________________ 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