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

[Roman Luštrik]

The biggest bottleneck was creating a SpatialPolygons object and
calculating gDistance. After computing the SP object outside the for loop,
the run time came down dramatically (see SO post).

Cheers,
Roman

On Mon, Feb 20, 2017 at 10:56 AM, Chris Reudenbach <
reudenbach at 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 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)
>>
>>
> --
> 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
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



-- 
In God we trust, all others bring data.

	[[alternative HTML version deleted]]