forked from r-spatial/stars
-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
267 lines (211 loc) · 10.1 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
---
output: github_document
---
```{r setup, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#",
fig.path = "man/figures/README-"
)
is_online = curl::has_internet()
```
# Spatiotemporal Arrays: Raster and Vector Datacubes
<!-- badges: start -->
[![R-CMD-check](https://github.com/r-spatial/stars/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/r-spatial/stars/actions/workflows/R-CMD-check.yaml)
[![CRAN](https://www.r-pkg.org/badges/version/stars)](https://cran.r-project.org/package=stars)
[![cran checks](https://badges.cranchecks.info/worst/stars.svg)](https://cran.r-project.org/web/checks/check_results_stars.html)
[![Downloads](https://cranlogs.r-pkg.org/badges/stars?color=brightgreen)](https://www.r-pkg.org/pkg/stars)
[![status](https://tinyverse.netlify.app/badge/stars)](https://CRAN.R-project.org/package=stars)
[![Codecov test coverage](https://codecov.io/gh/r-spatial/stars/branch/main/graph/badge.svg)](https://app.codecov.io/gh/r-spatial/stars?branch=main)
<!-- badges: end -->
Spatiotemporal data often comes in the form of dense arrays, with space and time being array dimensions. Examples include
- socio-economic or demographic data,
- environmental variables monitored at fixed stations,
- raster maps,
- time series of satellite images with multiple spectral bands,
- spatial simulations, and
- climate or weather model output.
This R package provides classes and methods for reading,
manipulating, plotting and writing such data cubes, to the extent
that there are proper formats for doing so.
## Raster and vector data cubes
The canonical data cube most of us have in mind is that where two
dimensions represent spatial raster dimensions, and the third time
(or band), as e.g. shown here:
```{r cube1,out.width='50%',echo=FALSE}
knitr::include_graphics("https://raw.githubusercontent.com/r-spatial/stars/main/images/cube1.png")
```
By data cubes however we also consider higher-dimensional cubes
(hypercubes) such as a five-dimensional cube where in addition to
time, spectral band and sensor form dimensions:
```{r cube2,out.width='50%',echo=FALSE}
knitr::include_graphics("https://raw.githubusercontent.com/r-spatial/stars/main/images/cube2.png")
```
or lower-dimensional cubes such as a raster image:
```{r plot1, message=FALSE}
library(dplyr)
library(stars)
tif = system.file("tif/L7_ETMs.tif", package = "stars")
read_stars(tif) |>
slice(index = 1, along = "band") |>
plot()
```
Raster data do not need to be regular and aligned with North/East,
and package `stars` supports besides _regular_ also _rotated_,
_sheared_, _rectilinear_ and _curvilinear_ rasters:
```{r plot2,echo=FALSE}
x = 1:5
y = 1:4
d = st_dimensions(x = x, y = y, .raster = c("x", "y"))
m = matrix(runif(20),5,4)
r1 = st_as_stars(r = m, dimensions = d)
r = attr(d, "raster")
r$affine = c(0.2, -0.2)
attr(d, "raster") = r
r2 = st_as_stars(r = m, dimensions = d)
r = attr(d, "raster")
r$affine = c(0.1, -0.3)
attr(d, "raster") = r
r3 = st_as_stars(r = m, dimensions = d)
x = c(1, 2, 3.5, 5, 6)
y = c(1, 1.5, 3, 3.5)
d = st_dimensions(x = x, y = y, .raster = c("x", "y"))
r4 = st_as_stars(r = m, dimensions = d)
grd = st_make_grid(cellsize = c(10,10), offset = c(-130,10), n = c(8,5), crs = st_crs(4326))
r5 = st_transform(grd, "+proj=laea +lon_0=-70 +lat_0=35")
par(mfrow = c(2,3))
r1 = st_make_grid(cellsize = c(1,1), n = c(5,4), offset = c(0,0))
plot(r1, main = "regular")
plot(st_geometry(st_as_sf(r2)), main = "rotated")
plot(st_geometry(st_as_sf(r3)), main = "sheared")
plot(st_geometry(st_as_sf(r4, as_points = FALSE)), main = "rectilinear")
plot(st_geometry((r5)), main = "curvilinear")
```
Vector data cubes arise when we do not have two regularly discretized spatial dimensions, but a single dimension that points to distinct spatial feature geometries, such as polygons (e.g. denoting administrative regions):
```{r cube3,out.width='50%',echo=FALSE}
knitr::include_graphics("https://raw.githubusercontent.com/r-spatial/stars/main/images/cube3.png")
```
or points (e.g. denoting sensor locations):
```{r cube4,out.width='50%',echo=FALSE}
knitr::include_graphics("https://raw.githubusercontent.com/r-spatial/stars/main/images/cube4.png")
```
NetCDF's CF-convention calls this a [discrete axis](https://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html#discrete-axis).
## NetCDF, GDAL
`stars` provides two functions to read data: `read_ncdf` and
`read_stars`, where the latter reads through GDAL. (In the future,
both will be integrated in `read_stars`.) For reading NetCDF files,
package `RNetCDF` is used, for reading through GDAL, package `sf`
provides the binary linking to GDAL.
For vector and raster operations, `stars` uses as much as possible
the routines available in GDAL and PROJ (e.g. `st_transform`,
`rasterize`, `polygonize`, `warp`). Read more about this in
the vignette on [vector-raster conversions, reprojection,
warping](https://r-spatial.github.io/stars/articles/stars5.html).
## Out-of-memory (on-disk) rasters
Package `stars` provides `stars_proxy` objects (currently only when
read through GDAL), which contain only the dimensions metadata and
pointers to the files on disk. These objects work lazily: reading
and processing data is postponed to the moment that pixels are
really needed (at plot time, or when writing to disk), and is done
at the lowest spatial resolution possible that still fulfills the
resolution of the graphics device. More
details are found in the [stars proxy
vignette](https://r-spatial.github.io/stars/articles/stars2.html).
The following methods are currently available for `stars_proxy` objects:
```{r methods}
methods(class = "stars_proxy")
```
## Raster and vector time series analysis example
In the following, a curvilinear grid with hourly precipitation values of
a hurricane is imported and the first 12 time steps are plotted:
```{r plot3}
prec_file = system.file("nc/test_stageiv_xyt.nc", package = "stars")
(prec = read_stars(gdal_subdatasets(prec_file)[[1]]))
# or: (prec = read_ncdf(prec_file, curvilinear = c("lon", "lat"), ignore_bounds = TRUE))
sf::read_sf(system.file("gpkg/nc.gpkg", package = "sf"), "nc.gpkg") |>
st_transform(st_crs(prec)) -> nc # transform from NAD27 to WGS84
nc_outline = st_union(st_geometry(nc))
plot_hook = function() plot(nc_outline, border = 'red', add = TRUE)
prec |>
slice(index = 1:12, along = "time") |>
plot(downsample = c(3, 3, 1), hook = plot_hook)
```
and next, intersected with with the counties of North Carolina, where
the maximum precipitation intensity was obtained per county, and plotted:
```{r plot4}
a = aggregate(prec, by = nc, FUN = max)
plot(a, max.plot = 23, border = 'grey', lwd = .5)
```
We can integrate over (reduce) time, for instance to find out _when_
the maximum precipitation occurred. The following code finds the time
index, and then the corresponding time value:
```{r plot5, fig.height=3}
index_max = function(x) ifelse(all(is.na(x)), NA, which.max(x))
b = st_apply(a, "geom", index_max)
b |> mutate(when = st_get_dimension_values(a, "time")[b$index_max]) |>
select(when) |>
plot(key.pos = 1, main = "time of maximum precipitation")
```
With package `cubble`, we can make a glyph map to see the magnitude and timings of county maximum precipitation:
```{r plot6, fig.height=3, message=FALSE, warning=FALSE}
library(cubble)
library(ggplot2)
a |> setNames("precip") |>
st_set_dimensions(2, name = "tm") |>
units::drop_units() |>
as_cubble(key = id, index = tm) -> a.cb
a.cb |>
face_temporal() |>
unfold(long, lat) |>
mutate(tm = as.numeric(tm)) |>
ggplot(aes(x_major = long, x_minor = tm, y_major = lat, y_minor = precip)) +
geom_sf(data = nc, inherit.aes = FALSE) +
geom_glyph_box(width = 0.3, height = 0.1) +
geom_glyph(width = 0.3, height = 0.1)
```
## Other packages for data cubes
### [`gdalcubes`](https://github.com/appelmar/gdalcubes)
Package `gdalcubes` can be used to create data cubes (or functions
from them) from image collections, sets of multi-band images with
varying
* spatial resolution
* spatial extent
* coordinate reference systems (e.g., spread over multiple UTM zones)
* observation times
and does this by resampling and/or aggregating over space and/or
time. It reuses GDAL VRT's and gdalwarp for spatial resampling and/or
warping, and handles temporal resampling or aggregation itself.
### [`ncdfgeom`](https://github.com/DOI-USGS/ncdfgeom)
`ncdfgeom` reads and writes vector data cubes from and to netcdf
files in a standards-compliant way.
### [`raster`](https://github.com/rspatial/raster/) and [`terra`](https://github.com/rspatial/terra/)
Packages `raster` and its successor, `terra` are powerful packages
for handling raster maps and stacks of raster maps both in memory
and on disk, but do not address
* non-raster time series,
* multi-attribute rasters time series
* rasters with mixed type attributes (e.g., numeric, logical, factor, POSIXct)
* rectilinear or curvilinear rasters
A list of `stars` commands matching
existing `raster` commands is found in this
[wiki](https://github.com/r-spatial/stars/wiki/How-%60raster%60-functions-map-to-%60stars%60-functions).
A list of translations in the opposite direction (from `stars` to
`raster` or `terra`) still needs to be made.
A comment on the differences between `stars` and `terra` is found
[here](https://github.com/r-spatial/stars/issues/633).
## Other `stars` resources:
* blog posts: [first](https://r-spatial.org/r/2017/11/23/stars1.html),
[second](https://r-spatial.org/r/2018/03/22/stars2.html),
[third](https://r-spatial.org/r/2018/03/23/stars3.html), and [newer blog posts](https://r-spatial.org/)
* [vignettes](https://r-spatial.github.io/stars/articles/)
* the original [R Consortium proposal](https://github.com/r-spatial/stars/blob/main/PROPOSAL.md).
### Acknowledgment
This project has been realized with financial
[support](https://www.r-consortium.org/blog/2017/04/03/q1-2017-isc-grants)
from the
<a href="https://r-consortium.org/all-projects/2017-group-1.html#stars-scalable-spatiotemporal-tidy-arrays-for-r">
<img src="https://r-consortium.org/images/RConsortium_Horizontal_Pantone.webp" width="300">
</a>
<!--
<img src="http://pebesma.staff.ifgi.de/RConsortium_Horizontal_Pantone.png" width="300">
-->