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

Roger Bivand Roger.Bivand at nhh.no
Sun Dec 5 19:59:09 CET 2010


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


More information about the R-sig-Geo mailing list