I succeeded with @eblonde in the discussion (problem with polygon smoothing in R) to smoothing polygons in R.


inter1= read.table("c:/inter2.csv", header=TRUE)
#add a category (required for later rasterizing/polygonizing)
inter1 <- cbind(inter1, 
                = 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),
                 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)
  proj4string = CRS("+init=epsg:4326")

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?

enter image description here


You can do a spatial intersection between your two datasets ("Al", "output") with the following code:

inter = gIntersection(Al, output, byid=T, drop_not_poly=T)
| improve this answer | |

Not the answer you're looking for? Browse other questions tagged or ask your own question.