# Optimizing spatial query in R

5 messages
Open this post in threaded view
|

## Optimizing spatial query in R

 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 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 ​ STRI post-doctoral fellow CarboForExpert​ ​ ​Skype: ervan-CH         [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Open this post in threaded view
|

## Re: Optimizing spatial query in R

 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 > 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) > -- 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 _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Open this post in threaded view
|

## Re: Optimizing spatial query in R

 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 <[hidden email]> 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 >> > 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 ​ STRI post-doctoral fellow CarboForExpert​ ​ ​Skype: ervan-CH         [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo