If the problem is memory handling in R (still a big headache to do much of GIS analysis in R), you should instead consider reading pieces of grids that you need e.g.:
info <- GDALinfo("dem50m.asc") info gridmaps01 = readGDAL("dem50m.asc", region.dim = round(c(info[["rows"]]/2,info[["columns"]]))) gridmaps02 = readGDAL("dem50m.asc", region.dim = round(c(info[["rows"]]/2,info[["columns"]])), offset=round(c(info[["rows"]]/2,0))) will split the original map into two tiles. You could practically read any 'tile'/subset of a grid using the 'region.dim' argument (rows and column coordinates). You could first prepare an empty grid topology, then overlay it over your polygons just to get the region dimension, and then read each piece of the grid that you need. In GIS terms, getting a raster from one grid topology to (any) other grid topology (also subgrids) is referred to 'resampling' (I do not think that this is yet implemented in any R package; a combination of akima and spTransform will do the trick but it is not really compact). A generic function to combine R+ILWIS to put a map to a different grid (e.g. to LatLon coordinates) would be: # write to ILWIS: > writeGDAL(dem25m[1], "dem25m.mpr", "ILWIS") # create a new grid: ... > gridparameters(geoarc) # resample the map (Bilinear) to the new geographic grid: > shell(cmd=paste(ILWIS, " crgrf geoarc.grf ",[EMAIL PROTECTED]@cells.dim[[2]]," ",[EMAIL PROTECTED]@cells.dim[[1]]," -crdsys=LatlonWGS84 -lowleft=(",[EMAIL PROTECTED]@cellcentre.offset[[1]],",",[EMAIL PROTECTED]@cellcentre.offset[[2]],") -pixsize=",[EMAIL PROTECTED]@cellsize[[1]],sep=""), wait=F) > shell(cmd=paste(ILWIS, "dem25m_ll_c.mpr = MapResample(dem25m_c.mpr, geoarc, > BiLinear)"), wait=F) > shell(cmd=paste(ILWIS, "open dem25m_ll_c.mpr -noask"), wait=F) see http://geomorphometry.org/R.asp hth, Tom Hengl http://spatial-analyst.net -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of Pieter Beck Sent: Thursday, October 23, 2008 9:15 PM To: r-sig-geo@stat.math.ethz.ch Subject: [R-sig-Geo] subsetting Dear all, I am a new user to the world of spatial applications of R and am wondering about this trivial question: What is the easiest way to subset a SpatialGridDataFrame, based on 4 corner coordinates? The reason I ask is that I want to do an overlay of a polygon file on a large raster dataset, to extract summary statistics from it. Both files are extremely large, however, and just applying overlay causes memory trouble. I was hoping to loop through the polygons, subset the SpatialGridDataFrame using the bounding box surrounding a polygon and then create the summary statistics one polygon at-a-time. To do this, I'd need to subset the SpatialGridDataFrame based on the bounding box corner coordinates. Thanks in advance for the help, Kind Regards, Pieter Beck [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo