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

2018-11-15 Thread Howard, Tim G (DEC) via R-sig-Geo
Mike, 
Thanks for the feedback - I'll check it out. 
Best, 
Tim

From: Michael Sumner [mailto:mdsum...@gmail.com] 
Sent: Wednesday, November 14, 2018 5:55 PM
To: Howard, Tim G (DEC) 
Cc: r-sig-geo@r-project.org
Subject: Re: [R-sig-Geo] stars::RasterIO using extent info?

FWIW I have a raster-delivering-front-end for RasterIO in this dev package:  

https://github.com/hypertidy/lazyraster

It uses the more obvious extent() idioms and will even use an open plot if 
nothing else is specified. 

(It uses an independent binding to GDAL in the vapour package). 

That might help, or not. It's on my list to find a sensible way for stars to 
leverage this obvious ease-of-use.  rgdal::readGDAL always had an interface to 
RasterIO, but only raster ever made use of that, and it's indirect - raster 
doesn't allow a mix of extent and resolution in its approach. 

Cheers, Mike. 

On Thu, 15 Nov 2018 at 02:27 Howard, Tim G (DEC) via R-sig-Geo 
<mailto:r-sig-geo@r-project.org> wrote:
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 <mailto:edzer.pebe...@uni-muenster.de>
To: mailto:r-sig-geo@r-project.org
Subject: Re: [R-sig-Geo] stars::RasterIO using extent info?
Message-ID: <mailto:9d037da7-dc9c-9886-d6fc-5864cf8b4...@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
>
> _

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

2018-11-14 Thread Michael Sumner
FWIW I have a raster-delivering-front-end for RasterIO in this dev package:

https://github.com/hypertidy/lazyraster

It uses the more obvious extent() idioms and will even use an open plot if
nothing else is specified.

(It uses an independent binding to GDAL in the vapour package).

That might help, or not. It's on my list to find a sensible way for stars
to leverage this obvious ease-of-use.  rgdal::readGDAL always had an
interface to RasterIO, but only raster ever made use of that, and it's
indirect - raster doesn't allow a mix of extent and resolution in its
approach.

Cheers, Mike.

On Thu, 15 Nov 2018 at 02:27 Howard, Tim G (DEC) via R-sig-Geo <
r-sig-geo@r-project.org> wrote:

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

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

2018-11-14 Thread Howard, Tim G (DEC) via R-sig-Geo
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 
To: r-sig-geo@r-project.org
Subject: Re: [R-sig-Geo] stars::RasterIO using extent info?
Message-ID: <9d037da7-dc9c-9886-d6fc-5864cf8b4...@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@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

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


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

2018-11-13 Thread Edzer Pebesma


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


pEpkey.asc
Description: application/pgp-keys
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


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

2018-11-13 Thread Howard, Tim G (DEC) via R-sig-Geo
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?


Thanks in advance for any tips. 
Tim

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo