[R-sig-Geo] Can I read only a portion of a raster? (With more code)

[Michael Sumner]

Use crop() once it's been "read":

r <- raster(system.file("external/test.grd", package="raster"))
my.ext <- extent(179000, 181000, 330000, 333000)
r <- crop(r, my.ext)

There's no "ext" argument to raster() for reading from a file. There
is that argument when you create a raster from scratch with no data,
as an alternative to xmin/ymin/xmax/ymax.

You can read a grid partially with rgdal, by either
1) using the index arguments to readGDAL (offset, region.dim) - but
you have to determine those index values yourself
2) using GDAL.open to create a handle to the file that can be indexed
with [,], but again you need to figure out the indexes

But, raster is pretty good at opening a file and not attempting to
read the whole thing. I would just try doing it as above with crop.

Cheers, Mike.




On Thu, Sep 26, 2013 at 4:35 PM, Ben Harrison
<harb at student.unimelb.edu.au> wrote:
> I thought I could use the extent of a raster to trim it, but it ignores the
> option. Is there another way to do it?
>
> library(raster)
>
> r <- raster(system.file("external/test.grd", package="raster"))
> print(r)
>
> # =-=-=-=- Results:
> class       : RasterLayer
> dimensions  : 115, 80, 9200  (nrow, ncol, ncell)
> resolution  : 40, 40  (x, y)
> extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
> coord. ref. : +init=epsg:28992
> +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812
> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079
> +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs
> data source :
> /home/blah/R/x86_64-pc-linux-gnu-library/3.0/raster/external/test.grd
> names       : test
> values      : 128.434, 1805.78  (min, max)
> # =-=-=-=-
>
> plot(r)
>
> my.ext <- extent(179000, 181000, 330000, 333000)
> # note: extent in relation to r's extent is smaller.
> q <- raster(system.file("external/test.grd", package="raster"),
>             ext=my.ext)
> print(q)
>
> # =-=-=-=- Results:
> class       : RasterLayer
> dimensions  : 115, 80, 9200  (nrow, ncol, ncell)
> resolution  : 40, 40  (x, y)
> extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
> coord. ref. : +init=epsg:28992
> +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812
> +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079
> +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs
> data source :
> /home/blah/R/x86_64-pc-linux-gnu-library/3.0/raster/external/test.grd
> names       : test
> values      : 128.434, 1805.78  (min, max)
> # =-=-=-=-
>
> plot(q)
> all.equal(target=r, current=q)
> [1] TRUE
>
> It reads in the entire file. How can I read in only the portion I need?
>
> Ben.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



-- 
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com