[R-sig-Geo] Clipping a large image without reading in the entire image

[Roger Bivand]

On Fri, 1 Jun 2007, Andrew Niccolai wrote:

> Greetings Spatial R users/developers,
> 
> I was curious as to whether or not it would be possible to describe a
> bounding box for a large image so that what is read into R would be just the
> contents of the image inside the bounding box.  Conceptually, it would be
> akin to "clipping" the image but without the image read into the program.
> My issue is that I have a very large 3-band orthophoto of several hundreds
> of hectares of forest (roughly 8MB image).  I am running an analysis of
> LIDAR data that requires the user to select a smaller area of a highly
> degraded lidar image.  I would like the bounding box of this LIDAR area to
> be used to "clip" a georegistered orthophoto that complements the LIDAR
> data.  However, the complete ortho image seems to crash my WindowsXP, R
> v2.4.1 setup. When I take the bounding box coordinates into ERDAS Imagine
> and subset the orthophoto the following code works fine to bring it into R:
> 
> ## Read in Tif file
> setwd(anypathway)
> tif.image <- "Pre_Subset_Image.tif"
> LIDAR.tif <- GDAL.open(tif.image)
> dim(LIDAR.tif)
> image(LIDAR.tif[,,1])

Have you looked at all of the arguments to the various rgdal functions? 
Here you are using the "[" method for GDALReadOnlyDataset objects, but

?"[,GDALReadOnlyDataset-method"

says

?getRasterData

which gives as mich detail in subsetting as is available in GDAL itself,
including: band = NULL, offset = c(0, 0), region.dim = dim(dataset),
output.dim = region.dim, interleave = c(0, 0). In your case offset= and 
region.dim= should do the trick. My reading of the code of the "[" method 
is that it will take these arguments and pass them through, so:

LIDAR.tif[,,1, offset=c(0,0), region.dim=c(50,50)]

should pick a 50x50 tile starting at (0,0), typically NW corner. This is 
how I understand Tim Keitt's reply.

Roger

> 
> ## This converts the rgdal class to an SpatialGridDataFrame class
> LIDAR.tif.1 <- LIDAR.tif[,]
> bbox(LIDAR.tif.1)
> image(LIDAR.tif.1,main="")
> mtext("Orthophoto of Preloaded Fusion Example", side=3, line = 3, font=3,
> cex=1.25)
> 
> 
> However, this would be impractical with the code that I am building to jump
> out of R and subset images in Imagine (plus very expensive for the end
> user).  Therefore, does anyone have any solutions for reading in an image
> that would work in theory like the following:
> 
> setwd(anypathway)
> tif.image <- "VERY_LARGE_IMAGE.tif"
> Bounded.tif <- GDAL.open(tif.image, bbox=c(UpL.X, LwR.X, UpL.Y, LwR.Y))
> 
> By the way, I am not set on GDAL.open as the solution but it seems to work
> well reading in Tifs & Jpegs.
> 
> Again, thanks in advance
> 
> 
> Andrew Niccolai
> Doctoral Candidate
> Yale School of Forestry
> 
> 

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