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.roud...@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.lon...@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,...){ >>> sp...@bbox[1,1]<-sp...@bbox[1,1]+exp_deg$l >>> sp...@bbox[1,2]<-sp...@bbox[1,2]+exp_deg$r >>> sp...@bbox[2,1]<-sp...@bbox[2,1]+exp_deg$b >>> sp...@bbox[2,2]<-sp...@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@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@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/ _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
