stars_proxy - HOWTO extract a time series at selected pixel
Hello,
I wish to extract all values, at a specific pixel, along a dimension using a stars_proxy object.
I deal with very large space-time datasets, and I need to extract "on-demand" all values along the time dimension for a specific pixel.
I know how to perform such operation when all the data are read in memory, but I am searching a way to perform such operation using stars_proxy objects, and without the need to read all the data in memory using st_as_stars().
plot.stars_proxy function is able to extract such data blazingly fast in the background, without reading all the data into memory. Therefore, I tried to imitate its back-end operations to extract quickly time series for a specific pixel.
Here below I post some trials that I made, but none of these methods allows to extract the original values that I obtain from a stars read-in-memory object (or through the call to st_as_stars( , downsample=0)).
Am I missing something obvious?
Thanks in advance for the help !
```
library(stars)
library(dplyr)
library(starsdata) # install.packages("starsdata", repos = "http://pebesma.staff.ifgi.de", type = "source")
Read sentinel data
granule = system.file("sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip", package = "starsdata")
s2 = paste0("SENTINEL2_L1C:/vsizip/", granule, "/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:10m:EPSG_32632")
p = read_stars(s2, proxy = TRUE) # read in memory
r = read_stars(s2, proxy = FALSE) # proxy object
names(p)
p
Plot is able to extract super quickly the data
- Requires at least 2 values in first 2 dimensions !!!
p1 <- p[i=TRUE, 2:3, 2:3,1:4]
plot(p1)
plot.stars_proxy = function(x, y, ..., downsample = get_downsample(dim(x))) {
x = st_as_stars(x, downsample = downsample, ...)
plot(x, ..., downsample = 0)
}
Extract time series on the fly
- Imitate plot.stars_proxy to extract the time series without reading the data into memory
row_ID = 10
col_ID = 20
Reference 1 - Extraction from a star object with all data already in memory
drop(pull(r[i=TRUE, row_ID, col_ID,1:4]))
r %>% slice(x, 10) %>% slice(y,20) %>% pull
Reference 2 - Extraction from a star proxy object
- Following the guidelines provided in stars vignettes
- Too slow ...
p3 <- p[i=TRUE, row_ID, col_ID,1:4]
p3 <- st_as_stars(p3)
drop(pull(p3))
Version 1 - Extraction from a star proxy object
- Trying to imitate plot.stars_proxy
- The downsample operation prevent to retrieve the correct values
downsample = stars:::get_downsample(dim(p))
x = st_as_stars(p, downsample = downsample)
drop(pull(x[i=TRUE, row_ID, col_ID, 1:4]))
Version 2 - Extraction from a star proxy object
- Trying to imitate plot.stars_proxy
- The downsample operation prevent to retrieve the correct values
p2 <- p[i=TRUE, row_ID, col_ID, 1:4]
downsample1 = stars:::get_downsample(dim(p2))
x1 = st_as_stars(p2, downsample = downsample1)
drop(pull(x1))
Version 3 (Error !!!)
p4 <- p[i=TRUE, row_ID, col_ID,1:4]
pull(p4) # Error !!!
Version 4 (Error !!!)
- In the call list, only the last call is recorded !!!
p5 <- p %>% slice(x, 10) %>% slice(y,20)
p5
pull(p5)
````
ghiggi
All 26 comments
Here is something that works, e.g. for retrieving band values for 100 randomly sampled locations from a proxy object:
s = st_sample(st_as_sfc(st_bbox(p)), 100)
a = aggregate(p, s, FUN = mean)
This takes over 20 seconds on my laptop, but it is not looking at the overviews (like plot, or downsample in general, would) but gets to the real 10m resolution. This means unzip and look into the highres .jp2 files.
I'd agree that for this common case, we'd need another verb - st_extract?
edzer
on 29 Apr 2020
Thanks! Nice workaround!
I post here below a solution, with a temporary st_extract() accepting row and column index as input arguments.
I think it could be very useful to have a st_extract verb. As input arguments would be nice to accept either xy_coords, rowcol_idx, or cell_number for compatibility with the raster package objects.
This function opens me a lot of doors for the interactive visualization, investigation, and intercomparison of space-time models output with shiny!
I guess that a verb called st_extract_patch could also be very useful for extracting image patches as well as train and deploy predictive models for large datasets.
````
library(stars)
library(dplyr)
library(starsdata) # install.packages("starsdata", repos = "http://pebesma.staff.ifgi.de", type = "source")
Read sentinel stars
granule = system.file("sentinel/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.zip", package = "starsdata")
s2 = paste0("SENTINEL2_L1C:/vsizip/", granule, "/S2A_MSIL1C_20180220T105051_N0206_R051_T32ULE_20180221T134037.SAFE/MTD_MSIL1C.xml:10m:EPSG_32632")
p = read_stars(s2, proxy = TRUE)
r = read_stars(s2, proxy = FALSE)
Define extraction
st_extract <- function(p, row_ID, col_ID) {
l_dims <- st_dimensions(p)
x_val = (col_ID-l_dims$x$from)l_dims$x$delta + l_dims$x$offset + l_dims$x$delta0.5
y_val = (row_ID-l_dims$y$from)l_dims$y$delta + l_dims$y$offset + l_dims$y$delta0.5
# if (isTRUE(cell_midpoints)) {
# x_val <- x_val - l_dims$x$delta0.5
# y_val <- y_val - l_dims$y$delta0.5
# }
s = st_sfc(st_point(c(x_val, y_val), dim="XY"))
return(drop(pull(aggregate(p, s, FUN = identity))))
}
Testing
row_ID = 1
col_ID = 1
st_extract(p, row_ID, col_ID)
drop(pull(r[i=TRUE, row_ID, col_ID,1:4]))
row_ID = 10
col_ID = 5
st_extract(p, row_ID, col_ID)
drop(pull(r[i=TRUE, col_ID, row_ID ,1:4]))
Check out of bounds case
row_ID = 10981
col_ID = 10981
st_extract(p, row_ID, col_ID)
drop(pull(r[i=TRUE, col_ID, row_ID ,1:4]))
```
ghiggi
on 29 Apr 2020
I think you can already get a patch by x[patch], or using st_crop. raster::extract seems to extract by spatial objects, not by row/col or cell ID. You can get row/col subsets by x[,row,col] also for ranges. Ah, I now see that raster::extract does cell numbers; would that really be a bonus?
edzer
on 29 Apr 2020
In raster you can extract all the values along the third dimension by doing x[cell_number].
When I have statistical models fitted to each pixel time series, I usually save such models:
- in named list, with the name corresponding to the ID (cell number) of the raster object.
- or in list columns of sf or dataframe objects, where each row of the data.frame corresponds to a pixel or spatial shape. In the case of gridded data, an ID column records the cell number of the raster object.
When processing the time series I often store the time series (in zoo format) as list-columns in (sf) data.frame objects (I have a dfTSlist (sf) data.frame and TSlist class for this), with a column providing the cell number of the raster object.
So I guess the bonus of cell numbers is that if someone want to use list or data.frame to save the models, it can simply put the cell numbers as the list_element names or create a column ID in the sf or data.frame object.... although it could still create and ID string of the style ID_<x>_<y> :)
Maybe thinking about other packages like spex and fasterize, if I remember correctly they return the cell number of the raster object ;)
I have a package that I am developing that performs all such data format conversion across sf-rasters-zoo-data.frame (dfTSlist, df-wide and df-long format) and now stars objects formats, keeping trace of measurement units and the time resolution of the time series; and distributing the computations parallelly across the cores using futures and/or parallel-plyr.
Regarding st_extract_patch, I was thinking more to a function taking as input arguments patch size and returning some sort of iterator (in python style) that allow sliding easily over all the image in a for loop and perform the operations: read - preprocess - fit and save the model and read-apply model-save data. It could be also very useful for people working on Keras and Tensorflow in R that currently use ad-hoc R6 structures to process image/volumetric data in batches.
ghiggi
on 29 Apr 2020
There is one difficulty: stars has arbitrary n-D arrays where space can be 1, 2 or 3 dimensions, and arbitrary many follow (or preceed). How would cell number reflect cells, and recycle in this space? For the other packages you mention the model is confined to a raster stack, where this is trivial (always x/y=raster, always recycle over the stack).
edzer
on 29 Apr 2020
I am not sure I fully understood your concerns.
The cell_number does not refer to the cell_index. The cell_number in raster objects is just defined on the spatial dimensions. Does stars distinguish between the spatial dimensions, the time dimension and the data dimensions through refsys right? Because the cell_number would be defined only on the spatial dimensions and is not an index that is extended to the time and data dimension (like a cell_index).
Using a cell number arguments, this would generalize automatically to 1D, 2D or 3D spatial dimensions. Oppositely, asking to specify the pixel coordinates would require specific arguments based on the space dimension (like xy for 2D or xyz for 3D) or a coord argument whose length must correspond to the spatial dimensions of the stars object.
Just as a thought, the creation of an attribute describing the dimensionality of the dataset could be useful to implement ad-hoc functions (i.e: space_dim: x, y z ; time_dim: t, data_dim: bands, realizazions_dim: ensembles). By specifying these dimensions when loading the data, the user may have access to specialized functions that will be part of separated packages.
ghiggi
on 1 May 2020
Right now, st_redimension lets you collapse spatial dimensions into one, then you can do your cell_number stuff on that collapsed dimension, and st_redimension back to the original one.
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
x1 = read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))
(r = st_redimension(x1, c(prod(dim(x1)[1:2]), 6)))
# stars object with 2 dimensions and 1 attribute
# attribute(s):
# L7_ETMs.tif
# Min. : 1.00
# 1st Qu.: 54.00
# Median : 69.00
# Mean : 68.91
# 3rd Qu.: 86.00
# Max. :255.00
# dimension(s):
# from to offset delta refsys point values
# X1 1 122848 NA NA NA NA NULL
# X2 1 6 NA NA NA NA NULL
# do your cell_number stuff as index on the first dimension;
# work result back to original dimension:
st_redimension(r, st_dimensions(x1))
# stars object with 3 dimensions and 1 attribute
# attribute(s):
# L7_ETMs.tif
# Min. : 1.00
# 1st Qu.: 54.00
# Median : 69.00
# Mean : 68.91
# 3rd Qu.: 86.00
# Max. :255.00
# dimension(s):
# from to offset delta refsys point values
# x 1 349 288776 28.5 UTM Zone 25, Southern Hem... FALSE NULL [x]
# y 1 352 9120761 -28.5 UTM Zone 25, Southern Hem... FALSE NULL [y]
# band 1 6 NA NA NA NA NULL
Is that another way of solving your problem?
edzer
on 1 May 2020
I considered to use st_redimension, but I guess that for extracting the values from the stars_proxy object I still have to pass from generating an sfc point object (with lat lon coordinates) and then use aggregate() %>% pull()
Additionally, when subsetting by coordinates, if these are in lon/lat, st_intersect now output the warning message although coordinates are longitude/latitude, st_intersectsassumes that they are planar. I guess that by using cell_number (or indexes of the spatial dimensions) the extraction would be much faster and would avoid the current inaccuracy related to the planar coordinates assumption of st_intersect!
ghiggi
on 1 May 2020
st_intersects using points to find matching cells in a regular grid and setting as_points=FALSE is very fast; there are no polygons involved. The warning is mostly harmless unless you grid covers poles or anti-meridians (and should disappear when as_points=FALSE is set).
edzer
on 1 May 2020
May I jump into this discussion. Fell free to kick me out if this is not related. I guess that I have a similar problem. My raster data is spatially coarse (say 32km x 32km), but temporally long. Thus, I use stars_proxy objects for those. Back in the days, I did interpolations from y (raster) to x (SpatialPoints) with
raster::extract(x, y, method = 'bilinear')
And I would love to do the same with x (stars_proxy), and y (sf or sfc with point geometries),
st_extract(x, y, method = 'bilinear')
But, I am not sure where to start. Now, following this discussion it seems to me that building st_extract for stars_proxy objects is not so easy because of their _arbitrary n-D arrays_.
Do I get this right?
meteosimon
on 5 May 2020
Thanks for jumping in @meteosimon! I think your use case is very common and needs a solid implementation in stars.
edzer
on 5 May 2020
Here (branch: extract) is a first shot at st_extract, which is now a wrapper around GDAL's warp function (no crs checking so far):
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif)
pts = st_sample(st_as_sfc(st_bbox(x)), 10)
st_extract(x, pts)
# Simple feature collection with 10 features and 6 fields
# geometry type: POINT
# dimension: XY
# bbox: xmin: 289089.3 ymin: 9111650 xmax: 297390.4 ymax: 9119211
# projected CRS: UTM Zone 25, Southern Hemisphere
# V1 V2 V3 V4 V5 V6 geometry
# 1 107 102 100 23 20 15 POINT (295007.7 9111791)
# 2 90 81 59 13 13 12 POINT (295245.2 9111650)
# 3 75 63 66 72 107 63 POINT (290656.7 9119211)
# 4 61 48 34 74 72 36 POINT (292504.8 9118447)
# 5 83 70 78 73 123 96 POINT (296360 9115012)
# 6 87 74 73 66 82 63 POINT (297390.4 9115722)
# 7 83 70 78 73 123 96 POINT (296364.8 9115010)
# 8 76 59 62 46 107 87 POINT (289200.5 9112814)
# 9 72 59 61 57 90 72 POINT (295940.9 9113339)
# 10 74 60 61 44 89 73 POINT (289089.3 9112768)
st_extract(x, pts, method = 'bilinear')
# Simple feature collection with 10 features and 6 fields
# geometry type: POINT
# dimension: XY
# bbox: xmin: 289089.3 ymin: 9111650 xmax: 297390.4 ymax: 9119211
# projected CRS: UTM Zone 25, Southern Hemisphere
# V1 V2 V3 V4 V5 V6
# 1 101.37649 95.87314 91.74907 20.48054 17.75872 14.08422
# 2 90.56324 82.04916 59.15646 12.91596 13.48591 11.12565
# 3 70.64030 59.12567 58.20448 76.13236 98.77301 55.15506
# 4 59.86362 46.33345 34.19751 70.95416 67.33197 34.06101
# 5 84.45757 72.80767 82.63710 70.98868 125.14644 97.82825
# 6 87.64106 74.25022 77.22578 62.84208 88.80307 69.23398
# 7 85.21119 72.72086 82.36779 71.40185 124.63238 97.69682
# 8 76.34825 60.02158 63.69492 47.77619 108.01000 88.85746
# 9 71.67777 58.56748 60.75430 58.13129 92.44106 72.81389
# 10 74.89028 60.75696 60.32932 45.36544 93.77626 78.18030
# geometry
# 1 POINT (295007.7 9111791)
# 2 POINT (295245.2 9111650)
# 3 POINT (290656.7 9119211)
# 4 POINT (292504.8 9118447)
# 5 POINT (296360 9115012)
# 6 POINT (297390.4 9115722)
# 7 POINT (296364.8 9115010)
# 8 POINT (289200.5 9112814)
# 9 POINT (295940.9 9113339)
# 10 POINT (289089.3 9112768)
Awesome!! I already assumed that there should be a gdal-way, the warper to a tiny cell seems reasonable. However, my test case data has 3 dimension (lon, lat, time), and thus nz within st_extract() is set to the length of the time dimension, here
nz = ifelse(length(dim(x)) == 2, 1, dim(x)[3])
Is there an easy way to identify the purpose of the dimensions, i.e., space, time, band, ... maybe using st_dimensions(x)$refsys or attr(st_dimensions(x), "raster")[["dimensions"]]?
meteosimon
on 5 May 2020
Ah, so what you want returned is a vector data cube, with two dimensions: point locations and time?Would you mind making an example dataset & script available?
edzer
on 5 May 2020
Sure, for instance, I do this,
tif = system.file("tif/L7_ETMs.tif", package = "stars")
x = read_stars(tif, proxy = TRUE)
pts <- st_sample(st_as_sfc(st_bbox(x)), 2)
xx <- st_extract(x, pts)
and get that
R> xx
Simple feature collection with 2 features and 6 fields
geometry type: POINT
dimension: XY
bbox: xmin: 289774.7 ymin: 9113201 xmax: 291165.8 ymax: 9113698
projected CRS: UTM Zone 25, Southern Hemisphere
band1 band2 band3 band4 band5 band6 geometry
1 62 50 36 89 70 32 POINT (291165.8 9113698)
2 84 73 84 55 133 112 POINT (289774.7 9113201)
but what I actually want is a stars object with 2 dimensions (geometry and band) and 1 attribute.
Do you think it makes sense to have an argument in st_extract() specifying the class of the return value?
meteosimon
on 5 May 2020
[EDITED:] It makes indeed more sense to return a stars object in case there is more than a one-layer raster. If there is a one-layer raster, an sf object is returned.
edzer
on 5 May 2020
Some timings against raster::extract (on a 32Gb RAM machine; not sure how meaningful these are):
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
f = "2019_30m_cdls.img"
x = read_stars(f)
dim(x)
# x y
# 153811 96523
pts = st_sample(st_as_sfc(st_bbox(x)), 100)
system.time(res1 <- st_extract(x, pts, method = 'near'))
# user system elapsed
# 2.188 0.055 2.311
system.time(res2 <- raster::extract(raster::raster(f), as(pts, "Spatial"), method = 'simple'))
# user system elapsed
# 1.157 1.769 5.769
all.equal(res2, res1[[1]])
# [1] TRUE
system.time(res1 <- st_extract(x, pts, method = 'bilinear'))
# user system elapsed
# 2.248 0.051 2.301
system.time(res2 <- raster::extract(raster::raster(f), as(pts, "Spatial"), method = 'bilinear'))
# Error: cannot allocate vector of size 7.7 Gb
# Timing stopped at: 3.185 4.49 7.691
# Execution halted
# Warning message:
# system call failed: Cannot allocate memory
for larger samples (1000+) the raster simple/near method is nearly 2 x faster; the bilinear method has succeeded with raster but much slower than st_extent() and with a massive memory footprint.
edzer
on 5 May 2020
I did two simiilar tests (timing_extract.pdf) for my data (ERA5 obtained via climate data store):
st_extract()vs.raster::extract()Result: Resulting values are not equal, asraster::extract()collapses the non-raster dimension.st_extract()produces data with the right dimensions. (... of course, I rather should have used aRasterStackorRasterBrick, not a plainrasterobject!)st_extract()forstarsobject va.stars_proxyobject. Result: While the values for thestarsobject look reasonable,st_extract()produces totally odd values for thestars_proxyobject.
Any idea why it fails for stars_proxy objects?
Link to my data on the Universitaet Innsbruck ftp: ftp://ftp.uibk.ac.at/private/c4031039_20200512_de55983e249ea5f3dafc38e096a22f33/ERA5_sfc_2018_2m_temperature.nc
meteosimon
on 6 May 2020
It doesn't fail, but gdalwarp works on cell values and ignores the offset and scale metadata fields. If you would transform the netcdf file such that the temperatures were unscaled, Float32 Kelvin values, everything would be fine.
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
rp = read_stars("ERA5_sfc_2018_2m_temperature.nc", proxy = TRUE)
s = st_sample(st_as_sfc(st_bbox(rp)), 100)
ep = st_extract(rp, s)
# [1] "/tmp/RtmpaeA5VM/file2116f48dbbd.tif"
ep$T = 276.5769284958387 + 0.001034976096027598 * ep$ERA5_sfc_2018_2m_temperature.nc
ep
# stars object with 2 dimensions and 2 attributes
# attribute(s):
# ERA5_sfc_2018_2m_temperature.nc T
# Min. :-29753 Min. :245.8
# 1st Qu.: 38 1st Qu.:276.6
# Median : 6592 Median :283.4
# Mean : 6635 Mean :283.4
# 3rd Qu.: 12968 3rd Qu.:290.0
# Max. : 32290 Max. :310.0
# dimension(s):
# from to offset delta refsys point
# sfc 1 100 NA NA NA TRUE
# time 1 8760 2018-01-01 UTC 1 hours POSIXct NA
# values
# sfc POINT (15.92234 46.47825),...,POINT (5.637295 50.55694)
# time NULL
Thank you @edzer for making this clear. I was not aware of this. I did some more testing, and things work fine, except when I have a stars_proxy object encapsulating multiple (netcdf) files. Then I run into a segmentation fault:
library(stars)
x <- c("avhrr-only-v2.19810901.nc", "avhrr-only-v2.19810902.nc")
file_list <- system.file(paste0("netcdf/", x), package = "starsdata")
y <- read_stars(file_list, quiet = TRUE, proxy = TRUE)
pt <- st_sample(st_as_sfc(st_bbox(y)), 2)
e <- st_extract(y, pt)
## *** caught segfault ***
## address (nil), cause 'memory not mapped'
##
## Traceback:
## 1: CPL_gdal_warper(source, destination, as.integer(resampling_method(options)), oo, doo)
## 2: sf::gdal_utils("warper", x[[1]], tmp, method)
## 3: st_extract.stars_proxy(y, pt)
## 4: st_extract(y, pt)
I am not sure how multiple files are resolved in a stars_proxy. Could this error indicate that R cannot find the way to the original data through a path of several GDAL virtual datasets?? (Just speculating!)
meteosimon
on 7 May 2020
The closest that works here is this:
library(stars)
x <- c("avhrr-only-v2.19810901.nc")
file_list <- system.file(paste0("netcdf/", x), package = "starsdata")
file_list = gdal_subdatasets(file_list)
y <- read_stars(file_list[[1]], quiet = TRUE, proxy = TRUE)
pt <- st_sample(st_as_sfc(st_bbox(y)), 2)
(e <- st_extract(y, pt))
i.e., give the subdataset to read_stars, and don't give it more than one.
Yes, read_stars _without_ proxy does a lot of stuff that we'd have to do manually when relying on bare GDAL routines behaving the same...
edzer
on 7 May 2020
This now works:
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
x <- c("avhrr-only-v2.19810901.nc", "avhrr-only-v2.19810902.nc", "avhrr-only-v2.19810903.nc")
file_list <- system.file(paste0("netcdf/", x), package = "starsdata")
(y <- read_stars(file_list[[1]], quiet = TRUE, proxy = TRUE))
# stars_proxy object with 4 attributes in files:
# $sst
# [1] "[...]/avhrr-only-v2.19810901.nc:sst"
# $anom
# [1] "[...]/avhrr-only-v2.19810901.nc:anom"
# $err
# [1] "[...]/avhrr-only-v2.19810901.nc:err"
# $ice
# [1] "[...]/avhrr-only-v2.19810901.nc:ice"
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 1440 0 0.25 NA NA NULL [x]
# y 1 720 90 -0.25 NA NA NULL [y]
# zlev 1 1 0 [m] NA NA NA NULL
# time 1 1 1981-09-01 UTC NA POSIXct NA NULL
(pt <- st_sample(st_as_sfc(st_bbox(y)), 2))
# Geometry set for 2 features
# geometry type: POINT
# dimension: XY
# bbox: xmin: 71.76336 ymin: -79.46659 xmax: 274.4065 ymax: 72.84809
# CRS: NA
# POINT (71.76336 72.84809)
# POINT (274.4065 -79.46659)
(e <- st_extract(y, pt))
# stars object with 4 dimensions and 1 attribute
# attribute(s):
# sst
# Min. : 0.00
# 1st Qu.: 0.00
# Median : 0.00
# Mean : 50.62
# 3rd Qu.: 36.50
# Max. :321.00
# dimension(s):
# from to offset delta refsys point
# sfc 1 2 NA NA NA TRUE
# zlev 1 1 0 [m] NA NA NA
# time 1 1 1981-09-01 UTC NA POSIXct NA
# band 1 4 NA NA NA NA
# values
# sfc POINT (71.7634 72.8481), POINT (274.406 -79.4666)
# zlev NULL
# time NULL
# band sst,...,ice
@meteosimon this was pretty much a complete rewrite of st_extract; I now dropped the ability to do bilinear interpolation, although I can see the use of this; in return we do get handling of complete proxy objects, ncdf time series, units and offset + scale.
edzer
on 6 Oct 2020
That puts back bilinear as (the only) interpolation option for st_extract.
edzer
on 7 Oct 2020
I now see:
library(stars)
# Loading required package: abind
# Loading required package: sf
# Linking to GEOS 3.8.0, GDAL 3.1.3, PROJ 7.1.1
x = c(
"avhrr-only-v2.19810901.nc",
"avhrr-only-v2.19810902.nc",
"avhrr-only-v2.19810903.nc",
"avhrr-only-v2.19810904.nc",
"avhrr-only-v2.19810905.nc",
"avhrr-only-v2.19810906.nc",
"avhrr-only-v2.19810907.nc",
"avhrr-only-v2.19810908.nc",
"avhrr-only-v2.19810909.nc"
)
file_list = system.file(paste0("netcdf/", x), package = "starsdata")
(y = read_stars(file_list, quiet = TRUE))
# stars object with 4 dimensions and 4 attributes
# attribute(s), summary of first 1e+05 cells:
# sst [掳*C] anom [掳*C] err [掳*C] ice [percent]
# Min. :-1.80 Min. :-4.69 Min. :0.110 Min. :0.010
# 1st Qu.:-1.19 1st Qu.:-0.06 1st Qu.:0.300 1st Qu.:0.730
# Median :-1.05 Median : 0.52 Median :0.300 Median :0.830
# Mean :-0.32 Mean : 0.23 Mean :0.295 Mean :0.766
# 3rd Qu.:-0.20 3rd Qu.: 0.71 3rd Qu.:0.300 3rd Qu.:0.870
# Max. : 9.36 Max. : 3.70 Max. :0.480 Max. :1.000
# NA's :13360 NA's :13360 NA's :13360 NA's :27377
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 1440 0 0.25 NA NA NULL [x]
# y 1 720 90 -0.25 NA NA NULL [y]
# zlev 1 1 0 [m] NA NA NA NULL
# time 1 9 1981-09-01 02:00:00 CEST 1 days POSIXct NA NULL
st_crs(y) = "OGC:CRS84"
pts = st_sample(st_as_sfc(st_bbox(y)), 10)
# although coordinates are longitude/latitude, st_intersects assumes that they are planar
# although coordinates are longitude/latitude, st_intersects assumes that they are planar
s1 = st_extract(y, pts)
(y = read_stars(file_list, quiet = TRUE, proxy = TRUE))
# stars_proxy object with 4 attributes in files:
# $sst
# [1] "[...]/avhrr-only-v2.19810901.nc:sst" "[...]/avhrr-only-v2.19810902.nc:sst"
# [3] "[...]/avhrr-only-v2.19810903.nc:sst" "[...]/avhrr-only-v2.19810904.nc:sst"
# [5] "[...]/avhrr-only-v2.19810905.nc:sst" "[...]/avhrr-only-v2.19810906.nc:sst"
# [7] "[...]/avhrr-only-v2.19810907.nc:sst" "[...]/avhrr-only-v2.19810908.nc:sst"
# [9] "[...]/avhrr-only-v2.19810909.nc:sst"
# $anom
# [1] "[...]/avhrr-only-v2.19810901.nc:anom"
# [2] "[...]/avhrr-only-v2.19810902.nc:anom"
# [3] "[...]/avhrr-only-v2.19810903.nc:anom"
# [4] "[...]/avhrr-only-v2.19810904.nc:anom"
# [5] "[...]/avhrr-only-v2.19810905.nc:anom"
# [6] "[...]/avhrr-only-v2.19810906.nc:anom"
# [7] "[...]/avhrr-only-v2.19810907.nc:anom"
# [8] "[...]/avhrr-only-v2.19810908.nc:anom"
# [9] "[...]/avhrr-only-v2.19810909.nc:anom"
# $err
# [1] "[...]/avhrr-only-v2.19810901.nc:err" "[...]/avhrr-only-v2.19810902.nc:err"
# [3] "[...]/avhrr-only-v2.19810903.nc:err" "[...]/avhrr-only-v2.19810904.nc:err"
# [5] "[...]/avhrr-only-v2.19810905.nc:err" "[...]/avhrr-only-v2.19810906.nc:err"
# [7] "[...]/avhrr-only-v2.19810907.nc:err" "[...]/avhrr-only-v2.19810908.nc:err"
# [9] "[...]/avhrr-only-v2.19810909.nc:err"
# $ice
# [1] "[...]/avhrr-only-v2.19810901.nc:ice" "[...]/avhrr-only-v2.19810902.nc:ice"
# [3] "[...]/avhrr-only-v2.19810903.nc:ice" "[...]/avhrr-only-v2.19810904.nc:ice"
# [5] "[...]/avhrr-only-v2.19810905.nc:ice" "[...]/avhrr-only-v2.19810906.nc:ice"
# [7] "[...]/avhrr-only-v2.19810907.nc:ice" "[...]/avhrr-only-v2.19810908.nc:ice"
# [9] "[...]/avhrr-only-v2.19810909.nc:ice"
# dimension(s):
# from to offset delta refsys point values x/y
# x 1 1440 0 0.25 NA NA NULL [x]
# y 1 720 90 -0.25 NA NA NULL [y]
# zlev 1 1 0 [m] NA NA NA NULL
# time 1 9 1981-09-01 02:00:00 CEST 1 days POSIXct NA NULL
st_crs(y) = "OGC:CRS84"
s2 = st_extract(y, pts)
all.equal(s1, s2)
# [1] TRUE
s2
# stars object with 3 dimensions and 4 attributes
# attribute(s):
# sst [掳*C] anom [掳*C] err [掳*C] ice [percent]
# Min. : 1.130 Min. :-2.2600 Min. :0.1100 Min. : NA
# 1st Qu.: 5.475 1st Qu.:-0.5050 1st Qu.:0.2300 1st Qu.: NA
# Median :10.580 Median :-0.1900 Median :0.3300 Median : NA
# Mean :12.514 Mean :-0.3662 Mean :0.3168 Mean :NaN
# 3rd Qu.:19.795 3rd Qu.: 0.0000 3rd Qu.:0.3900 3rd Qu.: NA
# Max. :29.610 Max. : 0.4600 Max. :0.5600 Max. : NA
# NA's :27 NA's :27 NA's :27 NA's :90
# dimension(s):
# from to offset delta refsys point
# geometry 1 10 NA NA WGS 84 TRUE
# zlev 1 1 0 [m] NA NA NA
# time 1 9 1981-09-01 02:00:00 CEST 1 days POSIXct NA
# values
# geometry POINT (56.28161 24.16192),...,POINT (97.80801 36.70394)
# zlev NULL
# time NULL
which I think closes this issue.
edzer
on 7 Oct 2020
Amazing!!! @edzer thanks a lot for sf, stars, ...
meteosimon
on 14 Oct 2020
Related issues
noamross
路
10Comments
michaeldorman
路
4Comments
michaeldorman
路
5Comments
rgzn
路
7Comments
dblodgett-usgs
路
16Comments