[R-sig-Geo] Problem with gIntersections {rgeos}
Colin Rundel
rundel at gmail.com
Fri Nov 1 17:05:09 CET 2013
This issue is occurring because your poly2 object contains a canopy polygon that is also in the poly object, as such when calculating the intersection the resulting geometry contains that canopy polygon as well as the intersection polygon which is not a valid geometry hence the error. You can either remove the offending canopy from poly, or more likely if you are attempting to work with a single set of canopy polygons you can do something like the following to obtain all the canopy intersections:
inter = gIntersects(poly,byid=TRUE)
inter[!lower.tri(inter)] = FALSE
ids = which(inter,arr.ind=TRUE)
results = list()
for(i in 1:nrow(ids))
{
results[[i]] = gIntersection(poly[ids[i,1],], poly[ids[i,2],])
}
plot(poly)
for(i in 1:length(results))
{
plot(results[[i]],col='red', add=TRUE)
}
-Colin
On Nov 1, 2013, at 9:23 AM, Ervan Rutishauser <er.rutishauser at gmail.com> wrote:
> Dear R users,
>
> I am having hard time while computing the intersections (overlaps) among
> tree crowns in a forest stand. I am using gIntersection() from 'rgeos'
> package, but I get the following error message:
> "Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id,
> "rgeos_intersection") :
> TopologyException: side location conflict at 55.270200771297759
> 44.641891473965337"
>
> It seems that the function does not allow 2 polygons (within the same
> layer) to be connected (touching each other). Is it a known problem? Any
> way to overcome this issue?
> Below is an example. Many thanks for your suggestion and help.
> Best regards,
> Ervan
>
> _____________________________________
>
> library(spatstat)
> library(rgeos)
> library(sp)
> x<- c(37.5, 42, 51, 63, 89.5, 110, 123, 142.5, 143.5, 141.5) # x
> coordinates
> y <- c(202, 27, 53, 39.5, 59, 25.5, 73.5, 98, 38, 5) # y coordinates
> crown <- c(11.8, 9.4, 9.4, 9.3, 9.1, 11.2, 9.2, 9.9, 11.1, 10.1) # the
> crown radius of the trees
> test <- data.frame(x=x,y=y,crown=crown)
> A <- list()
> for (i in 1: nrow(test)) {
> DISC <-
> disc(test$crown[i],centre=c(test$x[i],test$y[i]),npoly=50)[[4]][[1]][1:2]
> # compute the crown area based on crown radius
> DISC$x[50]<-DISC$x[1] # Polygons need to be closed -> closed the circle by
> replacing last values by initial one
> DISC$y[50]<-DISC$y[1]
> A[[i]] <- Polygons(list(Polygon(DISC)),i)
> }
> dat <- data.frame(ID=1:nrow(test),crown.radius=test$crown)
> poly <- SpatialPolygons(A,1:nrow(test))
> Poly <- SpatialPolygonsDataFrame(poly,dat)
> plot(poly)
>
> poly2 <- SpatialPolygons(A[2:3],1:2)
> plot(poly2,border=2,add=T)
> inter <- gIntersection(poly,poly2)
>
> points(55.27, 44.64,bg=3,pch=21)
>
> --
> _____________________________
>
> Ervan Rutishauser, PhD
> carboforexpert.ch <http://carboforexpert.ch/index.html>
> LinkedIn <http://www.linkedin.com/pub/ervan-rutishauser/36/463/a>
> Skype: ervan-CH
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
More information about the R-sig-Geo
mailing list