[R] Clip a contour with shapefile while using contourplot

[Paul Murrell]

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/