Brought up by @adrfantini in #48.
This is a first and very rough shot at proxy objects in stars. They contain only the dimension table and a pointer to the data (now: the file name, we could add the GDALDataset pointer if this helps perfomance-wise). Based on this we can do subsetting, cropping, and thinning; only when we plot those cells that we can see as pixels are read. The speedup for an Sentinel-2 image (4 10-m bands, 10000 x 10000) over complete reading and thinning in memory is a factor 100 or so (at thinning factor 43).
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3
tif = system.file("tif/L7_ETMs.tif", package = "stars")
s2 = "SENTINEL2_L1C:/vsizip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:10m:EPSG_32632"
p = read_stars(s2, proxy = TRUE)
p
# stars_proxy object with 1 attribute in file:
# $`MTD_MSIL1C.xml:10m:EPSG_32632`
# [1] "SENTINEL2_L1C:/vsizip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:10m:EPSG_32632"
# dimension(s):
# from to offset delta refsys point values
# x 1 10980 3e+05 10 +proj=utm +zone=32 +datum... NA NULL
# y 1 10980 6e+06 -10 +proj=utm +zone=32 +datum... NA NULL
# bands 1 4 NA NA NA NA NULL
st_crs(p)
# Coordinate Reference System:
# EPSG: 32632
# proj4string: "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
st_dimensions(p)
# from to offset delta refsys point values
# x 1 10980 3e+05 10 +proj=utm +zone=32 +datum... NA NULL
# y 1 10980 6e+06 -10 +proj=utm +zone=32 +datum... NA NULL
# bands 1 4 NA NA NA NA NULL
#p2 = p[,1:4000, 40:4060,]
p2 = p
stars:::get_downsample(dim(p2))
# [1] 43
png("s2.png", 1000, 1000)
system.time(plot(p2, rgb = 2:4))
# user system elapsed
# 0.970 0.208 0.942
system.time(plot(read_stars(s2), rgb = 2:4))
# user system elapsed
# 82.905 5.385 36.679
edzer
All 22 comments
The idea is that when we do for instance an st_apply, we can choose a read-write strategy to go through the array without completely reading, or keeping complete results, in memory.
Lazy evaluation means that an st_apply call would do nothing, only when we plot we would carry out the functions for those pixels shown, only when we do a write / download call would the st_apply be carried out to all the cells.
edzer
on 29 Aug 2018
I did not test yet and can't do so until monday, but this looks spot on. Does this already work with stars objects which are collages of multiple files?
adrfantini
on 30 Aug 2018
I think not.
edzer
on 30 Aug 2018
This now works with e.g. multiple netcdf files forming a time series, but not with spatial collages (mosaics?) of files; maybe GDAL VRT is the way to go here, if we work from disk?
See here for a new vignette showcasing the proxy objects, @mdsumner - comments very welcome, planning to do an r-spatial blog in this soon, and then make remote objects work over http connections (with R on both sides).
The vignette is now runnable if you download starsdata, a 1 Gb package in my drat repo (instructions in vignette).
edzer
on 10 Sep 2018
Ah, I also interfaced the GDAL resampling methods; see vignette.
edzer
on 10 Sep 2018
@appelmar
edzer
on 10 Sep 2018
Thanks! Awesome to have starsdata too
mdsumner
on 10 Sep 2018
Yes, and I'd be happy to add examples to that. Too bad it doesn't fit on github, with those bulky files (although, TBH I didn't try, so far).
edzer
on 10 Sep 2018
Might be an idea to bump sf and depend on the version? It caught me out briefly so may affect other early adopters
mdsumner
on 10 Sep 2018
I just noted that I hadn't committed my latest sf changes to github, you have to reinstall it once more. I now bumped the sf version, and required version of sf by stars.
edzer
on 10 Sep 2018
Stars happily reads in full and plots a 25m Tasmania DEM in less than 30s fwiw ;)
mdsumner
on 10 Sep 2018
Did you read in the foll dataset first, or plot the proxy object? I wonder because the S2 image took 1s to plot.
edzer
on 10 Sep 2018
Oh I meant full, before I realized my sf was old - I know the proxy will be instant!
mdsumner
on 10 Sep 2018
On Windows I see a problem, I can't even find where that "extent" code occurs? On Ubuntu all is well :)
library(stars)
tif = system.file("tif/L7_ETMs.tif", package = "stars")
rasterio = list(nXOff = 6, nYOff = 6, nXSize = 100, nYSize = 100, bands = c(1,3,4))
(x = read_stars(tif, RasterIO = rasterio))
Error in CPL_read_gdal(as.character(x), as.character(options), as.character(driver), :
Expecting a single value: [extent=0].
session_info fwiw
Session info ------------------------------------------------------------------------------------------------------------------------
setting value
version R version 3.5.1 Patched (2018-08-31 r75224)
system x86_64, mingw32
ui RStudio (1.1.456)
language (EN)
collate English_Australia.1252
tz Australia/Hobart
date 2018-09-11
Packages ----------------------------------------------------------------------------------------------------------------------------
package * version date source
abind * 1.4-5 2016-07-21 CRAN (R 3.5.0)
assertthat 0.2.0 2017-04-11 CRAN (R 3.5.1)
base * 3.5.1 2018-09-02 local
bindr 0.1.1 2018-03-13 CRAN (R 3.5.1)
bindrcpp 0.2.2 2018-03-29 CRAN (R 3.5.1)
class 7.3-14 2015-08-30 CRAN (R 3.5.1)
classInt 0.2-3 2018-04-16 CRAN (R 3.5.1)
compiler 3.5.1 2018-09-02 local
crayon 1.3.4 2017-09-16 CRAN (R 3.5.1)
datasets * 3.5.1 2018-09-02 local
DBI 1.0.0 2018-05-02 CRAN (R 3.5.1)
devtools 1.13.6 2018-06-27 CRAN (R 3.5.1)
digest 0.6.16 2018-08-22 CRAN (R 3.5.1)
dplyr 0.7.6 2018-06-29 CRAN (R 3.5.1)
e1071 1.7-0 2018-07-28 CRAN (R 3.5.1)
glue 1.3.0 2018-07-17 CRAN (R 3.5.1)
graphics * 3.5.1 2018-09-02 local
grDevices * 3.5.1 2018-09-02 local
grid 3.5.1 2018-09-02 local
magrittr 1.5 2014-11-22 CRAN (R 3.5.1)
memoise 1.1.0 2017-04-21 CRAN (R 3.5.1)
methods * 3.5.1 2018-09-02 local
ncdf4 1.16 2017-04-01 CRAN (R 3.5.0)
ncmeta 0.0.2 2018-02-13 CRAN (R 3.5.1)
pillar 1.3.0 2018-07-14 CRAN (R 3.5.1)
pkgconfig 2.0.2 2018-08-16 CRAN (R 3.5.1)
purrr 0.2.5 2018-05-29 CRAN (R 3.5.1)
R6 2.2.2 2017-06-17 CRAN (R 3.5.1)
Rcpp 0.12.18 2018-07-23 CRAN (R 3.5.1)
rlang 0.2.2 2018-08-16 CRAN (R 3.5.1)
RNetCDF 1.9-1 2017-10-08 CRAN (R 3.5.0)
rstudioapi 0.7 2017-09-07 CRAN (R 3.5.1)
sf * 0.7-0 2018-09-10 Github (r-spatial/sf@7c94ca6)
spData 0.2.9.3 2018-08-01 CRAN (R 3.5.1)
stars * 0.1-1 2018-09-10 local
stats * 3.5.1 2018-09-02 local
tibble 1.4.2 2018-01-22 CRAN (R 3.5.1)
tidyselect 0.2.4 2018-02-26 CRAN (R 3.5.1)
tools 3.5.1 2018-09-02 local
units 0.6-0 2018-06-09 CRAN (R 3.5.1)
utils * 3.5.1 2018-09-02 local
withr 2.1.2 2018-03-15 CRAN (R 3.5.1)
yaml 2.2.0 2018-07-25 CRAN (R 3.5.1)
CPL_read_gdal resides in sf. Have you tried reinstalling sf from github on Windows?
edzer
on 11 Sep 2018
Sorry, got caught two commits behind - it's fine!
mdsumner
on 11 Sep 2018
Is there any feasible way to make it behave more like raster, in which calls modifying the dimensions get immediately reflected on the dataset? One example will explain this better.
stars workflow:
library(stars)
(s = read_stars(somefile, proxy=TRUE))
#stars_proxy object with 1 attribute in file:
#$hmr.nc
#[1] "/home/clima-archive4-b/afantini/chym_simulations/analysis/rain/hmr/hmr.nc"
#
#dimension(s):
# from to offset delta refsys point values
#x 1 223 NA NA NA NA NULL [x]
#y 1 248 NA NA NA NA NULL [y]
#height 1 1 0 m NA NA NA NULL
#time 1 5114 2000-01-01 1 days POSIXct NA NULL
#affine parameters: NA NA
# Say I want to remove one redundant dimension, then slice on time, then check the dimensions:
(s %>% adrop)[,,,1:10] %>% st_dimensions
# from to offset delta refsys point values
#x 1 223 NA NA NA NA NULL [x]
#y 1 248 NA NA NA NA NULL [y]
#height 1 1 0 m NA NA NA NULL
#time 1 5114 2000-01-01 1 days POSIXct NA NULL
#affine parameters: NA NA
As you can see it's not showing the currently selected dimensions.
With raster:
library(raster)
(r = brick(somefile))
#class : RasterBrick
#dimensions : 248, 223, 55304, 5114 (nrow, ncol, ncell, nlayers)
#resolution : 1, 1 (x, y)
#extent : 0.5, 223.5, 0.5, 248.5 (xmin, xmax, ymin, ymax)
#coord. ref. : NA
#data source : /home/clima-archive4-b/afantini/chym_simulations/analysis/rain/hmr/hmr.nc
#names : X2000.01.01.00.00.00, X2000.01.02.00.00.00, X2000.01.03.00.00.00, X2000.01.04.00.00.00, X2000.01.05.00.00.00, X2000.01.06.00.00.00, X2000.01.07.00.00.00, X2000.01.08.00.00.00, X2000.01.09.00.00.00, X2000.01.10.00.00.00, X2000.01.11.00.00.00, X2000.01.12.00.00.00, X2000.01.13.00.00.00, X2000.01.14.00.00.00, X2000.01.15.00.00.00, ...
#Date/time : 2000-01-01 00:00:00, 2013-12-31 00:00:00 (min, max)
#varname : pr
#level : 1
r[[1:10]]
#class : RasterStack
#dimensions : 248, 223, 55304, 10 (nrow, ncol, ncell, nlayers)
#resolution : 1, 1 (x, y)
#extent : 0.5, 223.5, 0.5, 248.5 (xmin, xmax, ymin, ymax)
#coord. ref. : NA
#names : X2000.01.01.00.00.00, X2000.01.02.00.00.00, X2000.01.03.00.00.00, X2000.01.04.00.00.00, X2000.01.05.00.00.00, X2000.01.06.00.00.00, X2000.01.07.00.00.00, X2000.01.08.00.00.00, X2000.01.09.00.00.00, X2000.01.10.00.00.00
Seems the raster approach is a bit more user-friendly there.
adrfantini
on 18 Sep 2018
I think it is conceptually different, and it might not harm to understand the concepts behind what you are doing.
(s %>% adrop)[,,,1:10] %>% st_as_stars %>% st_dimensions
should tell you the dimensions of the data when fetched by st_as_stars.
edzer
on 30 Sep 2018
Did you remove read_stars_meta as a result of the introduction of this feature? Just to be sure.
adrfantini
on 3 Oct 2018
Yes.
edzer
on 3 Oct 2018
Ok, thanks.
adrfantini
on 3 Oct 2018
To verify https://github.com/r-spatial/stars/issues/63, otherwise to be closed.
edzer
on 10 Oct 2018
Related issues
kendonB
路
6Comments
noamross
路
10Comments
michaeldorman
路
11Comments
edzer
路
12Comments
adrfantini
路
10Comments