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

2017-02-20 Thread Tim Keitt 
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  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 5 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
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

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

2017-02-20 Thread 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 
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 5 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
​

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

2017-02-20 Thread 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 5 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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

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

2017-02-20 Thread Ervan Rutishauser 
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 5 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
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo