Method for: Creating contour maps of geocoded point shapefile attributes

Hello:

If one has a polygon shapefile of the USA and a point shapefile of USA 
cities (both in geographic coordinates), and the point shapefile 
includes demographic attributes for each city, what would be the best 
way to create a contour map of one of the attributes that overlays the 
polygon map?

I have been experimenting with using akima.interp (to generate a grid of 
interpolated values) and contour (to generate the contours), with mixed 
results. I think that I can work out a technique to get this done, but 
would like to know if this capability is embedded in any of the R 
spatial classes.

Not in the spatial classes, and the issue of Great Circle distances is 
there too. You could look at inverse distance weighting (in a small 
neighbourhood in gstat I think using GC distances (?idw, ?gstat and its 
predict method), or at kriging on GC distances in fields (?Krig, 
Distance="rdist.earth" - Tps does not seem to support GC distances).

This is just a sketch:

library(gstat)
library(fields)
data(ozone2)
o <- cbind(ozone2$lon.lat, ozone2$y[16,])
colnames(o) <- c("lon", "lat", "day16")
o <- as.data.frame(o)
o <- o[complete.cases(o),]
coordinates(o) <- c(1, 2)
proj4string(o) <- CRS("+proj=longlat")
grd <- GridTopology(c(-94,36), c(0.25,0.25), c(48,40))
SG <- SpatialGrid(grd, proj4string=CRS("+proj=longlat"))
p <- idw(day16 ~ 1, o, newdata=SG, nmax=16)
image(p[1], axes=TRUE)
points(o, pch=3)
pim <- as.image.SpatialGridDataFrame(p[1])
cl <- contourLines(pim)
library(maptools)
cL <- ContourLines2SLDF(cl)
proj4string(cL) <- CRS("+proj=longlat")
lines(cL)

Roger

Thanks for any insight into this!

Rick Reeves

Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no