Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

R is becoming pretty strong tool for handling and analyzing spatial data. I learned some useful things through questions like these at SO and thought it might be useful to have something simmilar, but more 'spatially' oriented.

Can you share some spatial R tips and tricks that you found useful?

share

locked by PolyGeo Jan 15 at 1:20

This question's answers are a collaborative effort: if you see something that can be improved, just edit the answer to improve it! No additional answers can be added here

5  
community wiki? – relet Nov 9 '10 at 9:52
1  
@relet: meta.gis.stackexchange.com/questions/191/… – radek Nov 10 '10 at 11:41
    
ack. Thanks. I believe I saw that one, but it's easy to forget. – relet Nov 10 '10 at 13:15
3  
Also: rspatialtips.org.uk – radek Apr 23 '11 at 12:51
    
I think this should be a community wiki too as this has become a list of sorts. – R.K. Oct 15 '12 at 7:08

11 Answers 11

up vote 28 down vote accepted

This is not so much a trick as it is spplot()'s nifty built-in functionality. spplot()’s ability to scale legend swatches (to match classification break ranges) serves as a useful pedagogical tool when discussing attribute data distribution and classification types. Combining cumulative distribution plots with the maps helps in this endeavor.

enter image description here enter image description here

The students only need to modify a few script parameters to explore classification types and data transformation effects. This is usually their first foray into R in what is a mostly ArcGIS centric course.

Here's a code snippet:

library(rgdal) # Loads SP package by default
NE = readOGR(".", "NewEngland") # Creates a SpatialPolygonsDataFrame class (sp)

library(classInt)
library(RColorBrewer)
pal = brewer.pal(7,"Greens")
brks.qt = classIntervals(NE$Frac_Bach, n = 7, style = "quantile")
brks.jk = classIntervals(NE$Frac_Bach, n = 7, style = "jenks")
brks.eq = classIntervals(NE$Frac_Bach, n = 7, style = "equal")

# Example of one of the map plots
spplot(NE, "Frac_Bach",at=brks.eq$brks,col.regions=pal, col="transparent",
       main = list(label="Equal breaks"))
# Example of one of the cumulative dist plots
plot(brks.eq,pal=pal,main="Equal Breaks")

Ref: Applied Spatial Data Analysis with R (R. Bivand, E Pebesma & V. Gomez-Rubio)

share
1  
very cool! That looks very useful. – djq Nov 29 '11 at 17:22
1  
Thank you very much! Very cool example indeed. – Andy W Nov 29 '11 at 19:38

I was pretty happy to find the dismo package with geocoding and google maps download:

library(dismo)
x <- geocode('110 George Street, Bathurst, NSW, Australia')
a <- x[5:8] + c(-0.001, 0.001, -0.001, 0.001)
e <- extent(as.numeric(a))
g <- gmap(e, type = "satellite")

plot(g)

That is in R 2.12.0 on Windows, it's trivial to install dismo and its dependencies there, not sure on other systems.

alt text

share
    
+1 - Excellent! – radek Nov 25 '10 at 9:28
1  
This looks very useful - however, I'm running into problems with the line e <- extent(x[4:7] + c(-0.001, 0.001, -0.001, 0.001)) ad I get an error message Error: c("x", "y") %in% names(x) is not all TRUE. x[4:7] seems fine though; any thoughts on what the problem might be? – djq Nov 25 '11 at 17:07
    
Yes you need a reproducible example – mdsumner Nov 25 '11 at 19:39
    
I'm trying to reproduce the example in this answer and it does not work. x <- geocode('110 George Street, Bathurst, NSW, Australia') returns ZERO_RESULTS for example, and when I use an example that returns a lat/long, the function e <- extent(x[4:7] + c(-0.001, 0.001, -0.001, 0.001)) also fails. – djq Nov 27 '11 at 21:18
2  
The following also works: e <- extent(as.numeric(x[4:7]) + c(-0.001, 0.001, -0.001, 0.001)) – snth Nov 28 '11 at 8:50

Here are some I have collected.

cran.r spatial, gmane blog, Time series analysis, geodatacenter

share

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

    
Thanks Brad. Useful stuff, especially geodacenter. – radek Nov 10 '10 at 11:34

Also not a trick but here are some resources/examples I have collected

An example of plotting small multiple maps of Areal data in R using the lattice package.

There are a few questions on StackOverflow asking about mapping and R, and here is one with a nice example. I would look at the other answers and the resources they give (as well as search for some more examples) on SO as well.

A different link to the same r-sig-geo group Brad already gave. It is very active, and Roger Bivand answers questions practically every day on the group. Both related to programming and statistical analysis.

Besides checking out the cran spatial page I would also suggest specifically checking out the Spatstat page maintained by Adrian Baddeley. Plenty of examples, a course, and a forthcoming e-book. (At the moment I have been going through the spatstat course, and I think it is a much gentler introduction than the Bivand book).

Not a free resource, but for anyone interested in R I would suggest you check out the Use R! series by Springer. The have a book Applied Spatial Data Analysis with R directly pertinent (also the book A Beginner's Guide to R is my suggested learning R book.)

A free e-book, A Practical Guide to Geostatistical Mapping (Hengl 2009), has examples of applied geostats in R, GRASS, and Google Earth (KML).

If I find anymore good examples I will continue to update (I hope other people post good examples as well!)

share
    
Thanks Andy. I like lattice example. Totally agree on Bivand et al. book - great resource. – radek Nov 10 '10 at 11:35

For raster analysis the raster package is extremely powerful. Beside the standard manual there a few vignettes to get started.

share

I'm not a PostGIS user, but after suggesting Voronoi polygons for a nearest neighbor question, I did a bit of searching. I found that with R, you can create Voronoi polygons for PostGIS. I'm impressed.

share
    
Thanks Kirk - neat trick indeed! – radek Nov 10 '10 at 11:37
    
I'm sure there are easier ways than what I am going to suggest, but you can create Tesselations in the spatstat package, and then convert that tess object to a sp object using this function provided by Adrian Baddeley. From an sp object you can export it to a shapefile if you want. – Andy W Nov 30 '11 at 14:42

See also here how to enjoy high quality statistic analysis in GRASS: http://grass.osgeo.org/wiki/R_statistics

share

I stumbled upon Spatial-Analyst.net. Very informative, comprehensive, and useful. More specific to this question and inline with some of the previous answers, see this page.

share

Example of Point pattern analysis:

#Load library
library(spatstat) 
#create some coordinates        
x=c(78,120,150,17,20,402) 
#prepare the window range      
y=c(70,103,100,205,200,301)
win=owin(range(x),range(y)) 
#create the point pattern
p <- ppp(x,y,window=win)
#Plot it
plot(p) 

Creates a point pattern and depict it. The spatstat package has a number of functions for analysing geographical data. Here are some spatstat tutorials:

share

With this function, you can easily make spatial joins, but only if all areas are filled by polygons.

library(rgeos)
library(sp) 
library(maptools)
library(rgdal)
library(sp)
xy.map <- readShapeSpatial("http://www.udec.cl/~jbustosm/points.shp")
manzana.map <- readShapeSpatial("http://www.udec.cl/~jbustosm/manzanas_from.shp" )

IntersectPtWithPoly <- function(x, y) { 
# Extracts values from a SpatialPolygonDataFrame with SpatialPointsDataFrame, and appends table (similar to 
# ArcGIS intersect)
# Args: 
#   x: SpatialPoints*Frame
#   y: SpatialPolygonsDataFrame
# Returns:
# SpatialPointsDataFrame with appended table of polygon attributes

  # Set up overlay with new column of join IDs in x
  z <- overlay(y, x)

  # Bind captured data to points dataframe
  x2 <- cbind(x, z)

  # Make it back into a SpatialPointsDataFrame 
  # Account for different coordinate variable names 
  if(("coords.x1" %in% colnames(x2)) & ("coords.x2" %in% colnames(x2))) {
    coordinates(x2) <- ~coords.x1 + coords.x2  
  } else if(("x" %in% colnames(x2)) & ("x" %in% colnames(x2))) {
    coordinates(x2) <- ~x + y 
  }

  # Reassign its projection if it has one
  if(is.na(CRSargs(x@proj4string)) == "FALSE") {
    x2@proj4string <- x@proj4string  
  }
  return(x2)
}


test<-IntersectPtWithPoly (xy.map,manzana.map)
share

Not sure if this qualifies as a "trick", but I'm a huge fan of the combination of the acs package (for selecting US Census data) and the leaflet package (for making interactive javascript maps that can be hosted online).

This tutorial does an excellent job illustrating the benefit of using these two packages together.

share

Not the answer you're looking for? Browse other questions tagged or ask your own question.