graticules for projected data

That really is a neat piece of code. Very interested - the result
looks very much like a ggplot2 plot, which I'm a big fan.

Cheers,

Pierre

2010/12/9 Josh London <josh.london at noaa.gov>:

Hello

One of the features I often find myself wanting within spplot is the ability to produce graticules (with corresponding axis labels) that represent latitude/longitude lines for projected data. I often find myself searching the R archives and google for examples and methods for accomplishing this, but have always come up empty.

This has been partially possible by using the gridlines() function within sp and then projecting those spatial lines to the desired projection. There are two drawbacks to this solution:

1) Often the extent of the projected gridlines is smaller than the data of interest so the lines stop short of the plot borders

2) There is no documented method I've found for producing the axis labels with degree units for each of the gridlines. gridat() is close but only supports unprojected data

So, I ventured forth and pulled together a few functions that attempt to provide a generalized solution in the hopes that a) there's not some existing, obvious, solution my endless searching has overlooked and b) other folks could take advantage of it and maybe expand/improve on it.

I've pulled everything together in a draft package, ProjectedGridlines, and it is available on GitHub (https://github.com/jmlondon/ProjectedGridlines). I've got some basic draft documentation in there for those interested.

I make no claims to the code being elegant or an approach that will work for everyone. The main purpose here was to get something that works and distribute it.

Below, I include the three functions and a reproducible example based on the north carolina dataset included with the maptools package. I've copied over some of the roxygen documentation, but those interested in more details should check the files on GitHub.

Cheers
Josh

GetGratLines<-function(spobj,proj_str,exp_deg,ndisc=500,...){
? ? ? ?spobj at bbox[1,1]<-spobj at bbox[1,1]+exp_deg$l
? ? ? ?spobj at bbox[1,2]<-spobj at bbox[1,2]+exp_deg$r
? ? ? ?spobj at bbox[2,1]<-spobj at bbox[2,1]+exp_deg$b
? ? ? ?spobj at bbox[2,2]<-spobj at bbox[2,2]+exp_deg$t
? ? ? ?gl<-spTransform(gridlines(spobj,ndisc=ndisc,...),CRS(proj_str))
? ? ? ?return(gl)
}

#' @param spobj an unrpojected sp object that is the basis of the final projected outcome
#' @param proj_str a proj4 string specifying the desired projection
#' @param exp_deg is a named list with values l,r,b,t specifying the amount in decimal
#' degrees each side (left, right, bottom, top) should be expanded to account for differences
#' in extent between the projected gridlines and the final projected extent of spobj
#' @param ... passing of parameters relevant to gridlines()

GetDegreeLabels<-function(spobj){
? ? ? ?ewl<-coordinates(gridlines(spobj)["EW"]@lines[[1]])
? ? ? ?nsl<-coordinates(gridlines(spobj)["NS"]@lines[[1]])
? ? ? ?ewLabels<-lapply(ewl,function(x) {return(x[1,1])})
? ? ? ?nsLabels<-lapply(nsl,function(x) {return(x[1,2])})
? ? ? ?return(list(EW_Labels=unlist(ewLabels),NS_Labels=unlist(nsLabels)))
}

#' @param spobj an unrpojected sp object that is the basis of the final projected outcome

GetLabelCoords<-function(gl,offset=list(ew=0,ns=0)){
? ? ? ?ewc<-coordinates(gl["EW"]@lines[[1]])
? ? ? ?nsc<-coordinates(gl["NS"]@lines[[1]])
? ? ? ?ewCoords<-lapply(ewc,function(x) {return(x[1,1])})
? ? ? ?nsCoords<-lapply(nsc,function(x) {return(x[1,2])})
? ? ? ?ewCoords<-lapply(ewCoords,function(x) {return(x+offset$ew)})
? ? ? ?nsCoords<-lapply(nsCoords,function(x) {return(x+offset$ns)})
? ? ? ?return(list(EW_Coords=unlist(ewCoords),NS_Coords=unlist(nsCoords)))
}

#' @param gl the SpatialLines object returned from GetGratLines
#' @param offset a named list with values ew and ns specifying an offset for all labels on the axis

And an example:

library(sp)
library(rgdal)
library(maptools)
library(RColorBrewer)

#read in the maptools north carolina shapefile, unprojected
d <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], proj4string=CRS("+proj=longlat +datum=NAD27"))

# create a dummy factor variable, with equal labels:
set.seed(31)
d$f = factor(sample(1:5,100,replace=T),labels=letters[1:5])

#specify a custom Albers equal centered on the coordinate data
#this works for points and polygons, but lines need to be converted to points
med.lat <- round(median(coordinates(d)[,2]))
med.lon <- round(median(coordinates(d)[,1]))
lat1<-bbox(d)[2,1]+diff(range(bbox(d)[2,]))/3
lat2<-bbox(d)[2,2]-diff(range(bbox(d)[2,]))/3

myProj <- paste("+proj=aea +lat_1=",lat1," +lat_2=",lat2," +lat_0=",med.lat,
? ? ? ? ? ? ? ?" +lon_0=", med.lon, "+x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m",sep="")

#project the data, the graticule lines. note the gridlines function calls the unprojected data
d.proj <- spTransform(d,CRS(myProj))
#grat.lines<-spTransform(gridlines(d,ndisc=500),CRS(myProj))
g<-GetGratLines(d,myProj,list(t=0.25,b=-0.25,l=-0.25,r=0.25))
#create the sp.layout component for the graticule lines
grat.plot<-list("sp.lines",g,lwd="1.5",col="white",first=TRUE)

x_grat_at<-GetLabelCoords(g,offset=list(ew=0))$EW_Coords

y_grat_at<-GetLabelCoords(g,offset=list(ns=0))$NS_Coords

x_grat<-list(draw=TRUE,
? ? ? ? ? ? ? ?at=x_grat_at,
? ? ? ? ? ? ? ?labels=as.character(GetDegreeLabels(d)$EW_Labels),
? ? ? ? ? ? ? ?tck=c(0,0),
? ? ? ? ? ? ? ?cex=0.75)

y_grat<-list(draw=TRUE,
? ? ? ? ? ? ? ?at=y_grat_at,
? ? ? ? ? ? ? ?labels=as.character(GetDegreeLabels(d)$NS_Labels),
? ? ? ? ? ? ? ?rot=90,
? ? ? ? ? ? ? ?tck=c(0,0),
? ? ? ? ? ? ? ?cex=0.75)

spplot(d.proj, c("f"),
? ? ? ? ? ? ? ?panel=function(...){
? ? ? ? ? ? ? ? ? ? ? ?panel.fill(col="gray90")
? ? ? ? ? ? ? ? ? ? ? ?panel.polygonsplot(...)
? ? ? ? ? ? ? ?},
? ? ? ? ? ? ? ?sp.layout=list(grat.plot),
? ? ? ? ? ? ? ?col.regions=brewer.pal(5, "Set3"),scales=list(x=x_grat,y=y_grat))





? ? ? ?[[alternative HTML version deleted]]