For SpatialPoints* and SpatialPixels* I'm using ggplot2: library(sp) library(ggplot2)
data(meuse) data(meuse.grid) coordinates(meuse) <- ~x+y coordinates(meuse.grid) <- ~x+y # ... some spatial processing plot_meuse <- ggplot(as.data.frame(meuse)) + geom_point(aes(x,y, colour=zinc)) + coord_equal() plot_meuse.grid <- ggplot(as.data.frame(meuse.grid)) + geom_tile(aes(x,y, fill=dist)) + coord_equal() # You can alo combine both geoms plot_combined <- ggplot() + geom_tile(data=as.data.frame(meuse.grid), aes(x,y, fill=factor(soil))) + geom_point(data=as.data.frame(meuse), aes(x,y, colour=zinc)) + coord_equal() HTH, Pierre 2010/12/14 Josh London <josh.lon...@noaa.gov>: > 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
