[R-sig-Geo] Clearing the area outside of a polygon defined on a grid.

peleve pitiur at gmail.com
Wed Nov 14 19:45:15 CET 2012


I have a result from Akima or Loesss interpolation using irregular data
defined inside a polygon. I only want the results inside the polygon.
Akima interpolates correctly inside the polygon but leaves a blue residue
outside the polygon (a result of Akima using a convex hull) that I want to
get rid of, and so far, have not been able to do so.  I only want the
information inside the polygon.
I've searched and searched, and read and read, but can't find a way to do
it.
A similar result with Loess, but it uses the whole grid.
Below the code of a highly simplified geometry that demonstrates the issue.
The first figure is what Akima produces. The second figure produced uses a
raster package approach using Akima's output but the data are flipped around
the y-axis.  It nearly seems to do the job, but there still is the outline
and contents of the convex hull that Akima does it's interpolation on.
How do I fix this, please?

library(raster)		# For raster
library(tgp)			# For interp.loess
library(fields)		# For image.plot
library(akima)		# For akima
library(sp)			# For spatial representation


#================================================================
# Definition of the polygon

x <- c(47, 44, 33, 20,19,30,16,4,16,34,45,52,60,60, 47)
y <- c(0,13,20,20,24,30,37,42,45,40,28,20,25,0, 0)
zero <- seq(0, 0, length=length(x))


# #
# # Definition of irregularly spaced -z- values
# # The points of the polygon have been add to force a zero constraint
# #
xx <- c(48,44,33,25,22,30,16,12,16,34,46,52,58,58,32,52,x)
yy <- c(1,14,22,22,22,28,38,42,43,36,21,18,20,4,30,12,y)
zz <- c(5, 6, 6, 4, 7, 5, 4, 3, 5, 6 ,8, 9, 6, 3, 5, 9, zero)


#================================================================
# Definition of the bounding box and related variables

xmin <- 0
xmax <- 60
ymin <- 0
ymax <- 50
dxy <- 0.2
nx <- (xmax-xmin)/dxy
ny <- (ymax-ymin)/dxy
xbox <- c(xmin, xmax, xmax, xmin)
ybox <- c(ymin, ymin, ymax, ymax)



#================================================================
# Akima interpolation/smoothing

x0 <- seq(xmin, xmax, length = nx)
y0 <- seq(ymin, ymax, length = ny)
akima <- interp(xx, yy, zz, xo=x0, yo=y0, duplicate="strip")
image.plot(akima)
polygon(x, y, lwd=2)
points(xx,yy)
contour(akima$x, akima$y, akima$z, add=TRUE, labcex=1.2)
lines(x,y, col="purple", lwd=5)
title("Akima and image.plot and contour")


#================================================================
# Data in spatial representation

xpoly <- c(x, x[1])
ypoly <- c(y, y[1])
xy <- cbind(xpoly, ypoly)
xypoly <- Polygon(xy, hole=TRUE)
sss1 <- Polygons(list(xypoly), "s1")

xxpoly <- c(xbox, xbox[1])
yypoly <- c(ybox, ybox[1])
# Define the map area polygon
maparea <- cbind(xxpoly, yypoly)
mappoly <- Polygon(maparea)
sss2 <- Polygons(list(mappoly), "s2")

SpP = SpatialPolygons(list(sss2,sss1), 1:2)

pol <- cbind(xpoly, ypoly)
addpoly <- SpatialPolygons(list(Polygons(list(Polygon(pol)), 1)))
addpoly <- as(addpoly, "SpatialPolygonsDataFrame")
addpoly at data[1,1] <- 10
r <- raster(akima$z, xmn=min(akima$x), xmx=max(akima$x), ymn=min(akima$y),
ymx=max(akima$y))
r2 <- rasterize(addpoly, r, field=1, update=TRUE, updateValue="NA")
#image(r2, main="")
plot(r2, main=" ")
contour(akima$x, akima$y, akima$z, add=TRUE, labcex=1.2)
lines(x,y, col="purple", lwd=5)
title("SpatialPolygonsDataFrame, Akima, plot and contour")






--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Clearing-the-area-outside-of-a-polygon-defined-on-a-grid-tp7581718.html
Sent from the R-sig-geo mailing list archive at Nabble.com.



More information about the R-sig-Geo mailing list