[R-sig-Geo] gIntersection error in rgeos package

[Roger Bivand]

On Mon, 24 Sep 2012, maogui.hu wrote:

> Hi all, I found a problem when using the "gIntersection" in rgeos for 
> two lines. The expected result is a long line shared by two lines. 
> However, the result I obtained is a short line (a part of the afore 
> mentioned long line). I have tried to change the setScale() value for 
> different operation precision. But the result is almost the same. I 
> think may be it is a bug of the function. Attached are the data, r code 
> and results. Can anybody help me to solve it? thanks very much.

Hi,

I think that this is a noding problem. sline1 has a long line along its 
northern boundary, which visually overlaps the segment in sline2 on the 
northern edge of the "bay" extending from the west. However, there is no 
node in sline1 matching the node where sline2 turns north-west, then 
north. The C interface that rgeos uses does not seem to offer noding; 
it is deeper in the C++ code.

What you are trying to do is to approximate the results you might get from 
a topological GIS, that is boundary line lengths between contiguous 
polygons. Usually, topological GIS will node line segments before 
cleaning. Curiously, your shapefiles contain LEFT_FID RIGHT_FID 
attributes, so somewhere the data have been in a topological GIS.

I did try:

row.names(sline2) <- "1"
slines12 <- rbind.SpatialLinesDataFrame(sline1, sline2)
spoly12 <- gPolygonize(slines12)
slines12a <- as(spoly12, "SpatialLines")
gI1 <- gIntersection(slines12a[1], slines12a[2])

but the result is unchanged, not surprisingly, as ?gPolygonize says that 
"Linestrings are expected to be fully noded, as such they must not cross 
and must touch only at endpoints."

On reading:

writeOGR(SpatialPolygonsDataFrame(spoly12, data=as(slines12,
  "data.frame")), ".", "spoly12", driver="ESRI Shapefile")

into GRASS, I see:

v.in.ogr dsn=/home/rsb/tmp/bigshape/rgeos/spoly12.shp layer=spoly12 
output=spoly12
Layer: spoly12
Counting polygons for 2 features...
Importing map 2 features...
...
-----------------------------------------------------
Cleaning polygons, result is not guaranteed!
-----------------------------------------------------
...
Building topology for vector map <spoly12_tmp>...
Building areas...
3 areas built
1 isles built
Attaching islands...
Number of nodes: 4
Number of primitives: 12
Number of points: 0
Number of lines: 0
Number of boundaries: 12
Number of centroids: 0
Number of areas: 3
Number of isles: 1
Number of areas without centroid: 3
-----------------------------------------------------
Find centroids for layer: spoly12
-----------------------------------------------------
Write centroids:
1 areas represent more (overlapping) features, because polygons overlap in 
input layer(s). Such areas are linked to more than 1 row in attribute 
table. The number of features for those areas is stored as category in 
layer 2
-----------------------------------------------------
2 input polygons
Total area: 1.84677E+08 (3 areas)
Overlapping area: 183.124 (1 areas)
Area without category: 0 (0 areas)
-----------------------------------------------------
Building topology for vector map <spoly12>...
Registering primitives...
8 primitives registered
177 vertices registered
Building areas...
3 areas built
1 isles built
Attaching islands...
Attaching centroids...
Number of nodes: 6
Number of primitives: 8
Number of points: 0
Number of lines: 0
Number of boundaries: 5
Number of centroids: 3
Number of areas: 3
Number of isles: 1

So rgeos/GEOS are being confused by the input topologies, in which lines 
that appear to coincide actually do not, and in which the line segments do 
not share the same end points. Reading the polygons back into R from 
GRASS, I get the same intersection length. So I'm afraid that the issue is 
in the input data, which contains a tiny overlap between the polygons, 
leading to a mismatch between the lines. The overlap is only 183 square 
metres, so is not visible on plots.

I think that the only way to repair the data is to add extra nodes, and 
snap them using a digitizer in a GIS, or possibly to edit the coordinates 
of sline1 manually to ensure that they match those of that part of sline2 
as closely as possible.

Hope this clarifies,

Roger

>
> require(maptools)
> require(rgeos)
>
> ws <- "E:/testdata/"
> sline1 <- readShapeLines(file.path(ws, "Line01.shp"))
> sline2 <- readShapeLines(file.path(ws, "Line02.shp"))
>
> Intersection <- function(line1, line2)
> {
> gL <- c()
> nLength <- 0
> for(s in 10^c(-2:8))
> {
> setScale(s)
> gL0 <- gIntersection(line1, line2)
> if(gLength(gL0) > nLength)
> {
> gL <- gL0
> nLength <- gLength(gL)
> }
> }
> setScale()
> return(gL)
> }
>
> par(mfrow=c(2,2))
> plot(sline1, col="red")
> plot(sline2, col="blue")
> plot(rbind.SpatialLines(sline1, sline2, makeUniqueIDs=TRUE), col=c("red","blue"))
> rslt <- Intersection(sline1, sline2)
> plot(rslt, col="green")
>
>
>
>
> maogui.hu

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no