I succeeded with @eblonde in the discussion (problem with polygon smoothing in R) to smoothing polygons in R.
Code:
inter1= read.table("c:/inter2.csv", header=TRUE)
#add a category (required for later rasterizing/polygonizing)
inter1 <- cbind(inter1,
cat
= rep(1L, nrow(inter1)),stringsAsFactors = FALSE)
#convert to spatial points
coordinates(inter1) = ~long + lat
#gridify your set of points
gridded(inter1) <- TRUE
#convert to raster
r <- raster(inter1)
#convert raster to polygons
sp = rasterToPolygons(r, dissolve = T)
#
# Splining a polygon.
#
# The rows of 'xy' give coordinates of the boundary vertices, in order.
# 'vertices' is the number of spline vertices to create.
# (Not all are used: some are clipped from the ends.)
# 'k' is the number of points to wrap around the ends to obtain
# a smooth periodic spline.
#
# Returns an array of points.
#
spline.poly <- function(xy, vertices, k=3, ...) {
# Assert: xy is an n by 2 matrix with n >= k.
# Wrap k vertices around each end.
n <- dim(xy)[1]
if (k >= 1) {
data <- rbind(xy[(n-k+1):n,], xy, xy[1:k, ])
} else {
data <- xy
}
# Spline the x and y coordinates.
data.spline <- spline(1:(n+2*k), data[,1], n=vertices, ...)
x <- data.spline$x
x1 <- data.spline$y
x2 <- spline(1:(n+2*k), data[,2], n=vertices, ...)$y
# Retain only the middle part.
cbind(x1, x2)[k < x & x <= n+k, ]
}
#addition transformation to distinguish well the set of polygons
polys <- slot(sp@polygons[[1]], "Polygons")
output <- SpatialPolygons(
Srl = lapply(1:length(polys),
function(x){
p <- polys[[x]]
#applying spline.poly function for smoothing polygon edges
px <- slot(polys[[x]], "coords")[,1]
py <- slot(polys[[x]], "coords")[,2]
bz <- spline.poly(slot(polys[[x]], "coords"),40, k=3)
bz <- rbind(bz, bz[1,])
slot(p, "coords") <- bz
# create Polygons object
poly <- Polygons(list(p), ID = x)
return(poly)
}),
proj4string = CRS("+init=epsg:4326")
)
#plot
par(mar = c(0, 0, 0, 0))
Al = readShapeLines("DZA_adm0.shp")
plot(Al, col = "gray")
#plot(sp, border = "gray", lwd = 2) #polygonize result
plot(output, border = "red", add = TRUE) #smoothed polygons
My problem now I need to be represented on a map, in fact I'm looking methodology to clear the area outside my shapefiles.
Is that possible to do in R?
