[R-sig-Geo] graticules for projected data

Hadley Wickham hadley at rice.edu
Mon Dec 13 19:45:31 CET 2010 

It should be possible to use RGEOS instead, but I'm not sure how far
along that is.
Hadley

On Thu, Dec 9, 2010 at 11:04 PM, Pierre Roudier
<pierre.roudier at gmail.com> wrote:
> Harley,
>
> In ggplot2 however, is there a way to plot shapefiles without gpclib?
> Its licence is a problem for any externaly-founded research
> unfortunately.
>
> Cheers,
>
> Pierre
>
> 2010/12/10 Hadley Wickham <hadley at rice.edu>:
>> 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/
>>
>
>



-- 
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
http://had.co.nz/

More information about the R-sig-Geo mailing list