Optimizing spatial query in R

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Optimizing spatial query in R

Ervan Rutishauser-3
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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Optimizing spatial query in R

Chris Reudenbach
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-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)
>

--
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
Reply | Threaded
Open this post in threaded view
|

Re: Optimizing spatial query in R

Ervan Rutishauser-3
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
>> <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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Optimizing spatial query in R

Roman Luštrik
In reply to this post by Chris Reudenbach
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 <
[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
>> <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
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



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

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Optimizing spatial query in R

Tim Keitt-3
In reply to this post by Ervan Rutishauser-3
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 <[hidden email]
> 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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo