Ok, fair enough that there's no magic involved. I've worked through the details with the small example as follows. The result is only a couple of cells different in each direction.
library(stars)
tif <- system.file("tif/L7_ETMs.tif", package = "stars")
x <- read_stars(tif)
x <- x[,,,1]
# calculate a circular polygon at the center of the raster
pol <- x %>% st_bbox() %>% st_as_sfc() %>% st_centroid() %>% st_buffer(500)
bb_pol <- st_bbox(pol)
xoff <- st_dimensions(x)$x$offset
xdelt <- st_dimensions(x)$x$delta
yoff <- st_dimensions(x)$y$offset
ydelt <- st_dimensions(x)$y$delta
cropXoff <- (bb_pol$xmin - xoff + xdelt)/xdelt
cropXsize <- (bb_pol$xmax - bb_pol$xmin)/xdelt
cropYoff <- (bb_pol$ymin - yoff + ydelt)/ydelt
cropYsize <- (bb_pol$ymax - bb_pol$ymin)/ydelt
# if ydelt is negative, get abs of ysize and move yoffset to the top
if(cropYsize < 0){
cropYsize <- abs(cropYsize)
cropYoff <- cropYoff - cropYsize
}
rasterio <- list(nXOff = cropXoff, nYOff = cropYoff, nXSize = cropXsize, nYSize = cropYsize, bands = c(1))
(z = read_stars(tif, RasterIO = rasterio))
plot(z)
# crop it if desired
plot(z[pol])
## compare to proxy method
x <- read_stars(tif, proxy = TRUE)
x <- x[,,,1]
y <- st_as_stars(x[pol])
plot(y)
# only a couple of cells off!
z
y
------------------------------
Date: Tue, 13 Nov 2018 17:26:16 +0100
From: Edzer Pebesma <edzer.pebesma at uni-muenster.de>
To: r-sig-geo at r-project.org
Subject: Re: [R-sig-Geo] stars::RasterIO using extent info?
Message-ID: <9d037da7-dc9c-9886-d6fc-5864cf8b4a17 at uni-muenster.de>
Content-Type: text/plain; charset="utf-8"
On 11/13/18 4:10 PM, Howard, Tim G (DEC) via R-sig-Geo wrote:
Dear list, I am exploring the different options for reading parts of large imagery object in stars, as discussed here: https://r-spatial.github.io/stars/articles/proxy.html My ultimate goal is to read into RAM only a clipped portion of a large raster (well, actually a raster stack, but taking baby steps here). My immediate question: the `RasterIO` option of read_stars defines cell offsets and cell counts (*Size). Is there a straightforward way to calculate these values given extent information? Reproducible example (mostly taken from here: https://www.r-spatial.org/r/2018/03/22/stars2.html): library(stars) tif <- system.file("tif/L7_ETMs.tif", package = "stars") x <- read_stars(tif) # read entire tif into ram x <- x[,,,1] #get just one layer for now # calculate a circular polygon at the center of the raster pol <- x %>% st_bbox() %>% st_as_sfc() %>% st_centroid() %>% st_buffer(500) plot(x) # interestingly, I don't think the circle is in the right place when plotted plot(st_geometry(pol), add = TRUE, border = "red") # this is what I'd like to be able to restrict to what is read in memory: plot(x[pol]) ## read only portion of tif using proxy object x <- read_stars(tif, proxy = TRUE) x <- x[,,,1] y <- st_as_stars(x[pol]) plot(y) # this is cropped to the extent (but not the circle - let's not worry about that right now) Question: can I do the equivalent with the RasterIO options in stars?? Said another way, instead of setting up the proxy, can I map my extent object (or bounding box) directly to the cell count values needed for RasterIO?
stars can do the math, and so can you; it is explained here: https://r-spatial.github.io/stars/articles/data_model.html stars uses some functions directly from GDAL which it doesn't expose to the user, but there is no magic going on here.
Thanks in advance for any tips. Tim
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081