[R-sig-Geo] Polygon dissolution or union issue

[Robert J. Hijmans]

What is the value of the interp function in Akima? I have noticed on
this list that people use it, but I do not why. Is the Tps (thin-plate
spline) method in Fields a sufficient (or better) alternative? In a
comparison using these data (sampling only a few points to 'train' the
interpolator, to make it non-trivial), Tps does _much_ better than
interp in Akima. But perhaps that is different with different
settings. or with different data. Or is interp more efficient with
large datasets? Or?

library(akima)
library(maptools)
gpclibPermit()
library(raster)
library(fields)

data(akima)
akimadf <- data.frame(x=akima$x, y=akima$y, z=akima$z)
dimgrid <- 50

tpcor=akcor=list()

for (i in 1:10) {
  s <- sample(50, 40)
  train <- akimadf[-s,]
  test <- akimadf[s,]
  ak.li <- interp(train$x, train$y, train$z, xo=seq(min(akima$x),
max(akima$x), length = dimgrid), yo=seq(min(akima$y), max(akima$y),
length = dimgrid),linear=FALSE, extrap=TRUE, duplicate = "mean")
  r <- raster(ak.li)
  tps <- Tps(train[,1:2], train[,'z'])
  p <- interpolate(r, tps)
  ak <- extract(r, test[,1:2])
  tp <- extract(p, test[,1:2])
  # correlation between predicted and withheld values
  akcor[i] <- cor(ak, test[,3])
  tpcor[i] <- cor(tp, test[,3])
}
mean(unlist(akcor))  # akima:interp
mean(unlist(tpcor))   # fields::Tps


A typical result:
> mean(unlist(akcor))
[1] 0.3952623
> mean(unlist(tpcor))
[1] 0.9307129
>



On Sun, Dec 5, 2010 at 10:59 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> On Sat, 4 Dec 2010, Robert J. Hijmans wrote:
>
>> Sébastien,
>>
>> You can fix this by rounding the coordinates of the polygons first:
>>
>> library(akima)
>> library(maptools)
>> gpclibPermit()
>
> Please note that both akima and gpclib have restricted licenses (excluding
> research for paying clients), so that this should not be suggested as a
> general solution. The equivalent unionSpatialPolygons() in the as yet
> unreleased maptools on R-forge will use rgeos (also on R-forge) instead -
> rgeos should reach CRAN this month. Replacing akima remains an open
> question, any offers?
>
> Roger
>
>
>> library(raster)
>> data(akima)
>> dimgrid <- 50
>> ak.li <- interp(akima$x, akima$y, akima$z, xo=seq(min(akima$x),
>> max(akima$x), length = dimgrid),yo=seq(min(akima$y), max(akima$y),
>> length = dimgrid),linear = TRUE, extrap=FALSE, duplicate = "mean")
>> r <- raster(ak.li)
>> pol <- rasterToPolygons(r, fun=function(x){x>20 & x<25})
>>
>> ### new code
>> digits <- 6
>> for (i in 1:length(pol at polygons)) {
>> # simple case with no multi-polygons
>>        pol at polygons[[i]]@Polygons[[1]]@coords <-
>> round(pol at polygons[[i]]@Polygons[[1]]@coords, digits)
>> }
>> ###
>>
>> union <- unionSpatialPolygons(pol, ID=rep(1, times=length(pol at polygons)))
>> plot(r)
>> plot(union, add=T)
>>
>>
>> Robert
>>
>> On Sat, Dec 4, 2010 at 10:00 PM, Sébastien Durand <v8extra at gmail.com>
>> wrote:
>>>
>>> Hello,
>>>
>>> Could somebody find me the solution to this issue!
>>>
>>> I cannot manage to merge or dissolve all the polygon  generated in the
>>> object "pol"
>>>
>>> Here is my example:
>>>
>>> library(akima)
>>> library(maptools)
>>> library(raster)
>>>
>>> data(akima)
>>> # define de dimension of grid
>>> dimgrid=50
>>>
>>> # Interpolate to regular grid
>>> ak.li=interp(akima$x, akima$y, akima$z, xo=seq(min(akima$x),
>>> max(akima$x), length = dimgrid),yo=seq(min(akima$y), max(akima$y), length =
>>> dimgrid),linear = TRUE, extrap=FALSE, duplicate = "mean")
>>>
>>> # Show interpolation
>>> image(ak.li)
>>> points(akima)
>>> with(akima, text(x, y, formatC(z,dig=2), adj = -0.1))
>>>
>>> # Create a raster with the data
>>> r <- raster(ak.li)
>>> plot(r)
>>> pol <- rasterToPolygons(r, fun=function(x){x>20 & x<25})
>>> plot(pol, add=T, col='red')
>>>
>>> r <- raster(nrow=length(ak.li$y), ncol=length(ak.li$x), crs=NA)
>>> tmp=ak.li$z
>>> tmp=apply(tmp, 1, rev)
>>> r[] <- tmp
>>>
>>> # Show the raster
>>> plot(r)
>>>
>>> # Convert the data found between ]20 and 25[ into polygons
>>> pol <- rasterToPolygons(r, fun=function(x){x>20 & x<25})
>>>
>>> # Show the polygons
>>> plot(pol, add=T, col='red')
>>>
>>> # Attempt to merge similar adjacent into one polygon
>>> union = unionSpatialPolygons(pol, ID=rep(1,
>>> times=length(pol at polygons)))
>>>
>>> # Show
>>> plot(union, add=TRUE, col="red")
>>>
>>> I beleive all polygons should be merged by that point but obviously I
>>> must have done something wrong.  All polygons were not merged into a single
>>> one and I don't know how to fix that!
>>>
>>> Any other solutions
>>>
>>> Help?
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Roger Bivand
> Economic Geography Section, Department of Economics, Norwegian School of
> Economics and Business Administration, 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