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

[Roger Bivand]

On Fri, 10 Mar 2006, Edzer J. Pebesma wrote:

> Roger Bivand wrote:
> 
> >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

pts1 <- SpatialPixelsDataFrame(pts, data=as(data.frame(BIR74=pts_BIR74), 
  "AttributeList"))

does it in one step, yes.

> >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
> >  
> >
> Two comments: first the step "grid the object" sounds obsolete
> somehow, second: the function you mention should perhaps
> enter sp, not maptools.

That's fair, I agree.

Thanks,

Roger

> --
> Edzer
> 

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