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
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo