Convert esri shape files to esri ascii grid from R

On Fri, 10 Mar 2006, Mulholland, Tom wrote:

 

I didn't know how to do this offhand, but a little research soon found a
way. So I think you need to search the existing resources before asking
your question as my basic search from R

   

RSiteSearch("esri grid") 
     

returned write.asciigrid as the first response. Since I routinely read
shapefiles into sp classes (readShapePoints in maptools) it's the
conversion of one to another.

The overview PDF file in the sp package gives an example on converting
points to grids ("5.2 Creating grids from points".)

It would appear that you just need to familarise yourself with what's in
sp and maptools. I would suggest looking through previous posts on this
list with regard to shapefiles which will give you an understanding of
some of the limitations that you might bump into. For example (I don't
know if it is still true) I have some point shapefiles that don't seem
to be currently supported.
   

Thanks, Tom, there are ways in sp/maptools, but not all of the 
possibilities have been explored - it turns out that there are (many) more 
than we thought of in developing sp classes.

Some notes from trying this out:

library(maptools)
nc <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], 
 IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
ncbb <- bbox(nc)
xymin <- ncbb[,1]
r_xy <- apply(ncbb, 1, diff)
resol <- 0.1 # here taken as ew=ns
ncells <- ceiling(r_xy/resol)
grd <- GridTopology(cellcentre.offset=xymin, cellsize=c(resol, resol), 
 cells.dim=ncells)
pts <- SpatialPoints(coordinates(grd), proj4string=CRS(proj4string(nc)))
pts_polys <- overlay(pts, nc)
pts_BIR74 <- nc$BIR74[pts_polys]
pts1 <- SpatialPointsDataFrame(pts, data=as(data.frame(BIR74=pts_BIR74), 
 "AttributeList"))
gridded(pts1) <- TRUE
pal <- colorRampPalette(c("wheat1", "red3"))(4)
brks <- quantile(nc$BIR74)
par(mfrow=c(3,1))
plot(nc, col=pal[findInterval(nc$BIR74, brks, all.inside=TRUE)], 
 axes=TRUE)
image(pts1, col=pal, breaks=brks, axes=TRUE)
plot(nc, add=TRUE)
fname <- tempfile()
writeAsciiGrid(pts1, fname) # note that this will not work now because
 # machine fuzz in the calculated resolution gets in the way, fixed in 
 # next maptools release, for now, the kludge is:
 # slot(slot(pts1, "grid"), "cellsize") <- rep(mean(slot(slot(pts1, 
 # "grid"), "cellsize")), 2)
rpts1 <- readAsciiGrid(fname, proj4string=CRS(proj4string(pts1)))
image(rpts1, col=pal, breaks=brks, axes=TRUE)
par(mfrow=c(1,1))

The steps are to make a GridTopology object from the input 
SpatialPolygons, get its coordinates, overlay() the coordinates on the 
polygons, retrieve the required data - could be multiple variables, adding 
it to the SpatialPoints object, grid the object, and write out.

In the new version of rgdal, you'll find that writing out rasters in other 
formats is well-supported, including their spatial reference system 
metadata - the same applies to reading vector data (but not yet writing 
vector data).

I guess that a function to make GridTopology objects from SpatialPolygons 
may enter maptools before long.

Comments, suggestions?

Roger