[R-sig-Geo] Convert esri shape files to esri ascii grid from R

[Roger Bivand]

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


> 
> Tom
> 
> > -----Original Message-----
> > From: r-sig-geo-bounces at stat.math.ethz.ch
> > [mailto:r-sig-geo-bounces at stat.math.ethz.ch]On Behalf Of Yong Li
> > Sent: Friday, 10 March 2006 7:18 AM
> > To: r-sig-geo at stat.math.ethz.ch
> > Subject: [R-sig-Geo] Convert esri shape files to esri ascii 
> > grid from R
> > 
> > 
> > Hi R spatial experts,
> > 
> > Has anybody done the conversion of esri shape files to esri ascii 
> > grid based on field values in R?
> > I knew it can be done easily in Arcview using avenue. But if it can 
> > be done in R, that will fly my daily work loads.
> > 
> > Regards
> > 
> > Yong
> > 
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
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