1

I'm trying

set97 <- gIntersection(setbp,d97)

and getting the error message:

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_intersection") : 
  TopologyException: found non-noded intersection between LINESTRING (533036 -314770, 533036 -314770) and LINESTRING (533036 -314780, 533036 -314770) at 533035.88110651996 -314769.97350772272

(setbp is a sector of a county - for census purposes - and d97 is the deforestation in brazilian Amazon up to 1997.) When I do

gIsValid(d97)
gIsValid(setbp)
gIntersects(setbp,d97)

I get TRUE for all three questions. So, what's wrong? I'll try gBuffer with 0 width (it solves some of these issues), but I'd like to understand what's going on.

EDIT: gBuffer didn't work. @mdsumner suggested me to do the gIntersection row by row, but up to now (after trying 3 different methods), it didn't work either (look at the comments). Stranger yet, the TopologyException is being caused by strange "strips" of geometry not detected by gIsValid, and worse yet, kilometers apart from the intersecting polygon (setbp above, setbpi in the comments). I'll try the function shown here, though it doesn't seem a rounding problem.

Rodrigo
  • 4,706
  • 6
  • 51
  • 94
  • 1
    According to @mdsumner in the SO thread below, "It means there is a line crossing another line and no intermediate coordinate records the intersection".. http://stackoverflow.com/questions/13662448/what-does-the-following-error-mean-topologyexception-found-non-nonded-intersec – geotheory Aug 02 '13 at 14:20
  • 1
    I don't know about that mdsumner guy's thing, but what is the class of d97 and setbp? Tell us class(d97) and class(setbp), if they are Spatial*DataFrame it might be worth doing the intersection row by row. – mdsumner Aug 02 '13 at 14:29
  • 1
    Thanks, geotheory! But in this case, why doesn't gIsValid accuse any problem? – Rodrigo Aug 02 '13 at 14:31
  • mdsumner, class command give the same for both: [1] "SpatialPolygonsDataFrame" attr(,"package") [1] "sp" But I did it exactly like that before with other polygons (hidrography and vegetation) and it worked fine! – Rodrigo Aug 02 '13 at 14:32
  • Any example on how doing this, mdsumner? Thank you! – Rodrigo Aug 02 '13 at 14:34
  • BTW, gBuffer with width=0 or 0.0000001 or 0.1 didn't work. (The help doesn't say, but this width is relative to distance in the proj4string? Since my data is in UTM, the distance is in meters?) – Rodrigo Aug 02 '13 at 14:38
  • 1
    for (i in j to q) for (k in fq2) gIntersection(setbp[i, ], d97[k,]) seriously though, try some things out for yourself – mdsumner Aug 02 '13 at 14:40
  • Thank you. Yes, I'm testing a lot of things by myself. Just asked because maybe you had an example or link at hand... – Rodrigo Aug 02 '13 at 17:46
  • So, inside the loop I tried "set97[length(set97)+1,] <- gIntersection(setbpi,d97[j,])", but it said "object of type 'S4' is not subsettable". Then I tried "rbind(set97,gIntersection(setbpi,d97[j,]))", but it said "invalid class “SpatialPolygons” object: non-unique Polygons ID slot values". Any help is appreciated. – Rodrigo Aug 02 '13 at 18:15
  • Now, following this [link](https://stat.ethz.ch/pipermail/r-sig-geo/2008-June/003725.html), I've used: "set97 <- SpatialPolygons(c(slot(set97,"polygons"),list(gIntersection(setbpi,d97[j,]))))", but it said: "no slot of name "area" for this object of class "SpatialPolygons"". Any ideas, @mdsumner, @geotheory? – Rodrigo Aug 02 '13 at 18:38
  • @mdsumner, I spent all day trying to do that, and achieved not too much. If it's easy too you, please tell me how to do it. Thank you. – Rodrigo Aug 02 '13 at 21:35
  • length(set97)+1 won't work because the length is the length! – mdsumner Aug 02 '13 at 22:10
  • "I tried 3 ways" is pointless unless you tell us exactly how you tried and give reproducible examples. We can't "tell you how to do it" if we don't have clear idea of what you are trying. What do the shapes look like? What do you intend to get from the intersection? Etc. etc. ad infinitum. – mdsumner Aug 02 '13 at 22:21
  • Thanks, @mdsumner. In the main message is written what the shapefiles are. One is small, part of a county. The other is big, pieces of deforestation all over brazilian Amazon, including some polygons inside the part of the county in question. The result I want is just the polygons of deforestation inside that part of a county. Should I send a link to the actual shapefiles? And the three examples are described in the comments, along with the error messages (what's also written in the main message). – Rodrigo Aug 03 '13 at 12:32
  • @mdsumner and others, here you can reproduce the error: http://stackoverflow.com/questions/18083363/gintersection-returns-topologyexception-when-gisvalid-in-r – Rodrigo Aug 06 '13 at 14:47

0 Answers0