[R-sig-Geo] subsetting

[Tomislav Hengl]

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 ",geoarc at grid@cells.dim[[2]],"
",geoarc at grid@cells.dim[[1]]," -crdsys=LatlonWGS84
-lowleft=(",geoarc at grid@cellcentre.offset[[1]],",",geoarc at grid@cellcentre.offset[[2]],")
-pixsize=",geoarc at grid@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: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of
Pieter Beck
Sent: Thursday, October 23, 2008 9:15 PM
To: r-sig-geo at 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 at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo