[R-sig-Geo] Creating polygons out of an irregular data set to represent specific classes

Sébastien Durand v8extra at gmail.com
Fri Dec 3 06:03:29 CET 2010


Ok here is my latest development on the topic! But I still have a question! At the end I still get segmented polygons!

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(nrow=length(ak.li$y), ncol=length(ak.li$x), crs=NA)
xmin(r) <- min(ak.li$x, na.rm=TRUE)
xmax(r) <- max(ak.li$x, na.rm=TRUE)
ymin(r) <- min(ak.li$y, na.rm=TRUE)
ymax(r) <- max(ak.li$y, na.rm=TRUE)
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")


All polygons were not merged into a single one and I don't know how to fix that!

Any help?


More information about the R-sig-Geo mailing list