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

[Edzer J. Pebesma]

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