[R-sig-Geo] stars::RasterIO using extent info?

[Howard, Tim G (DEC)]

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 using uni-muenster.de>
To: r-sig-geo using r-project.org
Subject: Re: [R-sig-Geo] stars::RasterIO using extent info?
Message-ID: <9d037da7-dc9c-9886-d6fc-5864cf8b4a17 using 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 using 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