Clip a contour with shapefile while using contourplot
Hi
Below is some code that does what I think you want by drawing a path
based on the map data. This does some grubby low-level work with the
'sp' objects that someone else may be able to tidy up
# The 21st polygon in 'hello' is the big outer boundary
# PLUS the 20 other inner holes
map <- as(hello, "SpatialPolygons")[21]
# Convert map outline to path information
polygonsToPath <- function(ps) {
# Turn the list of polygons into a single set of x/y
x <- do.call("c",
sapply(ps,
function(p) { p at coords[,1] }))
y <- do.call("c",
sapply(ps,
function(p) { p at coords[,2] }))
id.lengths <- sapply(ps, function(p) { nrow(p at coords) })
# Generate vertex set lengths
list(x=x, y=y, id.lengths=id.lengths)
}
path <- polygonsToPath(map at polygons[[1]]@Polygons)
# Generate rect surrounding the path
xrange <- range(path$x)
yrange <- range(path$y)
xbox <- xrange + c(-50000, 50000)
ybox <- yrange + c(-50000, 50000)
# Invert the path
invertpath <- list(x=c(xbox, rev(xbox), path$x),
y=c(rep(ybox, each=2), path$y),
id.lengths=c(4, path$id.lengths))
# Draw inverted map over contourplot
contourplot(Salinity ~ Eastings+Northings | Time, mydata,
cuts=30, pretty=TRUE,
panel=function(...) {
panel.contourplot(...)
grid.path(invertpath$x, invertpath$y,
id.lengths=invertpath$id.lengths,
default="native",
gp=gpar(col="green", fill="white"))
})
The final result is far from perfect, but I hope it might be of some help.
One issue is that most of the contour labels are obscured, though that
might be ameliorated by filling the inverted map with a semi-transparent
colour like rgb(1,1,1,.8).
Paul
On 15/02/13 08:58, Janesh Devkota wrote:
Hi, I have done the interpolation for my data and I was able to create the contours in multipanel with the help of Pascal. Now, I want to clip the contour with the shapefile. I want only the portion of contour to be displayed which falls inside the boundary of the shapefile. The data mydata.csv can be found on https://www.dropbox.com/s/khi7nv0160hi68p/mydata.csv The data for shapefile can be found on https://www.dropbox.com/sh/ztvmibsslr9ocmc/YOtiwB8p9p THe code I have used so far is as follows: # Load Libraries library(latticeExtra) library(sp) library(rgdal) library(lattice) library(gridExtra) #Read Shapefile hello <- readOGR("shape", layer="Export_Output_4") ## Project the shapefile to the UTM 16 zone proj4string(hello) <- CRS("+proj=utm +zone=16 +ellps=WGS84") ## Read Contour data mydata <- read.csv("mydata.csv") head(mydata ) summary(mydata) #Create a contourplot contourplot(Salinity ~ Eastings+Northings | Time, mydata, cuts=30,pretty=TRUE) Thank you so much. I would welcome any other ways to do this aside from contourplot and lattice. Best Regards, Janesh [[alternative HTML version deleted]]
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Dr Paul Murrell Department of Statistics The University of Auckland Private Bag 92019 Auckland New Zealand 64 9 3737599 x85392 paul at stat.auckland.ac.nz http://www.stat.auckland.ac.nz/~paul/