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

[Robert J. Hijmans]

Sébastien,

You can fix this by rounding the coordinates of the polygons first:

library(akima)
library(maptools)
gpclibPermit()
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
>