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@coords[,1] }))
y <- do.call("c",
sapply(ps,
function(p) { p@coords[,2] }))
id.lengths <- sapply(ps, function(p) { nrow(p@coords) })
# Generate vertex set lengths
list(x=x, y=y, id.lengths=id.lengths)
}
path <- polygonsToPath(map@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@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
p...@stat.auckland.ac.nz
http://www.stat.auckland.ac.nz/~paul/
______________________________________________
R-help@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.
