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

[Ervan Rutishauser]

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 <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
>
>


-- 
_____________________________
Ervan Rutishauser, PhD
​ <http://carboforexpert.ch/>STRI post-doctoral fellow
CarboForExpert​
​ <http://carboforexpert.ch/>
​Skype: ervan-CH

	[[alternative HTML version deleted]]