[R-sig-Geo] Problems with gDelaunayTriangulation and neighbors.

[Roger Bivand]

On Wed, 18 Dec 2013, Christian Jansson wrote:

> Hi,
>
> I'm very happy that rgeos was updated to also do delaunay triangulation,
> but I have problems to match the results with the points.
>
> I need to to what is in the example-code in gDelaunayTriangulation-help.
> With meuse-data it works, but not with my own data.

Thanks for a reproducible example.

Well, the Meuse bank data are projected, yours are not. So you also run 
into the GEOS scaling issue. Your input longitude for the third point is:

> print(coordinates(tmp.coords[3]), digits=15)
                    lon               lat
[1,] -69.9492587302774 -17.9166152585801

while after having been through the GEOS triangulation operator, you get 
for example:

> print(coordinates(out1[22]), digits=20)
[[1]]
[[1]][[1]]
                           x                      y
[1,] -69.951718409999998016 -17.917070869999999871
[2,] -69.949258729999996831 -17.916615260000000376

where -69.949258729999996831 and -69.9492587302774 do not touch at least 
at the standard scale.

The rgeos function is there simply to expose what GEOS now offers, hence 
the code example on the help page (which you used without referring to it) 
is there to show the hoops you have to jump through to get GEOS to deliver 
triangles or egdes identified by point IDs. As you've seen, deldir does a 
better job at this, but trying alternatives is always sensible.

>
> How can the error in the following code be avoided? So I get out5.
>

Project the geographical coordinates?

Hope this clarifies,

Roger

> -------------
> library(sp)
> library(rgeos)
>
> tmp.coords <- structure(list(lon = c(-69.9601058046182, -69.9307142312949,
> -69.9492587302774, -69.938119947523, -69.9599320928521, -69.958492947506,
> -69.9517184088961, -69.9515127466649, -69.9608290339626,
> -69.9463054201695), lat = c(-17.922645866545, -17.9157740845672,
> -17.9166152585801, -17.9329659455957, -17.9270430312235, -17.933297069203,
> -17.9170708735245, -17.916779044717, -17.9405611407101,
> -17.9098645380403)), .Names = c("lon", "lat"), class = "data.frame",
> row.names = c(NA, -10L))
> coordinates(tmp.coords) <- c("lon", "lat")
> plot(gDelaunayTriangulation(tmp.coords))
> points(tmp.coords)
>
> out <- gDelaunayTriangulation(tmp.coords, onlyEdges=TRUE)
> lns <- slot(slot(out, "lines")[[1]], "Lines")
> out1 <- SpatialLines(lapply(seq(along=lns), function(i)
> Lines(list(lns[[i]]), ID=as.character(i))))
> out2 <- lapply(1:length(out1), function(i) which(gTouches(tmp.coords,
> out1[i], byid=TRUE)))
> out3 <- do.call("rbind", out2)
> o <- order(out3[,1], out3[,2])
> # The previous line generates "Error in out3[, 1] : subscript out of bounds"
>
> out4 <- out3[o,]
> out5 <- data.frame(from=out4[,1], to=out4[,2], weight=1)
> -------------
>
> Thanks.
>
> Best Regards,
> Chris
>
> 	[[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
>

-- 
Roger Bivand
Department of Economics, 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