[R-sig-Geo] graticules for projected data

[Josh London]

Hadley

As you can tell from the final example, I really appreciate the look and feel of ggplot2 and use it for most of my non-map plots. However, I have, maybe wrongly, presumed fortify only works for spatial polygons and lines. I often have a need to map spatial points, pixels or grids and so have stuck with spplot for most of my mapping.

cheers
Josh
 
On Dec 9, 2010, at 3:36 PM, Hadley Wickham wrote:

> ggplot2 + coord_map() already does graticules in a similar manner, as
> far as I can tell.
> 
> Hadley
> 
> On Thu, Dec 9, 2010 at 2:22 PM, Pierre Roudier <pierre.roudier at gmail.com> wrote:
>> 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]]
>>> 
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> 
>> 
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> 
> 
> 
> 
> -- 
> Assistant Professor / Dobelman Family Junior Chair
> Department of Statistics / Rice University
> http://had.co.nz/