Skip to content

Clearing the area outside of a polygon defined on a grid.

6 messages · Robert J. Hijmans, peleve, Roger Bivand

peleve
# 
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.
peleve
# 
Thank you very much.  That worked, mostly!?!  But I see the methodology now.
I had been thinking about masks, but I was asking the wrong question I now
realize.
In my real application, using your approach, I still have some unwanted bits
left outside of the polygon.  I traced it as a left over piece from
contouring the data... contours are vector objects. Your suggestion works
for the raster object,  and worked for the Akima result (outside the polygon
and inside the convex polygon that Akima works in).
I now realize that what I need is a hole in a polygon, outside colored
white, inside transparent, or a polygon with everything outside white and
inside transparent. The hole is the outline of the area with the contours
within the bounding box.
Any thoughts on how to do that?
Thank you again.


I think you can do:

x <- raster(akima)
y <- mask( x, SpP[2] )
plot(y)






--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Clearing-the-area-outside-of-a-polygon-defined-on-a-grid-tp7581718p7581721.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
peleve
# 
Thanks Robert,
I was able to do it with the following using the 'gpclib" package:

p2 <- as(xy, "gpc.poly")
p1 <- as(maparea, "gpc.poly")
plot(setdiff(p1, p2), poly.args = list(col = 3), add = TRUE)

Did exactly what I was looking for.  Pretty simple, actually.

Thanks for inspiring me.

Pete




--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Clearing-the-area-outside-of-a-polygon-defined-on-a-grid-tp7581718p7581727.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
Roger Bivand
#
On Thu, 15 Nov 2012, peleve wrote:

            

      
      
      
    

        
For completeness, please use rgeos rather than gpclib, which we are 
advising against strongly, because of the unclear and very restrictive 
licence of its included C code. rgeos contains gpclib classes and methods 
using GEOS rather than GPC under the hood, and should give identical 
results. It would be very helpful if you could confirm that the same code 
run in rgeos gives the same results. Once you get to rgeos, you may then 
find that there are efficiency gains in using rgeos wrappers to GEOS 
functions directly.

Roger