Title: | Spatiotemporal Arrays, Raster and Vector Data Cubes |
---|---|
Description: | Reading, manipulating, writing and plotting spatiotemporal arrays (raster and vector data cubes) in 'R', using 'GDAL' bindings provided by 'sf', and 'NetCDF' bindings by 'ncmeta' and 'RNetCDF'. |
Authors: | Edzer Pebesma [aut, cre] |
Maintainer: | Edzer Pebesma <[email protected]> |
License: | Apache License |
Version: | 0.6-7 |
Built: | 2024-09-02 16:21:44 UTC |
Source: | https://github.com/r-spatial/stars |
evaluate whether cube values are in a given set
## S4 method for signature 'stars' x %in% table
## S4 method for signature 'stars' x %in% table
x |
data cube value |
table |
values of the set |
spatially or temporally aggregate stars object, returning a data cube with lower spatial or temporal resolution
## S3 method for class 'stars' aggregate( x, by, FUN, ..., drop = FALSE, join = st_intersects, as_points = any(st_dimension(by) == 2, na.rm = TRUE), rightmost.closed = FALSE, left.open = FALSE, exact = FALSE )
## S3 method for class 'stars' aggregate( x, by, FUN, ..., drop = FALSE, join = st_intersects, as_points = any(st_dimension(by) == 2, na.rm = TRUE), rightmost.closed = FALSE, left.open = FALSE, exact = FALSE )
x |
object of class |
by |
object of class |
FUN |
aggregation function, such as |
... |
arguments passed on to |
drop |
logical; ignored |
join |
function; function used to find matches of |
as_points |
see st_as_sf: shall raster pixels be taken as points, or small square polygons? |
rightmost.closed |
see findInterval |
left.open |
logical; used for time intervals, see findInterval and cut.POSIXt |
exact |
logical; if |
aggregate, st_interpolate_aw, st_extract, https://github.com/r-spatial/stars/issues/317
# aggregate time dimension in format Date tif = system.file("tif/L7_ETMs.tif", package = "stars") t1 = as.Date("2018-07-31") x = read_stars(c(tif, tif, tif, tif), along = list(time = c(t1, t1+1, t1+2, t1+3)))[,1:30,1:30] st_get_dimension_values(x, "time") x_agg_time = aggregate(x, by = t1 + c(0, 2, 4), FUN = max) # aggregate time dimension in format Date - interval by_t = "2 days" x_agg_time2 = aggregate(x, by = by_t, FUN = max) st_get_dimension_values(x_agg_time2, "time") #TBD: #x_agg_time - x_agg_time2 # aggregate time dimension in format POSIXct x = st_set_dimensions(x, 4, values = as.POSIXct(c("2018-07-31", "2018-08-01", "2018-08-02", "2018-08-03")), names = "time") by_t = as.POSIXct(c("2018-07-31", "2018-08-02")) x_agg_posix = aggregate(x, by = by_t, FUN = max) st_get_dimension_values(x_agg_posix, "time") #TBD: # x_agg_time - x_agg_posix aggregate(x, "2 days", mean) if (require(ncmeta, quietly = TRUE)) { # Spatial aggregation, see https://github.com/r-spatial/stars/issues/299 prec_file = system.file("nc/test_stageiv_xyt.nc", package = "stars") prec = read_ncdf(prec_file, curvilinear = c("lon", "lat")) prec_slice = dplyr::slice(prec, index = 17, along = "time") nc = sf::read_sf(system.file("gpkg/nc.gpkg", package = "sf"), "nc.gpkg") nc = st_transform(nc, st_crs(prec_slice)) agg = aggregate(prec_slice, st_geometry(nc), mean) plot(agg) } # example of using a function for "by": aggregate by month-of-year d = c(10, 10, 150) a = array(rnorm(prod(d)), d) # pure noise times = Sys.Date() + seq(1, 2000, length.out = d[3]) m = as.numeric(format(times, "%m")) signal = rep(sin(m / 12 * pi), each = prod(d[1:2])) # yearly period s = (st_as_stars(a) + signal) %>% st_set_dimensions(3, values = times) f = function(x, format = "%B") { months = format(as.Date(paste0("01-", 1:12, "-1970")), format) factor(format(x, format), levels = months) } agg = aggregate(s, f, mean) plot(agg)
# aggregate time dimension in format Date tif = system.file("tif/L7_ETMs.tif", package = "stars") t1 = as.Date("2018-07-31") x = read_stars(c(tif, tif, tif, tif), along = list(time = c(t1, t1+1, t1+2, t1+3)))[,1:30,1:30] st_get_dimension_values(x, "time") x_agg_time = aggregate(x, by = t1 + c(0, 2, 4), FUN = max) # aggregate time dimension in format Date - interval by_t = "2 days" x_agg_time2 = aggregate(x, by = by_t, FUN = max) st_get_dimension_values(x_agg_time2, "time") #TBD: #x_agg_time - x_agg_time2 # aggregate time dimension in format POSIXct x = st_set_dimensions(x, 4, values = as.POSIXct(c("2018-07-31", "2018-08-01", "2018-08-02", "2018-08-03")), names = "time") by_t = as.POSIXct(c("2018-07-31", "2018-08-02")) x_agg_posix = aggregate(x, by = by_t, FUN = max) st_get_dimension_values(x_agg_posix, "time") #TBD: # x_agg_time - x_agg_posix aggregate(x, "2 days", mean) if (require(ncmeta, quietly = TRUE)) { # Spatial aggregation, see https://github.com/r-spatial/stars/issues/299 prec_file = system.file("nc/test_stageiv_xyt.nc", package = "stars") prec = read_ncdf(prec_file, curvilinear = c("lon", "lat")) prec_slice = dplyr::slice(prec, index = 17, along = "time") nc = sf::read_sf(system.file("gpkg/nc.gpkg", package = "sf"), "nc.gpkg") nc = st_transform(nc, st_crs(prec_slice)) agg = aggregate(prec_slice, st_geometry(nc), mean) plot(agg) } # example of using a function for "by": aggregate by month-of-year d = c(10, 10, 150) a = array(rnorm(prod(d)), d) # pure noise times = Sys.Date() + seq(1, 2000, length.out = d[3]) m = as.numeric(format(times, "%m")) signal = rep(sin(m / 12 * pi), each = prod(d[1:2])) # yearly period s = (st_as_stars(a) + signal) %>% st_set_dimensions(3, values = times) f = function(x, format = "%B") { months = format(as.Date(paste0("01-", 1:12, "-1970")), format) factor(format(x, format), levels = months) } agg = aggregate(s, f, mean) plot(agg)
Coerce stars object into a Raster raster or brick
Coerce stars object into a terra SpatRaster
from |
object to coerce |
If the stars object has more than three dimensions, all dimensions higher than the third will be collapsed into the third dimensions. If the stars object has only an x/y raster but multiple attributes, these are merged first, then put in a raster brick.
If the stars object has more than three dimensions, all dimensions higher than the third will be collapsed into the third dimensions. If the stars object has only an x/y raster but multiple attributes, these are merged first, then put in a SpatRaster.
RasterLayer or RasterBrick
SpatRaster
These are the monthly observational data used for BCSD downscaling. See: <https://gdo-dcp.ucllnl.org/downscaled_cmip_projections/dcpInterface.html#About> for more information." ; "Atmospheric Temperature, Air Temperature Atmosphere, Precipitation, Rain, Maximum Daily Temperature, Minimum Daily Temperature" ;
bcsd_obs
bcsd_obs
An object of class stars_proxy
(inherits from stars
) of dimension 81 x 33 x 12.
combine multiple stars objects, or combine multiple attributes in a single stars object into a single array
## S3 method for class 'stars_proxy' c( ..., along = NA_integer_, along_crs = FALSE, try_hard = FALSE, nms = names(list(...)), tolerance = sqrt(.Machine$double.eps) ) ## S3 method for class 'stars' c( ..., along = NA_integer_, try_hard = FALSE, nms = names(list(...)), tolerance = sqrt(.Machine$double.eps) )
## S3 method for class 'stars_proxy' c( ..., along = NA_integer_, along_crs = FALSE, try_hard = FALSE, nms = names(list(...)), tolerance = sqrt(.Machine$double.eps) ) ## S3 method for class 'stars' c( ..., along = NA_integer_, try_hard = FALSE, nms = names(list(...)), tolerance = sqrt(.Machine$double.eps) )
... |
object(s) of class |
along |
integer; see read_stars |
along_crs |
logical; if |
try_hard |
logical; if |
nms |
character; vector with array names |
tolerance |
numeric; values used in all.equal to compare dimension values combine those that dimensions matching to the first array |
An error is raised when attempting to combine arrays with different
measurement units into a single array. If this was intentded, drop_units
can be used to remove units of a stars
object before merging.
a single stars
object with merged (binded) arrays.
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) (new = c(x, x)) c(new) # collapses two arrays into one with an additional dimension c(x, x, along = 3)
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) (new = c(x, x)) c(new) # collapses two arrays into one with an additional dimension c(x, x, along = 3)
plot contours of a stars object
## S3 method for class 'stars' contour(x, ...)
## S3 method for class 'stars' contour(x, ...)
x |
object of class |
... |
other parameters passed on to contour |
this uses the R internal contour algorithm, which (by default) plots contours; st_contour uses the GDAL contour algorithm that returns contours as simple features.
d = st_dimensions(x = 1:ncol(volcano), y = 1:nrow(volcano)) r = st_as_stars(t(volcano)) r = st_set_dimensions(r, 1, offset = 0, delta = 1) r = st_set_dimensions(r, 2, offset = 0, delta = -1) plot(r, reset = FALSE) contour(r, add = TRUE)
d = st_dimensions(x = 1:ncol(volcano), y = 1:nrow(volcano)) r = st_as_stars(t(volcano)) r = st_set_dimensions(r, 1, offset = 0, delta = 1) r = st_set_dimensions(r, 2, offset = 0, delta = -1) plot(r, reset = FALSE) contour(r, add = TRUE)
cut methods for stars objects
## S3 method for class 'array' cut(x, breaks, ...) ## S3 method for class 'matrix' cut(x, breaks, ...) ## S3 method for class 'stars' cut(x, breaks, ...)
## S3 method for class 'array' cut(x, breaks, ...) ## S3 method for class 'matrix' cut(x, breaks, ...) ## S3 method for class 'stars' cut(x, breaks, ...)
x |
see cut |
breaks |
see cut |
... |
see cut |
R's factor
only works for vectors, not for arrays or matrices. This is a work-around (or hack?) to keep the factor levels generated by cut
and use them in plots.
an array or matrix with a levels
attribute; see details
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) cut(x, c(0, 50, 100, 255)) cut(x[,,,1], c(0, 50, 100, 255)) plot(cut(x[,,,1], c(0, 50, 100, 255))) tif = system.file("tif/L7_ETMs.tif", package = "stars") x1 = read_stars(tif) (x1_cut = cut(x1, breaks = c(0, 50, 100, Inf))) # shows factor in summary plot(x1_cut[,,,c(3,6)]) # propagates through [ and plot
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) cut(x, c(0, 50, 100, 255)) cut(x[,,,1], c(0, 50, 100, 255)) plot(cut(x[,,,1], c(0, 50, 100, 255))) tif = system.file("tif/L7_ETMs.tif", package = "stars") x1 = read_stars(tif) (x1_cut = cut(x1, breaks = c(0, 50, 100, Inf))) # shows factor in summary plot(x1_cut[,,,c(3,6)]) # propagates through [ and plot
dplyr verbs for stars objects; package dplyr needs to be loaded before these methods can be used for stars objects.
filter.stars(.data, ...) filter.stars_proxy(.data, ...) mutate.stars(.data, ...) mutate.stars_proxy(.data, ...) transmute.stars(.data, ...) transmute.stars_proxy(.data, ...) select.stars(.data, ...) select.stars_proxy(.data, ...) rename.stars(.data, ...) rename.stars_proxy(.data, ...) pull.stars(.data, var = -1) pull.stars_proxy(.data, ...) as.tbl_cube.stars(x, ...) slice.stars(.data, along, index, ..., drop = length(index) == 1) slice.stars_proxy(.data, along, index, ...) replace_na.stars(data, replace, ...) replace_na.stars_proxy(data, ...)
filter.stars(.data, ...) filter.stars_proxy(.data, ...) mutate.stars(.data, ...) mutate.stars_proxy(.data, ...) transmute.stars(.data, ...) transmute.stars_proxy(.data, ...) select.stars(.data, ...) select.stars_proxy(.data, ...) rename.stars(.data, ...) rename.stars_proxy(.data, ...) pull.stars(.data, var = -1) pull.stars_proxy(.data, ...) as.tbl_cube.stars(x, ...) slice.stars(.data, along, index, ..., drop = length(index) == 1) slice.stars_proxy(.data, along, index, ...) replace_na.stars(data, replace, ...) replace_na.stars_proxy(data, ...)
.data |
object of class |
... |
see filter |
var |
see pull |
x |
object of class |
along |
name or index of dimension to which the slice should be applied |
index |
integer value(s) for this index |
drop |
logical; drop dimensions that only have a single index? |
data |
data set to work on |
replace |
see replace_na: list with variable=value pairs, where value is the replacement value for NA's |
tif = system.file("tif/L7_ETMs.tif", package = "stars") x1 = read_stars(tif) if (require(dplyr, quietly = TRUE)) { x1 %>% slice("band", 2:3) x1 %>% slice("x", 50:100) }
tif = system.file("tif/L7_ETMs.tif", package = "stars") x1 = read_stars(tif) if (require(dplyr, quietly = TRUE)) { x1 %>% slice("band", 2:3) x1 %>% slice("x", 50:100) }
expand the dimension values into a list
expand_dimensions(x, ...) ## S3 method for class 'dimensions' expand_dimensions(x, ..., max = FALSE, center = NA)
expand_dimensions(x, ...) ## S3 method for class 'dimensions' expand_dimensions(x, ..., max = FALSE, center = NA)
x |
object of class 'stars' or 'dimensions' |
... |
ignored |
max |
logical; if 'TRUE' return the max (end) values of the dimensions intervals |
center |
logical; if 'TRUE' return the center values of intervals, otherwise return offset (start) of intervals; if 'NA' (default) return centers for x/y dimensions, offsets for all others |
ggplot geom for stars objects
geom_stars( mapping = NULL, data = NULL, ..., downsample = 0, sf = FALSE, na.action = na.pass ) theme_stars(...)
geom_stars( mapping = NULL, data = NULL, ..., downsample = 0, sf = FALSE, na.action = na.pass ) theme_stars(...)
mapping |
see geom_raster |
data |
see geom_raster |
... |
see geom_raster |
downsample |
downsampling rate: e.g. 3 keeps rows and cols 1, 4, 7, 10 etc.; a value of 0 does not downsample; can be specified for each dimension, e.g. |
sf |
logical; if |
na.action |
function; if |
geom_stars
returns (a call to) either geom_raster, geom_tile, or geom_sf, depending on the raster or vector geometry; for the first to, an aes call is constructed with the raster dimension names and the first array as fill variable. Further calls to coord_equal and facet_wrap are needed to control aspect ratio and the layers to be plotted; see examples. If a stars
array contains hex color values, and no fill
parameter is given, the color values are used as fill color; see the example below.
If visual artefacts occur (Moiré-Effekt), then see the details section of plot.stars
system.file("tif/L7_ETMs.tif", package = "stars") %>% read_stars() -> x if (require(ggplot2, quietly = TRUE)) { ggplot() + geom_stars(data = x) + coord_equal() + facet_wrap(~band) + theme_void() + scale_x_discrete(expand=c(0,0))+ scale_y_discrete(expand=c(0,0)) # plot rgb composite: st_as_stars(L7_ETMs)[,,,1:3] |> st_rgb() -> x # x contains colors as pixel values ggplot() + geom_stars(data = x) }
system.file("tif/L7_ETMs.tif", package = "stars") %>% read_stars() -> x if (require(ggplot2, quietly = TRUE)) { ggplot() + geom_stars(data = x) + coord_equal() + facet_wrap(~band) + theme_void() + scale_x_discrete(expand=c(0,0))+ scale_y_discrete(expand=c(0,0)) # plot rgb composite: st_as_stars(L7_ETMs)[,,,1:3] |> st_rgb() -> x # x contains colors as pixel values ggplot() + geom_stars(data = x) }
Probably containing the six 30 m bands:
Band 1 Visible (0.45 - 0.52 µm) 30 m
Band 2 Visible (0.52 - 0.60 µm) 30 m
Band 3 Visible (0.63 - 0.69 µm) 30 m
Band 4 Near-Infrared (0.77 - 0.90 µm) 30 m
Band 5 Short-wave Infrared (1.55 - 1.75 µm) 30 m
Band 7 Mid-Infrared (2.08 - 2.35 µm) 30 m
L7_ETMs
L7_ETMs
An object of class stars_proxy
(inherits from stars
) of dimension 349 x 352 x 6.
create an intervals object, assuming left-closed and right-open intervals
make_intervals(start, end)
make_intervals(start, end)
start |
vector with start values, or 2-column matrix with start and end values in column 1 and 2, respectively |
end |
vector with end values |
Read or write data using GDAL's multidimensional array API
read_mdim( filename, variable = character(0), ..., options = character(0), raster = NULL, offset = integer(0), count = integer(0), step = integer(0), proxy = FALSE, debug = FALSE, bounds = TRUE, curvilinear = NA ) write_mdim( x, filename, driver = detect.driver(filename), ..., root_group_options = character(0), options = character(0), as_float = TRUE )
read_mdim( filename, variable = character(0), ..., options = character(0), raster = NULL, offset = integer(0), count = integer(0), step = integer(0), proxy = FALSE, debug = FALSE, bounds = TRUE, curvilinear = NA ) write_mdim( x, filename, driver = detect.driver(filename), ..., root_group_options = character(0), options = character(0), as_float = TRUE )
filename |
name of the source or destination file or data source |
variable |
name of the array to be read; if '"?"', a list of array names is returned, with group name as list element names. |
... |
ignored |
options |
character; driver specific options regarding the opening (read_mdim) or creation (write_mdim) of the dataset |
raster |
names of the raster variables (default: first two dimensions) |
offset |
integer; zero-based offset for each dimension (pixels) of sub-array to read, defaults to 0 for each dimension(requires sf >= 1.0-9) |
count |
integer; size for each dimension (pixels) of sub-array to read (default: read all); a value of NA will read the corresponding dimension entirely; counts are relative to the step size (requires sf >= 1.0-9) |
step |
integer; step size for each dimension (pixels) of sub-array to read; defaults to 1 for each dimension (requires sf >= 1.0-9) |
proxy |
logical; return proxy object? |
debug |
logical; print debug info? |
bounds |
logical or character: if |
curvilinear |
control reading curvilinear (geolocation) coordinate arrays; if |
x |
stars object |
driver |
character; driver name |
root_group_options |
character; driver specific options regarding the creation of the root group |
as_float |
logical; if |
it is assumed that the first two dimensions are easting and northing
gdal_utils, in particular util mdiminfo
to query properties of a file or data source containing arrays
set.seed(135) m = matrix(runif(10), 2, 5) names(dim(m)) = c("stations", "time") times = as.Date("2022-05-01") + 1:5 pts = st_as_sfc(c("POINT(0 1)", "POINT(3 5)")) s = st_as_stars(list(Precipitation = m)) |> st_set_dimensions(1, values = pts) |> st_set_dimensions(2, values = times) nc = tempfile(fileext=".nc") if (compareVersion(sf_extSoftVersion()["GDAL"], "3.4.0") > -1) { write_mdim(s, nc) # try ncdump on the generated file print(read_mdim(nc)) }
set.seed(135) m = matrix(runif(10), 2, 5) names(dim(m)) = c("stations", "time") times = as.Date("2022-05-01") + 1:5 pts = st_as_sfc(c("POINT(0 1)", "POINT(3 5)")) s = st_as_stars(list(Precipitation = m)) |> st_set_dimensions(1, values = pts) |> st_set_dimensions(2, values = times) nc = tempfile(fileext=".nc") if (compareVersion(sf_extSoftVersion()["GDAL"], "3.4.0") > -1) { write_mdim(s, nc) # try ncdump on the generated file print(read_mdim(nc)) }
merge attributes into a dimension, or split a dimension over attributes
## S3 method for class 'stars' split(x, f = length(dim(x)), drop = TRUE, ...) ## S3 method for class 'stars' merge(x, y, ..., name = "attributes")
## S3 method for class 'stars' split(x, f = length(dim(x)), drop = TRUE, ...) ## S3 method for class 'stars' merge(x, y, ..., name = "attributes")
x |
object of class |
f |
the name or index of the dimension to split; by default the last dimension |
drop |
ignored |
... |
if defined, the first unnamed argument is used for dimension values, if not defined, attribute names are used for dimension values |
y |
needs to be missing |
name |
name for the new dimension |
split.stars works on the first attribute, and will give an error when more than one attribute is present
merge merges attributes of a stars object into a new dimension; split splits a dimension over attributes
Ops functions for stars objects, including comparison, product and divide, add, subtract
## S3 method for class 'stars' Ops(e1, e2) ## S3 method for class 'stars' Math(x, ...) ## S3 method for class 'stars_proxy' Ops(e1, e2) ## S3 method for class 'stars_proxy' Math(x, ...)
## S3 method for class 'stars' Ops(e1, e2) ## S3 method for class 'stars' Math(x, ...) ## S3 method for class 'stars_proxy' Ops(e1, e2) ## S3 method for class 'stars_proxy' Math(x, ...)
e1 |
object of class |
e2 |
object of class |
x |
object of class stars |
... |
parameters passed on to the Math functions |
if e1
or e2
is is a numeric vector, or e2
has less or smaller dimensions than e1
, then e2
is recycled
such that it fits e1
, using usual R array recycling rules. The user
needs to make sure this is sensible; it may be needed to use aperm
to permutate dimensions first.
object of class stars
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) x * x x / x x + x x + 10 all.equal(x * 10, 10 * x) tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) a = sqrt(x) b = log(x, base = 10)
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) x * x x / x x + x x + 10 all.equal(x * 10, 10 * x) tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) a = sqrt(x) b = log(x, base = 10)
plot stars object, with subplots for each level of first non-spatial dimension, and customization of legend key
## S3 method for class 'nc_proxy' plot(x, y, ..., downsample = get_downsample(dim(x)), max_times = 16) ## S3 method for class 'stars' plot( x, y, ..., join_zlim = TRUE, main = make_label(x, 1), axes = FALSE, downsample = TRUE, nbreaks = 11, breaks = "quantile", col = grey(1:(nbreaks - 1)/nbreaks), key.pos = get_key_pos(x, ...), key.width = kw_dflt(x, key.pos), key.length = 0.618, key.lab = main, reset = TRUE, box_col = NA, center_time = FALSE, hook = NULL, mfrow = NULL, compact = TRUE ) ## S3 method for class 'stars' image( x, ..., band = 1, attr = 1, asp = NULL, rgb = NULL, maxColorValue = ifelse(inherits(rgb, "data.frame"), 255, max(x[[attr]], na.rm = TRUE)), xlab = if (!axes) "" else names(d)[1], ylab = if (!axes) "" else names(d)[2], xlim = st_bbox(extent)$xlim, ylim = st_bbox(extent)$ylim, text_values = FALSE, text_color = "black", axes = FALSE, interpolate = FALSE, as_points = FALSE, key.pos = NULL, logz = FALSE, key.width = kw_dflt(x, key.pos), key.length = 0.618, add.geom = NULL, border = NA, useRaster = isTRUE(dev.capabilities()$rasterImage == "yes"), extent = x ) ## S3 method for class 'stars_proxy' plot(x, y, ..., downsample = get_downsample(dim(x)))
## S3 method for class 'nc_proxy' plot(x, y, ..., downsample = get_downsample(dim(x)), max_times = 16) ## S3 method for class 'stars' plot( x, y, ..., join_zlim = TRUE, main = make_label(x, 1), axes = FALSE, downsample = TRUE, nbreaks = 11, breaks = "quantile", col = grey(1:(nbreaks - 1)/nbreaks), key.pos = get_key_pos(x, ...), key.width = kw_dflt(x, key.pos), key.length = 0.618, key.lab = main, reset = TRUE, box_col = NA, center_time = FALSE, hook = NULL, mfrow = NULL, compact = TRUE ) ## S3 method for class 'stars' image( x, ..., band = 1, attr = 1, asp = NULL, rgb = NULL, maxColorValue = ifelse(inherits(rgb, "data.frame"), 255, max(x[[attr]], na.rm = TRUE)), xlab = if (!axes) "" else names(d)[1], ylab = if (!axes) "" else names(d)[2], xlim = st_bbox(extent)$xlim, ylim = st_bbox(extent)$ylim, text_values = FALSE, text_color = "black", axes = FALSE, interpolate = FALSE, as_points = FALSE, key.pos = NULL, logz = FALSE, key.width = kw_dflt(x, key.pos), key.length = 0.618, add.geom = NULL, border = NA, useRaster = isTRUE(dev.capabilities()$rasterImage == "yes"), extent = x ) ## S3 method for class 'stars_proxy' plot(x, y, ..., downsample = get_downsample(dim(x)))
x |
object of class |
y |
ignored |
... |
further arguments: for |
downsample |
logical or numeric; if |
max_times |
integer; maximum number of time steps to attempt to plot. |
join_zlim |
logical; if |
main |
character; subplot title prefix; use |
axes |
logical; should axes and box be added to the plot? |
nbreaks |
number of color breaks; should be one more than number of colors. If missing and |
breaks |
numeric vector with actual color breaks, or a |
col |
colors to use for grid cells, or color palette function |
key.pos |
numeric; side to plot a color key: 1 bottom, 2 left, 3 top, 4 right; set to |
key.width |
amount of space reserved for width of the key (labels); relative or absolute (using lcm) |
key.length |
amount of space reserved for length of the key (labels); relative or absolute (using lcm) |
key.lab |
character; label for color key in case of multiple subplots, use |
reset |
logical; if |
box_col |
color for box around sub-plots; use |
center_time |
logical; if |
hook |
NULL or function; hook function that will be called on every sub-plot; see examples. |
mfrow |
length-2 integer vector with nrows, ncolumns of a composite plot, to override the default layout |
compact |
logical; place facets compactly (TRUE), or spread over the plotting device area? |
band |
integer; which band (dimension) to plot |
attr |
integer; which attribute to plot |
asp |
numeric; aspect ratio of image |
rgb |
integer; specify three bands to form an rgb composite. Experimental: rgb color table; see Details. |
maxColorValue |
numeric; passed on to rgb |
xlab |
character; x axis label |
ylab |
character; y axis label |
xlim |
x axis limits |
ylim |
y axis limits |
text_values |
logical; print values as text on image? |
text_color |
character; color for printed text values |
interpolate |
logical; when using rasterImage (rgb), should pixels be interpolated? |
as_points |
logical; for curvilinear or sheared grids: parameter passed on to st_as_sf, determining whether raster cells will be plotted as symbols (fast, approximate) or small polygons (slow, exact) |
logz |
logical; if |
add.geom |
object of class |
border |
color used for cell borders (only in case |
useRaster |
logical; use the rasterImage capabilities of the graphics device? |
extent |
object which has a |
when plotting a subsetted stars_proxy
object, the default value for argument downsample
will not be computed correctly, and has to be set manually.
Downsampling: a value for downsample
of 0: no downsampling, 1: after every dimension value (pixel/line/band), one value is skipped (half of the original resolution), 2: after every dimension value, 2 values are skipped (one third of the original resolution), etc. If downsample
is TRUE
or a length 1 numeric vector, downsampling is only applied to the raster [x] and [y] dimensions.
To remove unused classes in a categorical raster, use the droplevels function.
When bitmaps show visual artefacts (Moiré effects), make sure that device png is used rather than ragg::agg_png
as the latter uses antialiasing for filled polygons which causes this; see also https://github.com/r-spatial/stars/issues/573 .
use of an rgb color table is experimental; see https://github.com/r-spatial/mapview/issues/208
when plotting a subsetted stars_proxy
object, the default value for argument downsample
will not be computed correctly, and has to be set manually.
st_bbox(L7_ETMs) |> st_as_sfc() |> st_centroid() |> st_coordinates() -> pt hook1 = function() { text(pt[,"X"], pt[,"Y"], "foo", col = 'orange', cex = 2) } plot(L7_ETMs, hook = hook1) x = st_set_dimensions(L7_ETMs, 3, paste0("B_", 1:6)) hook2 = function(..., row, col, nr, nrow, ncol, value, bbox) { str = paste0("row ", row, "/", nrow, ", col ", col, "/", ncol, "\nnr: ", nr, " value: ", value) bbox |> st_as_sfc() |> st_centroid() |> st_coordinates() -> pt text(pt[,"X"], pt[,"Y"], str, col = 'red', cex = 2) } plot(x, hook = hook2, col = grey(c(.2,.25,.3,.35))) if (isTRUE(dev.capabilities()$rasterImage == "yes")) { lc = read_stars(system.file("tif/lc.tif", package = "stars")) levels(lc[[1]]) = abbreviate(levels(lc[[1]]), 6) # so it's not only legend plot(lc, key.pos=4) } tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) image(x, col = grey((3:9)/10)) if (isTRUE(dev.capabilities()$rasterImage == "yes")) { image(x, rgb = c(1,3,5)) # false color composite }
st_bbox(L7_ETMs) |> st_as_sfc() |> st_centroid() |> st_coordinates() -> pt hook1 = function() { text(pt[,"X"], pt[,"Y"], "foo", col = 'orange', cex = 2) } plot(L7_ETMs, hook = hook1) x = st_set_dimensions(L7_ETMs, 3, paste0("B_", 1:6)) hook2 = function(..., row, col, nr, nrow, ncol, value, bbox) { str = paste0("row ", row, "/", nrow, ", col ", col, "/", ncol, "\nnr: ", nr, " value: ", value) bbox |> st_as_sfc() |> st_centroid() |> st_coordinates() -> pt text(pt[,"X"], pt[,"Y"], str, col = 'red', cex = 2) } plot(x, hook = hook2, col = grey(c(.2,.25,.3,.35))) if (isTRUE(dev.capabilities()$rasterImage == "yes")) { lc = read_stars(system.file("tif/lc.tif", package = "stars")) levels(lc[[1]]) = abbreviate(levels(lc[[1]]), 6) # so it's not only legend plot(lc, key.pos=4) } tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) image(x, col = grey((3:9)/10)) if (isTRUE(dev.capabilities()$rasterImage == "yes")) { image(x, rgb = c(1,3,5)) # false color composite }
Compute principle components of stars object
## S3 method for class 'stars_proxy' prcomp(x, ..., downsample = 0) ## S3 method for class 'stars' prcomp(x, ..., quiet = FALSE)
## S3 method for class 'stars_proxy' prcomp(x, ..., downsample = 0) ## S3 method for class 'stars' prcomp(x, ..., quiet = FALSE)
x |
object of class 'stars' or 'stars_proxy' |
... |
see prcomp |
downsample |
see st_as_stars |
quiet |
logical; if 'TRUE', suppress message that PCs will be computed on last dimension; see details |
if 'x' has only one attribute, principle components will be computed in the space of the last dimension of 'x' to predict PC scores into a 'stars' object, use predict.stars; see example below
object of class 'prcomp', see prcomp
l7 = split(st_as_stars(L7_ETMs), 3) # use bands as features l7 |> prcomp() |> plot() l7 |> prcomp() |> predict(l7, model = _) |> merge() |> plot()
l7 = split(st_as_stars(L7_ETMs), 3) # use bands as features l7 |> prcomp() |> plot() l7 |> prcomp() |> predict(l7, model = _) |> merge() |> plot()
Predict values, given a model object, for a stars or stars_proxy object
## S3 method for class 'stars_proxy' predict(object, model, ...) ## S3 method for class 'stars' predict(object, model, ..., drop_dimensions = FALSE)
## S3 method for class 'stars_proxy' predict(object, model, ...) ## S3 method for class 'stars' predict(object, model, ..., drop_dimensions = FALSE)
object |
object of class 'stars' |
model |
model object of a class that has a predict method; check with 'methods(class = class(object))' |
... |
arguments passed on to this predict method |
drop_dimensions |
logical; if 'TRUE', remove dimensions (coordinates etc) from 'data.frame' with predictors |
separate predictors in object need to be separate attributes in object; in case they are e.g. in a band dimension, use 'split(object)'
print stars or dimensions object
## S3 method for class 'dimensions' as.data.frame( x, ..., digits = max(3, getOption("digits") - 3), usetz = TRUE, stars_crs = getOption("stars.crs") %||% 28, all = FALSE ) ## S3 method for class 'dimensions' print(x, ...) ## S3 method for class 'stars' print(x, ..., n = 1e+05, abbrev = 30)
## S3 method for class 'dimensions' as.data.frame( x, ..., digits = max(3, getOption("digits") - 3), usetz = TRUE, stars_crs = getOption("stars.crs") %||% 28, all = FALSE ) ## S3 method for class 'dimensions' print(x, ...) ## S3 method for class 'stars' print(x, ..., n = 1e+05, abbrev = 30)
x |
object of class stars or of class dimensions |
... |
passed on to |
digits |
number of digits to print numbers |
usetz |
logical; used to format |
stars_crs |
maximum width of string for CRS objects |
all |
logical; if |
n |
when prod(dim(x)) > 10 * n, the first n cells are used for attribute summary statistics |
abbrev |
number of characters to abbreviate attribute names to |
Read data from a file (or source) using the NetCDF library directly.
read_ncdf( .x, ..., var = NULL, ncsub = NULL, curvilinear = character(0), eps = sqrt(.Machine$double.eps), ignore_bounds = FALSE, make_time = TRUE, make_units = TRUE, proxy = NULL, downsample = 0 )
read_ncdf( .x, ..., var = NULL, ncsub = NULL, curvilinear = character(0), eps = sqrt(.Machine$double.eps), ignore_bounds = FALSE, make_time = TRUE, make_units = TRUE, proxy = NULL, downsample = 0 )
.x |
NetCDF file or source as a character vector or an nc_proxy object. |
... |
ignored |
var |
variable name or names (they must be on matching grids) |
ncsub |
matrix of start, count columns (see Details) |
curvilinear |
length two character named vector with names of variables holding longitude and latitude values for all raster cells. 'stars' attempts to figure out appropriate curvilinear coordinates if they are not supplied. |
eps |
numeric; dimension value increases are considered identical when they differ less than |
ignore_bounds |
logical; should bounds values for dimensions, if present, be ignored? |
make_time |
if |
make_units |
if |
proxy |
logical; if |
downsample |
integer; number of cells to omit between samples along each dimension.
e.g. |
The following logic is applied to coordinates. If any coordinate axes have regularly spaced coordinate variables they are reduced to the offset/delta form with 'affine = c(0, 0)', otherwise the values of the coordinates are stored and used to define a rectilinear grid.
If the data has two or more dimensions and the first two are regular they are nominated as the 'raster' for plotting.
If the curvilinear
argument is used it specifies the 2D arrays
containing coordinate values for the first two dimensions of the data read. It is currently
assumed that the coordinates are 2D and that they relate to the first two dimensions in
that order.
If var
is not set the first set of variables on a shared grid is used.
start
and count
columns of ncsub must correspond to the variable dimension (nrows)
and be valid index using var.get.nc
convention (start is 1-based). If the count value
is NA
then all steps are included. Axis order must match that of the variable/s being read.
f <- system.file("nc/reduced.nc", package = "stars") if (require(ncmeta, quietly = TRUE)) { read_ncdf(f) read_ncdf(f, var = c("anom")) read_ncdf(f, ncsub = cbind(start = c(1, 1, 1, 1), count = c(10, 12, 1, 1))) } if (require(ncmeta, quietly = TRUE)) { #' precipitation data in a curvilinear NetCDF prec_file = system.file("nc/test_stageiv_xyt.nc", package = "stars") prec = read_ncdf(prec_file, curvilinear = c("lon", "lat"), ignore_bounds = TRUE) } ##plot(prec) ## gives error about unique breaks ## remove NAs, zeros, and give a large number ## of breaks (used for validating in detail) qu_0_omit = function(x, ..., n = 22) { x = units::drop_units(na.omit(x)) c(0, quantile(x[x > 0], seq(0, 1, length.out = n))) } if (require(dplyr, quietly = TRUE)) { prec_slice = slice(prec, index = 17, along = "time") plot(prec_slice, border = NA, breaks = qu_0_omit(prec_slice[[1]]), reset = FALSE) nc = sf::read_sf(system.file("gpkg/nc.gpkg", package = "sf"), "nc.gpkg") plot(st_geometry(nc), add = TRUE, reset = FALSE, col = NA) }
f <- system.file("nc/reduced.nc", package = "stars") if (require(ncmeta, quietly = TRUE)) { read_ncdf(f) read_ncdf(f, var = c("anom")) read_ncdf(f, ncsub = cbind(start = c(1, 1, 1, 1), count = c(10, 12, 1, 1))) } if (require(ncmeta, quietly = TRUE)) { #' precipitation data in a curvilinear NetCDF prec_file = system.file("nc/test_stageiv_xyt.nc", package = "stars") prec = read_ncdf(prec_file, curvilinear = c("lon", "lat"), ignore_bounds = TRUE) } ##plot(prec) ## gives error about unique breaks ## remove NAs, zeros, and give a large number ## of breaks (used for validating in detail) qu_0_omit = function(x, ..., n = 22) { x = units::drop_units(na.omit(x)) c(0, quantile(x[x > 0], seq(0, 1, length.out = n))) } if (require(dplyr, quietly = TRUE)) { prec_slice = slice(prec, index = 17, along = "time") plot(prec_slice, border = NA, breaks = qu_0_omit(prec_slice[[1]]), reset = FALSE) nc = sf::read_sf(system.file("gpkg/nc.gpkg", package = "sf"), "nc.gpkg") plot(st_geometry(nc), add = TRUE, reset = FALSE, col = NA) }
read raster/array dataset from file or connection
read_stars( .x, sub = TRUE, ..., options = character(0), driver = character(0), quiet = FALSE, NA_value = NA_real_, along = NA_integer_, RasterIO = list(), proxy = getOption("stars.n_proxy") %||% 1e+08, curvilinear = character(0), normalize_path = TRUE, RAT = character(0), tolerance = 1e-10, exclude = "", shorten = TRUE )
read_stars( .x, sub = TRUE, ..., options = character(0), driver = character(0), quiet = FALSE, NA_value = NA_real_, along = NA_integer_, RasterIO = list(), proxy = getOption("stars.n_proxy") %||% 1e+08, curvilinear = character(0), normalize_path = TRUE, RAT = character(0), tolerance = 1e-10, exclude = "", shorten = TRUE )
.x |
character vector with name(s) of file(s) or data source(s) to be read, or a function that returns such a vector |
sub |
character, integer or logical; name, index or indicator of sub-dataset(s) to be read |
... |
passed on to st_as_stars if |
options |
character; opening options |
driver |
character; driver to use for opening file. To override fixing for subdatasets and autodetect them as well, use |
quiet |
logical; print progress output? |
NA_value |
numeric value to be used for conversion into NA values; by default this is read from the input file |
along |
length-one character or integer, or list; determines how several arrays are combined, see Details. |
RasterIO |
list with named parameters for GDAL's RasterIO, to further control the extent, resolution and bands to be read from the data source; see details. |
proxy |
logical; if |
curvilinear |
length two character vector with names of subdatasets holding longitude and latitude values for all raster cells, or named length 2 list holding longitude and latitude matrices; the names of this list should correspond to raster dimensions referred to |
normalize_path |
logical; if |
RAT |
character; raster attribute table column name to use as factor levels |
tolerance |
numeric; passed on to all.equal for comparing dimension parameters. |
exclude |
character; vector with category value(s) to exclude |
shorten |
logical or character; if |
In case .x
contains multiple files, they will all be read and combined with c.stars. Along which dimension, or how should objects be merged? If along
is set to NA
it will merge arrays as new attributes if all objects have identical dimensions, or else try to merge along time if a dimension called time
indicates different time stamps. A single name (or positive value) for along
will merge along that dimension, or create a new one if it does not already exist. If the arrays should be arranged along one of more dimensions with values (e.g. time stamps), a named list can passed to along
to specify them; see example.
RasterIO
is a list with zero or more of the following named arguments:
nXOff
, nYOff
(both 1-based: the first row/col has offset value 1),
nXSize
, nYSize
, nBufXSize
, nBufYSize
, bands
, resample
.
See https://gdal.org/doxygen/classGDALDataset.html for their meaning;
bands
is an integer vector containing the band numbers to be read (1-based: first band is 1).
Note that if nBufXSize
or nBufYSize
are specified for downsampling an image,
resulting in an adjusted geotransform. resample
reflects the resampling method and
has to be one of: "nearest_neighbour" (the default),
"bilinear", "cubic", "cubic_spline", "lanczos", "average", "mode", or "Gauss".
Data that are read into memory (proxy=FALSE
) are read into a numeric (double) array, except for categorical variables which are read into an numeric (integer) array of class factor
.
object of class stars
tif = system.file("tif/L7_ETMs.tif", package = "stars") (x1 = read_stars(tif)) (x2 = read_stars(c(tif, tif))) (x3 = read_stars(c(tif, tif), along = "band")) (x4 = read_stars(c(tif, tif), along = "new_dimensions")) # create 4-dimensional array x1o = read_stars(tif, options = "OVERVIEW_LEVEL=1") t1 = as.Date("2018-07-31") # along is a named list indicating two dimensions: read_stars(c(tif, tif, tif, tif), along = list(foo = c("bar1", "bar2"), time = c(t1, t1+2))) m = matrix(1:120, nrow = 12, ncol = 10) dim(m) = c(x = 10, y = 12) # named dim st = st_as_stars(m) attr(st, "dimensions")$y$delta = -1 attr(st, "dimensions")$y$offset = 12 st tmp = tempfile(fileext = ".tif") write_stars(st, tmp) (red <- read_stars(tmp)) read_stars(tmp, RasterIO = list(nXOff = 1, nYOff = 1, nXSize = 10, nYSize = 12, nBufXSize = 2, nBufYSize = 2))[[1]] (red <- read_stars(tmp, RasterIO = list(nXOff = 1, nYOff = 1, nXSize = 10, nYSize = 12, nBufXSize = 2, nBufYSize = 2))) red[[1]] # cell values of subsample grid: ## Not run: plot(st, reset = FALSE, axes = TRUE, ylim = c(-.1,12.1), xlim = c(-.1,10.1), main = "nBufXSize & nBufYSize demo", text_values = TRUE) plot(st_as_sfc(red, as_points = TRUE), add = TRUE, col = 'red', pch = 16) plot(st_as_sfc(st_as_stars(st), as_points = FALSE), add = TRUE, border = 'grey') plot(st_as_sfc(red, as_points = FALSE), add = TRUE, border = 'green', lwd = 2) ## End(Not run) file.remove(tmp)
tif = system.file("tif/L7_ETMs.tif", package = "stars") (x1 = read_stars(tif)) (x2 = read_stars(c(tif, tif))) (x3 = read_stars(c(tif, tif), along = "band")) (x4 = read_stars(c(tif, tif), along = "new_dimensions")) # create 4-dimensional array x1o = read_stars(tif, options = "OVERVIEW_LEVEL=1") t1 = as.Date("2018-07-31") # along is a named list indicating two dimensions: read_stars(c(tif, tif, tif, tif), along = list(foo = c("bar1", "bar2"), time = c(t1, t1+2))) m = matrix(1:120, nrow = 12, ncol = 10) dim(m) = c(x = 10, y = 12) # named dim st = st_as_stars(m) attr(st, "dimensions")$y$delta = -1 attr(st, "dimensions")$y$offset = 12 st tmp = tempfile(fileext = ".tif") write_stars(st, tmp) (red <- read_stars(tmp)) read_stars(tmp, RasterIO = list(nXOff = 1, nYOff = 1, nXSize = 10, nYSize = 12, nBufXSize = 2, nBufYSize = 2))[[1]] (red <- read_stars(tmp, RasterIO = list(nXOff = 1, nYOff = 1, nXSize = 10, nYSize = 12, nBufXSize = 2, nBufYSize = 2))) red[[1]] # cell values of subsample grid: ## Not run: plot(st, reset = FALSE, axes = TRUE, ylim = c(-.1,12.1), xlim = c(-.1,10.1), main = "nBufXSize & nBufYSize demo", text_values = TRUE) plot(st_as_sfc(red, as_points = TRUE), add = TRUE, col = 'red', pch = 16) plot(st_as_sfc(st_as_stars(st), as_points = FALSE), add = TRUE, border = 'grey') plot(st_as_sfc(red, as_points = FALSE), add = TRUE, border = 'green', lwd = 2) ## End(Not run) file.remove(tmp)
redimension array, or collapse attributes into a new dimension
## S3 method for class 'stars_proxy' st_redimension( x, new_dims = st_dimensions(x), along = list(new_dim = names(x)), ... ) st_redimension(x, new_dims, along, ...) ## S3 method for class 'stars' st_redimension( x, new_dims = st_dimensions(x), along = setNames(list(names(x)), name), ..., name = "new_dim" )
## S3 method for class 'stars_proxy' st_redimension( x, new_dims = st_dimensions(x), along = list(new_dim = names(x)), ... ) st_redimension(x, new_dims, along, ...) ## S3 method for class 'stars' st_redimension( x, new_dims = st_dimensions(x), along = setNames(list(names(x)), name), ..., name = "new_dim" )
x |
object of class |
new_dims |
target dimensions: either a ‘dimensions' object or an integer vector with the dimensions’ sizes |
along |
named list with new dimension name and values |
... |
ignored |
name |
character name of the new dimension |
st_apply apply a function to array dimensions: aggregate over space, time, or something else
## S3 method for class 'stars' st_apply( X, MARGIN, FUN, ..., CLUSTER = NULL, PROGRESS = FALSE, FUTURE = FALSE, rename = TRUE, .fname, single_arg = has_single_arg(FUN, list(...)) || can_single_arg(FUN), keep = FALSE )
## S3 method for class 'stars' st_apply( X, MARGIN, FUN, ..., CLUSTER = NULL, PROGRESS = FALSE, FUTURE = FALSE, rename = TRUE, .fname, single_arg = has_single_arg(FUN, list(...)) || can_single_arg(FUN), keep = FALSE )
X |
object of class |
MARGIN |
see apply; index number(s) or name(s) of the dimensions over which |
FUN |
see apply and see Details. |
... |
arguments passed on to |
CLUSTER |
cluster to use for parallel apply; see makeCluster |
PROGRESS |
logical; if |
FUTURE |
logical;if |
rename |
logical; if |
.fname |
function name for the new attribute name (if one or more
dimensions are reduced) or the new dimension (if a new dimension is created);
if missing, the name of |
single_arg |
logical; if |
keep |
logical; if |
FUN is a function which either operates on a single object, which will be the data of each iteration step over dimensions MARGIN, or a function that has as many arguments as there are elements in such an object. See the NDVI examples below. The second form can be VERY much faster e.g. when a trivial function is not being called for every pixel, but only once (example).
The heuristics for the default of single_arg
work often, but not always; try
setting this to the right value when st_apply
gives an error.
object of class stars
with accordingly reduced number of dimensions;
in case FUN
returns more than one value, a new dimension is created carrying
the name of the function used; see the examples. Following the logic of
apply, This new dimension is put before the
other dimensions; use aperm to rearrange this, see last example.
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) st_apply(x, 1:2, mean) # mean band value for each pixel st_apply(x, c("x", "y"), mean) # equivalent to the above st_apply(x, 3, mean) # mean of all pixels for each band ## Not run: st_apply(x, "band", mean) # equivalent to the above st_apply(x, 1:2, range) # min and max band value for each pixel fn_ndvi1 = function(x) (x[4]-x[3])/(x[4]+x[3]) # ONE argument: will be called for each pixel fn_ndvi2 = function(red,nir) (nir-red)/(nir+red) # n arguments: will be called only once ndvi1 = st_apply(x, 1:2, fn_ndvi1) # note that we can select bands 3 and 4 in the first argument: ndvi2 = st_apply(x[,,,3:4], 1:2, fn_ndvi2) all.equal(ndvi1, ndvi2) # compute the (spatial) variance of each band; https://github.com/r-spatial/stars/issues/430 st_apply(x, 3, function(x) var(as.vector(x))) # as.vector is required! # to get a progress bar also in non-interactive mode, specify: if (require(pbapply)) { # install it, if FALSE pboptions(type = "timer") } st_apply(x, 1:2, range) # dimension "range" is first; rearrange by: st_apply(x, 1:2, range) %>% aperm(c(2,3,1)) ## End(Not run)
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) st_apply(x, 1:2, mean) # mean band value for each pixel st_apply(x, c("x", "y"), mean) # equivalent to the above st_apply(x, 3, mean) # mean of all pixels for each band ## Not run: st_apply(x, "band", mean) # equivalent to the above st_apply(x, 1:2, range) # min and max band value for each pixel fn_ndvi1 = function(x) (x[4]-x[3])/(x[4]+x[3]) # ONE argument: will be called for each pixel fn_ndvi2 = function(red,nir) (nir-red)/(nir+red) # n arguments: will be called only once ndvi1 = st_apply(x, 1:2, fn_ndvi1) # note that we can select bands 3 and 4 in the first argument: ndvi2 = st_apply(x[,,,3:4], 1:2, fn_ndvi2) all.equal(ndvi1, ndvi2) # compute the (spatial) variance of each band; https://github.com/r-spatial/stars/issues/430 st_apply(x, 3, function(x) var(as.vector(x))) # as.vector is required! # to get a progress bar also in non-interactive mode, specify: if (require(pbapply)) { # install it, if FALSE pboptions(type = "timer") } st_apply(x, 1:2, range) # dimension "range" is first; rearrange by: st_apply(x, 1:2, range) %>% aperm(c(2,3,1)) ## End(Not run)
Convert stars object into an sf object
## S3 method for class 'stars' st_as_sfc(x, ..., as_points, which = seq_len(prod(dim(x)[1:2]))) ## S3 method for class 'stars' st_as_sf( x, ..., as_points = FALSE, merge = FALSE, na.rm = TRUE, use_integer = is.logical(x[[1]]) || is.integer(x[[1]]), long = FALSE, connect8 = FALSE ) ## S3 method for class 'stars_proxy' st_as_sf(x, ..., downsample = 0)
## S3 method for class 'stars' st_as_sfc(x, ..., as_points, which = seq_len(prod(dim(x)[1:2]))) ## S3 method for class 'stars' st_as_sf( x, ..., as_points = FALSE, merge = FALSE, na.rm = TRUE, use_integer = is.logical(x[[1]]) || is.integer(x[[1]]), long = FALSE, connect8 = FALSE ) ## S3 method for class 'stars_proxy' st_as_sf(x, ..., downsample = 0)
x |
object of class |
... |
ignored |
as_points |
logical; should cells be converted to points or to polygons? See details. |
which |
linear index of cells to keep (this argument is not recommended to be used) |
merge |
logical; if |
na.rm |
logical; should missing valued cells be removed, or also be converted to features? |
use_integer |
(relevant only if |
long |
logical; if |
connect8 |
logical; if |
downsample |
see st_as_stars |
If merge
is TRUE
, only the first attribute is converted into an sf
object. If na.rm
is FALSE
, areas with NA
values are also written out as polygons. Note that the resulting polygons are typically invalid, and use st_make_valid to create valid polygons out of them.
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) x = x[,1:100,1:100,6] # subset of a band with lower values in it x[[1]][x[[1]] < 30] = NA # set lower values to NA x[[1]] = x[[1]] < 100 # make the rest binary x (p = st_as_sf(x)) # removes NA areas (p = st_as_sf(x[,,,1], merge = TRUE)) # glues polygons together all(st_is_valid(p)) # not all valid, see details plot(p, axes = TRUE) (p = st_as_sf(x, na.rm = FALSE, merge = TRUE)) # includes polygons with NA values plot(p, axes = TRUE)
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) x = x[,1:100,1:100,6] # subset of a band with lower values in it x[[1]][x[[1]] < 30] = NA # set lower values to NA x[[1]] = x[[1]] < 100 # make the rest binary x (p = st_as_sf(x)) # removes NA areas (p = st_as_sf(x[,,,1], merge = TRUE)) # glues polygons together all(st_is_valid(p)) # not all valid, see details plot(p, axes = TRUE) (p = st_as_sf(x, na.rm = FALSE, merge = TRUE)) # includes polygons with NA values plot(p, axes = TRUE)
convert objects into a stars object
## S3 method for class 'cubble_df' st_as_stars(.x, ..., check_times = FALSE) ## S3 method for class 'ncdfgeom' st_as_stars(.x, ..., sf_geometry = NA) ## S3 method for class 'OpenStreetMap' st_as_stars(.x, ..., as_col = FALSE) ## S3 method for class 'stars_proxy' st_as_stars( .x, ..., downsample = 0, url = attr(.x, "url"), envir = parent.frame() ) ## S3 method for class 'data.frame' st_as_stars(.x, ..., dims = coords, xy, y_decreasing = TRUE, coords = 1:2) ## S3 method for class 'Raster' st_as_stars(.x, ..., att = 1, ignore_file = FALSE) ## S3 method for class 'SpatRaster' st_as_stars( .x, ..., ignore_file = FALSE, as_attributes = all(terra::is.factor(.x)) ) ## S3 method for class 'sf' st_as_stars(.x, ..., dims = attr(.x, "sf_column")) st_as_stars(.x, ...) ## S3 method for class 'list' st_as_stars(.x, ..., dimensions = NULL) ## Default S3 method: st_as_stars(.x = NULL, ..., raster = NULL) ## S3 method for class 'stars' st_as_stars(.x, ..., curvilinear = NULL, crs = st_crs("OGC:CRS84")) ## S3 method for class 'bbox' st_as_stars( .x, ..., nx, ny, dx = dy, dy = dx, xlim = .x[c("xmin", "xmax")], ylim = .x[c("ymin", "ymax")], values = 0, n = 64800, pretty = FALSE, inside = FALSE, nz, proxy = FALSE ) ## S3 method for class 'xts' st_as_stars(.x, ..., dimensions, name = "attr")
## S3 method for class 'cubble_df' st_as_stars(.x, ..., check_times = FALSE) ## S3 method for class 'ncdfgeom' st_as_stars(.x, ..., sf_geometry = NA) ## S3 method for class 'OpenStreetMap' st_as_stars(.x, ..., as_col = FALSE) ## S3 method for class 'stars_proxy' st_as_stars( .x, ..., downsample = 0, url = attr(.x, "url"), envir = parent.frame() ) ## S3 method for class 'data.frame' st_as_stars(.x, ..., dims = coords, xy, y_decreasing = TRUE, coords = 1:2) ## S3 method for class 'Raster' st_as_stars(.x, ..., att = 1, ignore_file = FALSE) ## S3 method for class 'SpatRaster' st_as_stars( .x, ..., ignore_file = FALSE, as_attributes = all(terra::is.factor(.x)) ) ## S3 method for class 'sf' st_as_stars(.x, ..., dims = attr(.x, "sf_column")) st_as_stars(.x, ...) ## S3 method for class 'list' st_as_stars(.x, ..., dimensions = NULL) ## Default S3 method: st_as_stars(.x = NULL, ..., raster = NULL) ## S3 method for class 'stars' st_as_stars(.x, ..., curvilinear = NULL, crs = st_crs("OGC:CRS84")) ## S3 method for class 'bbox' st_as_stars( .x, ..., nx, ny, dx = dy, dy = dx, xlim = .x[c("xmin", "xmax")], ylim = .x[c("ymin", "ymax")], values = 0, n = 64800, pretty = FALSE, inside = FALSE, nz, proxy = FALSE ) ## S3 method for class 'xts' st_as_stars(.x, ..., dimensions, name = "attr")
.x |
object to convert |
... |
in case |
check_times |
logical; should we check that the time stamps of all time series are identical? |
sf_geometry |
sf data.frame with geometry and attributes to be added to stars object. Must have same number of rows as timeseries instances. |
as_col |
logical; return rgb numbers (FALSE) or (character) color values (TRUE)? |
downsample |
integer: if larger than 0, downsample with this rate (number of pixels to skip in every row/column); if length 2, specifies downsampling rate in x and y. |
url |
character; URL of the stars endpoint where the data reside |
envir |
environment to resolve objects in |
dims |
the column names or indices that form the cube dimensions |
xy |
the x and y raster dimension names or indices; only takes effect after |
y_decreasing |
logical; if TRUE, (numeric) y values get a negative delta (decrease with increasing index) |
coords |
same as dims, for symmetry with st_as_sf |
att |
see factorValues; column in the RasterLayer's attribute table |
ignore_file |
logical; if |
as_attributes |
logical; if |
dimensions |
object of class dimensions |
raster |
character; the names of the dimensions that denote raster dimensions |
curvilinear |
only for creating curvilinear grids: named length 2 list holding longitude and latitude matrices or stars arrays, or the names of the corresponding attributes in |
crs |
object of class |
nx |
integer; number of cells in x direction; see details |
ny |
integer; number of cells in y direction; see details |
dx |
numeric or object of class units; cell size in x direction; see details |
dy |
numeric or object of class units; cell size in y direction; see details |
xlim |
length 2 numeric vector with extent (min, max) in x direction |
ylim |
length 2 numeric vector with extent (min, max) in y direction |
values |
value(s) to populate the raster values with |
n |
the (approximate) target number of grid cells |
pretty |
logical; should cell coordinates have pretty values? |
inside |
logical; should all cells entirely fall inside the bbox, potentially not covering it completely ( |
nz |
integer; number of cells in z direction; if missing no z-dimension is created. |
proxy |
logical; should a |
name |
character; attribute name for array from an |
For the ncdfgeom
method: objects are point-timeseries with optional line or polygon geometry for each timeseries specified with the sf_geometry
parameter. See ncdfgeom for more about this NetCDF-based format for geometry and timeseries.
If xy
is not specified and the first two dimensions in dims
are both numeric,
then it is set to these two dimensions.
The st_as_stars
method for sf
objects without any additional
arguments returns a one-dimensional data cube with a dimension for the simple features
geometries, and all remaining attributes as data cube attributes. When used with
further arguments, the method for data.frame
s is called.
if curvilinear
is a list with stars
objects with longitude and latitude values, its coordinate reference system is typically not that of the latitude and longitude values. If curvilinear
contains the names of two arrays in .x
, then these are removed from the returned object.
For the bbox
method: if pretty
is TRUE
, raster cells may extend the coordinate range of .x
on all sides. If in addition to nx
and ny
, dx
and dy
are also missing, these are set to a single value computed as sqrt(diff(xlim)*diff(ylim)/n)
.
If nx
and ny
are missing and values
is a matrix, the number of columns and rows of the matrix are taken.
Otherwise, if nx
and ny
are missing, they are computed as the (ceiling, floor, or rounded to integer value) of the ratio of the (x or y) range divided by (dx or dy), depending on the value of inside
. Positive dy
will be made negative. Further named arguments (...
) are passed on to pretty
. If dx
or dy
are units
objects, their value is converted to the units of st_crs(.x)
(only when sf >= 1.0-7).
for the xts
methods, if dimensions
are provided, time has to be the first dimension.
if (require(plm, quietly = TRUE)) { data(Produc, package = "plm") st_as_stars(Produc) } if (require(dplyr, quietly = TRUE)) { # https://stackoverflow.com/questions/77368957/ spatial_dim <- st_sf( ID = 1:3, geometry = list( st_polygon(list( cbind(c(0, 1, 1, 0, 0), c(0, 0, 1, 1, 0)) )), st_polygon(list( cbind(c(1, 2, 2, 1, 1), c(0, 0, 1, 1, 0)) )), st_polygon(list( cbind(c(2, 3, 3, 2, 2), c(0, 0, 1, 1, 0)) )) ) ) weekdays_dim <- data.frame(weekdays = c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")) hours_dim <- data.frame(hours = c("8am", "11am", "4pm", "11pm")) sf_dta <- spatial_dim |> cross_join(weekdays_dim)|> cross_join(hours_dim) |> mutate(population = rnorm(n(), mean = 1000, sd = 200)) |> select(everything(), geometry) st_as_stars(sf_dta, dims = c("weekdays", "hours", "geometry")) } demo(nc, echo=FALSE,ask=FALSE) st_as_stars(nc) st_as_stars(st_drop_geometry(nc), dims = "NAME") data.frame(expand.grid(x=1:5, y = 1:5), z = rnorm(25)) |> st_as_stars() nc = st_read(system.file("gpkg/nc.gpkg", package="sf")) st_as_stars(nc)
if (require(plm, quietly = TRUE)) { data(Produc, package = "plm") st_as_stars(Produc) } if (require(dplyr, quietly = TRUE)) { # https://stackoverflow.com/questions/77368957/ spatial_dim <- st_sf( ID = 1:3, geometry = list( st_polygon(list( cbind(c(0, 1, 1, 0, 0), c(0, 0, 1, 1, 0)) )), st_polygon(list( cbind(c(1, 2, 2, 1, 1), c(0, 0, 1, 1, 0)) )), st_polygon(list( cbind(c(2, 3, 3, 2, 2), c(0, 0, 1, 1, 0)) )) ) ) weekdays_dim <- data.frame(weekdays = c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday")) hours_dim <- data.frame(hours = c("8am", "11am", "4pm", "11pm")) sf_dta <- spatial_dim |> cross_join(weekdays_dim)|> cross_join(hours_dim) |> mutate(population = rnorm(n(), mean = 1000, sd = 200)) |> select(everything(), geometry) st_as_stars(sf_dta, dims = c("weekdays", "hours", "geometry")) } demo(nc, echo=FALSE,ask=FALSE) st_as_stars(nc) st_as_stars(st_drop_geometry(nc), dims = "NAME") data.frame(expand.grid(x=1:5, y = 1:5), z = rnorm(25)) |> st_as_stars() nc = st_read(system.file("gpkg/nc.gpkg", package="sf")) st_as_stars(nc)
If the object has been cropped without normalization, then the indices return
are relative to the original uncropped extent. See st_crop
st_cells(x, sf)
st_cells(x, sf)
x |
object of class |
sf |
object of class |
set.seed(1345) st_bbox(L7_ETMs) |> st_as_sfc() |> st_sample(10) -> pts (x <- st_cells(L7_ETMs, pts)) # get the pixel values (first band only): st_as_stars(L7_ETMs)[[1]][x] # get pixel values for all bands: st_as_stars(L7_ETMs) |> split() |> sapply(`[`, x) # compare with st_extract(): st_as_stars(L7_ETMs) |> split() |> st_extract(pts)
set.seed(1345) st_bbox(L7_ETMs) |> st_as_sfc() |> st_sample(10) -> pts (x <- st_cells(L7_ETMs, pts)) # get the pixel values (first band only): st_as_stars(L7_ETMs)[[1]][x] # get pixel values for all bands: st_as_stars(L7_ETMs) |> split() |> sapply(`[`, x) # compare with st_extract(): st_as_stars(L7_ETMs) |> split() |> st_extract(pts)
Compute contour lines or sets
st_contour( x, na.rm = TRUE, contour_lines = FALSE, breaks = classInt::classIntervals(na.omit(as.vector(x[[1]])))$brks )
st_contour( x, na.rm = TRUE, contour_lines = FALSE, breaks = classInt::classIntervals(na.omit(as.vector(x[[1]])))$brks )
x |
object of class |
na.rm |
logical; should missing valued cells be removed, or also be converted to features? |
contour_lines |
logical; if |
breaks |
numerical; values at which to "draw" contour levels |
this function requires GDAL >= 2.4.0
for polygonizing rasters following grid boundaries, see st_as_sf with arguments as_points=FALSE
and merge=TRUE
; contour plots contour lines using R's native algorithm (which also plots contour levels)
retrieve coordinates for raster or vector cube cells
## S3 method for class 'stars' st_coordinates(x, ..., add_max = FALSE, center = TRUE) ## S3 method for class 'stars' as.data.frame(x, ..., add_max = FALSE, center = NA, add_coordinates = TRUE) as_tibble.stars(.x, ..., add_max = FALSE, center = NA)
## S3 method for class 'stars' st_coordinates(x, ..., add_max = FALSE, center = TRUE) ## S3 method for class 'stars' as.data.frame(x, ..., add_max = FALSE, center = NA, add_coordinates = TRUE) as_tibble.stars(.x, ..., add_max = FALSE, center = NA)
x |
object of class |
... |
ignored |
add_max |
logical; if |
center |
logical; (only if |
add_coordinates |
logical; if 'TRUE', columns with dimension values preceed the array values, otherwise they are omitted |
.x |
object to be converted to a tibble |
crop a stars object
## S3 method for class 'mdim' st_crop(x, y, ...) ## S3 method for class 'stars_proxy' st_crop( x, y, ..., crop = TRUE, epsilon = sqrt(.Machine$double.eps), collect = TRUE ) ## S3 method for class 'stars' st_crop( x, y, ..., crop = TRUE, epsilon = sqrt(.Machine$double.eps), as_points = all(st_dimension(y) == 2, na.rm = TRUE), normalize = FALSE )
## S3 method for class 'mdim' st_crop(x, y, ...) ## S3 method for class 'stars_proxy' st_crop( x, y, ..., crop = TRUE, epsilon = sqrt(.Machine$double.eps), collect = TRUE ) ## S3 method for class 'stars' st_crop( x, y, ..., crop = TRUE, epsilon = sqrt(.Machine$double.eps), as_points = all(st_dimension(y) == 2, na.rm = TRUE), normalize = FALSE )
x |
object of class |
y |
object of class |
... |
ignored |
crop |
logical; if |
epsilon |
numeric; factor to shrink the bounding box of |
collect |
logical; if |
as_points |
logical; only relevant if |
normalize |
logical; if |
for raster x
, st_crop
selects cells that intersect with y
.
For intersection, are raster cells interpreted as points or as small polygons?
If y
is of class stars
, x
raster cells are interpreted as points; if y
is of class bbox
, x
cells are interpreted as cells (small polygons). Otherwise, if as_points
is not given, cells are interpreted as points if y
has a two-dimensional geometry.
l7 = read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) d = st_dimensions(l7) # area around cells 3:10 (x) and 4:11 (y): offset = c(d[["x"]]$offset, d[["y"]]$offset) res = c(d[["x"]]$delta, d[["y"]]$delta) bb = st_bbox(c(xmin = offset[1] + 2 * res[1], ymin = offset[2] + 11 * res[2], xmax = offset[1] + 10 * res[1], ymax = offset[2] + 3 * res[2]), crs = st_crs(l7)) l7[bb] # equivalent: st_crop(l7, bb) plot(l7[,1:13,1:13,1], reset = FALSE) image(l7[bb,,,1], add = TRUE, col = sf.colors()) plot(st_as_sfc(bb), add = TRUE, border = 'green', lwd = 2) # slightly smaller bbox: bb = st_bbox(c(xmin = offset[1] + 2.1 * res[1], ymin = offset[2] + 10.9 * res[2], xmax = offset[1] + 9.9 * res[1], ymax = offset[2] + 3.1 * res[2]), crs = st_crs(l7)) l7[bb] plot(l7[,1:13,1:13,1], reset = FALSE) image(l7[bb,,,1], add = TRUE, col = sf.colors()) plot(st_as_sfc(bb), add = TRUE, border = 'green', lwd = 2) # slightly larger bbox: bb = st_bbox(c(xmin = offset[1] + 1.9 * res[1], ymin = offset[2] + 11.1 * res[2], xmax = offset[1] + 10.1 * res[1], ymax = offset[2] + 2.9 * res[2]), crs = st_crs(l7)) l7[bb] plot(l7[,1:13,1:13,1], reset = FALSE) image(l7[bb,,,1], add = TRUE, col = sf.colors()) plot(st_as_sfc(bb), add = TRUE, border = 'green', lwd = 2) # half a cell size larger bbox: bb = st_bbox(c(xmin = offset[1] + 1.49 * res[1], ymin = offset[2] + 11.51 * res[2], xmax = offset[1] + 10.51 * res[1], ymax = offset[2] + 2.49 * res[2]), crs = st_crs(l7)) l7[bb] plot(l7[,1:13,1:13,1], reset = FALSE) image(l7[bb,,,1], add = TRUE, col = sf.colors()) plot(st_as_sfc(bb), add = TRUE, border = 'green', lwd = 2)
l7 = read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) d = st_dimensions(l7) # area around cells 3:10 (x) and 4:11 (y): offset = c(d[["x"]]$offset, d[["y"]]$offset) res = c(d[["x"]]$delta, d[["y"]]$delta) bb = st_bbox(c(xmin = offset[1] + 2 * res[1], ymin = offset[2] + 11 * res[2], xmax = offset[1] + 10 * res[1], ymax = offset[2] + 3 * res[2]), crs = st_crs(l7)) l7[bb] # equivalent: st_crop(l7, bb) plot(l7[,1:13,1:13,1], reset = FALSE) image(l7[bb,,,1], add = TRUE, col = sf.colors()) plot(st_as_sfc(bb), add = TRUE, border = 'green', lwd = 2) # slightly smaller bbox: bb = st_bbox(c(xmin = offset[1] + 2.1 * res[1], ymin = offset[2] + 10.9 * res[2], xmax = offset[1] + 9.9 * res[1], ymax = offset[2] + 3.1 * res[2]), crs = st_crs(l7)) l7[bb] plot(l7[,1:13,1:13,1], reset = FALSE) image(l7[bb,,,1], add = TRUE, col = sf.colors()) plot(st_as_sfc(bb), add = TRUE, border = 'green', lwd = 2) # slightly larger bbox: bb = st_bbox(c(xmin = offset[1] + 1.9 * res[1], ymin = offset[2] + 11.1 * res[2], xmax = offset[1] + 10.1 * res[1], ymax = offset[2] + 2.9 * res[2]), crs = st_crs(l7)) l7[bb] plot(l7[,1:13,1:13,1], reset = FALSE) image(l7[bb,,,1], add = TRUE, col = sf.colors()) plot(st_as_sfc(bb), add = TRUE, border = 'green', lwd = 2) # half a cell size larger bbox: bb = st_bbox(c(xmin = offset[1] + 1.49 * res[1], ymin = offset[2] + 11.51 * res[2], xmax = offset[1] + 10.51 * res[1], ymax = offset[2] + 2.49 * res[2]), crs = st_crs(l7)) l7[bb] plot(l7[,1:13,1:13,1], reset = FALSE) image(l7[bb,,,1], add = TRUE, col = sf.colors()) plot(st_as_sfc(bb), add = TRUE, border = 'green', lwd = 2)
create an array with dimension values
st_dim_to_attr(x, which = seq_along(dim(x)))
st_dim_to_attr(x, which = seq_along(dim(x)))
x |
object of class |
which |
integer; indices of the dimensions to address (default: all) |
stars
object with dimension values as attributes
tif = system.file("tif/L7_ETMs.tif", package = "stars") x1 = read_stars(tif) (x = st_dim_to_attr(x1)) plot(x) (x = st_dim_to_attr(x1, 2:3)) plot(x) (x= st_dim_to_attr(x1, 3)) plot(x)
tif = system.file("tif/L7_ETMs.tif", package = "stars") x1 = read_stars(tif) (x = st_dim_to_attr(x1)) plot(x) (x = st_dim_to_attr(x1, 2:3)) plot(x) (x= st_dim_to_attr(x1, 3)) plot(x)
get dimensions from stars object
st_dimensions(.x, ...) ## S3 method for class 'stars' st_dimensions(.x, ...) st_dimensions(x) <- value ## S3 replacement method for class 'stars' st_dimensions(x) <- value ## S3 replacement method for class 'stars_proxy' st_dimensions(x) <- value ## S3 replacement method for class 'list' st_dimensions(x) <- value ## S3 method for class 'array' st_dimensions(.x, ...) ## Default S3 method: st_dimensions( .x, ..., .raster, affine = c(0, 0), cell_midpoints = FALSE, point = FALSE ) st_set_dimensions( .x, which, values = NULL, point = NULL, names = NULL, xy, ... ) st_get_dimension_values(.x, which, ..., where = NA, max = FALSE, center = NA)
st_dimensions(.x, ...) ## S3 method for class 'stars' st_dimensions(.x, ...) st_dimensions(x) <- value ## S3 replacement method for class 'stars' st_dimensions(x) <- value ## S3 replacement method for class 'stars_proxy' st_dimensions(x) <- value ## S3 replacement method for class 'list' st_dimensions(x) <- value ## S3 method for class 'array' st_dimensions(.x, ...) ## Default S3 method: st_dimensions( .x, ..., .raster, affine = c(0, 0), cell_midpoints = FALSE, point = FALSE ) st_set_dimensions( .x, which, values = NULL, point = NULL, names = NULL, xy, ... ) st_get_dimension_values(.x, which, ..., where = NA, max = FALSE, center = NA)
.x |
object to retrieve dimensions information from |
... |
further arguments |
x |
object of class |
value |
new object of class |
.raster |
length 2 character array with names (if any) of the raster dimensions |
affine |
numeric; specify parameters of the affine transformation |
cell_midpoints |
logical; if |
point |
logical; does the pixel value (measure) refer to a point (location) value or to an pixel (area) summary value? |
which |
integer or character; index or name of the dimension to be changed |
values |
values for this dimension (e.g. |
names |
character; vector with new names for all dimensions, or with the single new name for the dimension indicated by |
xy |
length-2 character vector; (new) names for the |
where |
character, one of 'start', 'center' or 'end'. Set to NA (default) to ignore and use |
max |
logical; if |
center |
logical; if |
dimensions can be specified in two ways. The simplest is to pass a vector with numeric values for a numeric dimension, or character values for a categorical dimension. Parameter cell_midpoints
is used to specify whether numeric values refer to the offset (start) of a dimension interval (default), or to the center; the center case is only available for regular dimensions. For rectilinear numeric dimensions, one can specify either a vector with cell borders (start values), or a data.frame with two columns named "start" and "end", with the respective interval start and end values. In the first case, the end values are computed from the start values by assuming the last two intervals have equal width.
the dimensions
attribute of x
, of class dimensions
x = read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) # Landsat 7 ETM+ band semantics: https://landsat.gsfc.nasa.gov/the-enhanced-thematic-mapper-plus/ # set bands to values 1,2,3,4,5,7: (x1 = st_set_dimensions(x, "band", values = c(1,2,3,4,5,7), names = "band_number", point = TRUE)) # set band values as bandwidth rbind(c(0.45,0.515), c(0.525,0.605), c(0.63,0.69), c(0.775,0.90), c(1.55,1.75), c(2.08,2.35)) %>% units::set_units("um") -> bw # or: units::set_units(µm) -> bw # set bandwidth midpoint: (x2 = st_set_dimensions(x, "band", values = 0.5 * (bw[,1]+bw[,2]), names = "bandwidth_midpoint", point = TRUE)) # set bandwidth intervals: (x3 = st_set_dimensions(x, "band", values = make_intervals(bw), names = "bandwidth")) m = matrix(1:20, nrow = 5, ncol = 4) dim(m) = c(x = 5, y = 4) # named dim (s = st_as_stars(m)) st_get_dimension_values(s, 'x', where = "start") st_get_dimension_values(s, 'x', center = FALSE) st_get_dimension_values(s, 'x', where = "center") st_get_dimension_values(s, 'x', center = TRUE) st_get_dimension_values(s, 'x', where = "end") st_get_dimension_values(s, 'x', max = TRUE)
x = read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) # Landsat 7 ETM+ band semantics: https://landsat.gsfc.nasa.gov/the-enhanced-thematic-mapper-plus/ # set bands to values 1,2,3,4,5,7: (x1 = st_set_dimensions(x, "band", values = c(1,2,3,4,5,7), names = "band_number", point = TRUE)) # set band values as bandwidth rbind(c(0.45,0.515), c(0.525,0.605), c(0.63,0.69), c(0.775,0.90), c(1.55,1.75), c(2.08,2.35)) %>% units::set_units("um") -> bw # or: units::set_units(µm) -> bw # set bandwidth midpoint: (x2 = st_set_dimensions(x, "band", values = 0.5 * (bw[,1]+bw[,2]), names = "bandwidth_midpoint", point = TRUE)) # set bandwidth intervals: (x3 = st_set_dimensions(x, "band", values = make_intervals(bw), names = "bandwidth")) m = matrix(1:20, nrow = 5, ncol = 4) dim(m) = c(x = 5, y = 4) # named dim (s = st_as_stars(m)) st_get_dimension_values(s, 'x', where = "start") st_get_dimension_values(s, 'x', center = FALSE) st_get_dimension_values(s, 'x', where = "center") st_get_dimension_values(s, 'x', center = TRUE) st_get_dimension_values(s, 'x', where = "end") st_get_dimension_values(s, 'x', max = TRUE)
downsample a stars or stars_proxy object either by skipping rows, columns and bands, or by computing a single value (e.g. the mean) from the sub-tiles involved
st_downsample(x, n, ...) ## S3 method for class 'stars' st_downsample(x, n, ..., offset = 0, FUN) ## S3 method for class 'stars_proxy' st_downsample(x, n, ...)
st_downsample(x, n, ...) ## S3 method for class 'stars' st_downsample(x, n, ..., offset = 0, FUN) ## S3 method for class 'stars_proxy' st_downsample(x, n, ...)
x |
object of class stars or stars_proxy |
n |
integer; for each dimension the number of pixels/lines/bands etc that will be skipped; see Details. |
... |
arguments passed on to |
offset |
integer; offset(s) for downsampling, in pixels, starting at the offset of
each dimension; should be smaller or equal to |
FUN |
function; if given, downsampling will apply FUN to each of the the subtiles |
If all n == 0, no downsampling takes place; if it is 1, every second row/column/band is skipped, if it is 2, every second+third row/column/band are skipped, etc.
Downsampling a stars_proxy
object returns a stars
object, is
equivalent to calling st_as_stars(x, downsample = 2)
, and only downsamples
the first two (x and y) dimensions.
Downsampled regular rasters keep their dimension offsets, have a cell size (delta) that is n[i]+1 times larger, and may result in a (slightly) different extent.
Note that terra's aggregate with fact=2
corresponds to
st_downsample(x, n = 1, FUN = mean)
: fact
is one larger than n
.
(m = matrix(1:121, 11, 11)) (s = st_as_stars(m)) st_downsample(s, 1) st_downsample(s, 1)[[1]] st_downsample(s, 1, offset = 1) st_downsample(s, 1, offset = 1)[[1]] st_downsample(s, 1, offset = c(0,1)) st_downsample(s, 1, offset = c(0,1))[[1]] st_downsample(s, 1, FUN = mean) st_downsample(s, 1, FUN = mean)[[1]] st_downsample(s, 1, offset = 1, FUN = mean) st_downsample(s, 1, offset = c(0,1), FUN = mean)[[1]]
(m = matrix(1:121, 11, 11)) (s = st_as_stars(m)) st_downsample(s, 1) st_downsample(s, 1)[[1]] st_downsample(s, 1, offset = 1) st_downsample(s, 1, offset = 1)[[1]] st_downsample(s, 1, offset = c(0,1)) st_downsample(s, 1, offset = c(0,1))[[1]] st_downsample(s, 1, FUN = mean) st_downsample(s, 1, FUN = mean)[[1]] st_downsample(s, 1, offset = 1, FUN = mean) st_downsample(s, 1, offset = c(0,1), FUN = mean)[[1]]
Extract cell values at point locations
st_extract(x, ...) ## S3 method for class 'stars' st_extract( x, at, ..., bilinear = FALSE, time_column = attr(at, "time_column") %||% attr(at, "time_col"), interpolate_time = bilinear, FUN = mean )
st_extract(x, ...) ## S3 method for class 'stars' st_extract( x, at, ..., bilinear = FALSE, time_column = attr(at, "time_column") %||% attr(at, "time_col"), interpolate_time = bilinear, FUN = mean )
x |
object of class |
... |
passed on to aggregate.stars when geometries are not exclusively POINT geometries |
at |
object of class |
bilinear |
logical; use bilinear interpolation rather than nearest neighbour? |
time_column |
character or integer; name or index of a column with time or date values that will be matched to values of the first temporal dimension (matching classes |
interpolate_time |
logical; should time be interpolated? if FALSE, time instances are matched using the coinciding or the last preceding time in the data cube. |
FUN |
function used to aggregate pixel values when geometries of |
points outside the raster are returned as NA
values. For
large sets of points for which extraction is needed, passing a matrix as
to at
may be much faster than passing an sf
or sfc
object.
if at
is of class matrix
, a matrix with extracted values is returned;
otherwise: if x
has more dimensions than only x and y (raster), an
object of class stars
with POINT geometries replacing x and y raster
dimensions, if this is not the case, an object of sf
with extracted values.
tif = system.file("tif/L7_ETMs.tif", package = "stars") r = read_stars(tif) pnt = st_sample(st_as_sfc(st_bbox(r)), 10) st_extract(r, pnt) st_extract(r, pnt) %>% st_as_sf() st_extract(r[,,,1], pnt) st_extract(r, st_coordinates(pnt)) # "at" is a matrix: return a matrix
tif = system.file("tif/L7_ETMs.tif", package = "stars") r = read_stars(tif) pnt = st_sample(st_as_sfc(st_bbox(r)), 10) st_extract(r, pnt) st_extract(r, pnt) %>% st_as_sf() st_extract(r[,,,1], pnt) st_extract(r, st_coordinates(pnt)) # "at" is a matrix: return a matrix
get or set the geotransform, or rotation matrix
st_geotransform(x, ...) st_geotransform(x) <- value ## S3 replacement method for class 'stars' st_geotransform(x) <- value
st_geotransform(x, ...) st_geotransform(x) <- value ## S3 replacement method for class 'stars' st_geotransform(x) <- value
x |
object of class stars or dimensions |
... |
ignored |
value |
length 6 numeric vector, or 2 x 2 (scaled) rotation matrix |
# using the "classical" rotation matrix, see https://en.wikipedia.org/wiki/Rotation_matrix : rot = function(theta, dxdy = c(1., -1.)) { th = theta / 180 * pi matrix(c(cos(th), sin(th), -sin(th), cos(th)), 2, 2) %*% matrix(c(dxdy[2], 0, 0, dxdy[1]), 2, 2) } l = st_downsample(st_as_stars(L7_ETMs), 9) # save time in plotting st_geotransform(l) = rot(20, c(28.5, 28.5)) # clockwise, 20 degrees, scale by cell size plot(l[,,,1]) m = rot(20, c(1, 2)) g = expand.grid(x = 0:4, y = 0:4) plot(g[1:2], asp = 1) text(g[,1], g[,2], labels = seq_along(g[,1]), pos = 4) g = t(m %*% t(as.matrix(g))) points(g, col = 'red') text(g[,1], g[,2], labels = seq_along(g[,1]), pos = 4, col = 'red') m = matrix(1:20, 4) s0 = st_as_stars(m) s = s0 # dy > 0, clockwise rotation: st_geotransform(s) = rot(10, c(1,1)) plot(s0, reset = FALSE) plot(s, add = TRUE) # dy < 0, counter clockwise rotation, + expansion in x-direction: layout(1) s0 = st_as_stars(st_bbox(s0), dx = 1) s0$values = 1:20 s0 plot(s0, reset = FALSE) s = s0 st_geotransform(s) = rot(10, c(2,1)) plot(s, add = TRUE)
# using the "classical" rotation matrix, see https://en.wikipedia.org/wiki/Rotation_matrix : rot = function(theta, dxdy = c(1., -1.)) { th = theta / 180 * pi matrix(c(cos(th), sin(th), -sin(th), cos(th)), 2, 2) %*% matrix(c(dxdy[2], 0, 0, dxdy[1]), 2, 2) } l = st_downsample(st_as_stars(L7_ETMs), 9) # save time in plotting st_geotransform(l) = rot(20, c(28.5, 28.5)) # clockwise, 20 degrees, scale by cell size plot(l[,,,1]) m = rot(20, c(1, 2)) g = expand.grid(x = 0:4, y = 0:4) plot(g[1:2], asp = 1) text(g[,1], g[,2], labels = seq_along(g[,1]), pos = 4) g = t(m %*% t(as.matrix(g))) points(g, col = 'red') text(g[,1], g[,2], labels = seq_along(g[,1]), pos = 4, col = 'red') m = matrix(1:20, 4) s0 = st_as_stars(m) s = s0 # dy > 0, clockwise rotation: st_geotransform(s) = rot(10, c(1,1)) plot(s0, reset = FALSE) plot(s, add = TRUE) # dy < 0, counter clockwise rotation, + expansion in x-direction: layout(1) s0 = st_as_stars(st_bbox(s0), dx = 1) s0$values = 1:20 s0 plot(s0, reset = FALSE) s = s0 st_geotransform(s) = rot(10, c(2,1)) plot(s, add = TRUE)
spatial intersect predicate for stars and sfc object
## S3 method for class 'stars' st_intersects(x, y, sparse = TRUE, ..., as_points = NA, transpose = FALSE)
## S3 method for class 'stars' st_intersects(x, y, sparse = TRUE, ..., as_points = NA, transpose = FALSE)
x |
object of class stars |
y |
object that has an 'st_geometry' method: of class 'sf' or 'sfc', or 'stars' object with an 'sfc' dimension |
sparse |
logical; if TRUE, return the a sparse logical matrix (object of class 'sgbp'), if FALSE, return a logical matrix |
... |
ignored, or passed on to 'st_intersects.sf' for curvilinear grids |
as_points |
logical, should grid cells be considered as points (TRUE) or polygons (FALSE)? Default: FALSE and warning emitted |
transpose |
logical; should the transpose of the 'sgbp' object be returned? |
curvilinear grids are always converted to polygons, so points on grid boundaries may intersect with two cells touched; for other grids each cell boundary or corner belongs only to one cell.
'sgbp' object if sparse = TRUE, logical matrix otherwise
Spatially join a stars and an 'sf' object
## S3 method for class 'stars' st_join( x, y, join = st_intersects, ..., what = "left1", as_points = NA, warn = TRUE )
## S3 method for class 'stars' st_join( x, y, join = st_intersects, ..., what = "left1", as_points = NA, warn = TRUE )
x |
object of class stars |
y |
object of class sf, or one that can be coerced into that by st_as_sf |
join |
the join function, which should return an sgbp object; see details |
... |
arguments that will be passed on to the join function |
what |
"left1", "right" or "inner"; see details |
as_points |
logical; controls whether grid cells in |
warn |
logical; if TRUE, warn on 1-to-many matches when |
When there is more than one match to a single x value, the first matching record from y is taken (and if warn
is TRUE a warning is raised). If what is "inner", an object of class sf
with all matching records of x and y.
If what is "left1", an object of class stars with the (first) value of y at spatial instances of x
build mosaic (composite) of several spatially disjoint stars objects
st_mosaic(.x, ...) ## S3 method for class 'stars' st_mosaic( .x, ..., dst = tempfile(fileext = file_ext), options = c("-vrtnodata", "-9999", "-srcnodata", "nan"), file_ext = ".tif" ) ## S3 method for class 'character' st_mosaic( .x, ..., dst = tempfile(fileext = file_ext), options = c("-vrtnodata", "-9999"), file_ext = ".tif" ) ## S3 method for class 'stars_proxy' st_mosaic( .x, ..., dst = tempfile(fileext = file_ext), options = c("-vrtnodata", "-9999"), file_ext = ".tif" )
st_mosaic(.x, ...) ## S3 method for class 'stars' st_mosaic( .x, ..., dst = tempfile(fileext = file_ext), options = c("-vrtnodata", "-9999", "-srcnodata", "nan"), file_ext = ".tif" ) ## S3 method for class 'character' st_mosaic( .x, ..., dst = tempfile(fileext = file_ext), options = c("-vrtnodata", "-9999"), file_ext = ".tif" ) ## S3 method for class 'stars_proxy' st_mosaic( .x, ..., dst = tempfile(fileext = file_ext), options = c("-vrtnodata", "-9999"), file_ext = ".tif" )
.x |
object of class stars, or character vector with input dataset names |
... |
further input stars objects |
dst |
character; destination file name; this will be a VRT file with references to the source file(s), see details |
options |
character; options to the gdalbuildvrt command |
file_ext |
character; file extension, determining the format used to write to (".tif" implies GeoTIFF) |
the gdal function buildvrt builds a mosaic of input images; these input images can be multi-band, but not higher-dimensional data cubes or stars objects with multiple attributes; note that for the 'stars' method, the 'dst' file may contain references to temporary files that are going to be removed at termination of the R session.
uses gdal_utils to internally call buildvrt
; no executables external to R are called.
the stars method returns a stars object with the composite of the input; the character
method returns the file name of the file with the mosaic; see also the GDAL documentation of gdalbuildvrt
x = read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) x1 = x[,100:200,100:200,] x2 = x[,150:300,150:300,] plot(st_mosaic(x1, x2))
x = read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) x1 = x[,100:200,100:200,] x2 = x[,150:300,150:300,] plot(st_mosaic(x1, x2))
get the raster type (if any) of a stars object
st_raster_type(x, dimension = character(0))
st_raster_type(x, dimension = character(0))
x |
object of class |
dimension |
optional: numbers or names of dimension(s) to get per-dimension type |
categories "curvilinear"
and "affine"
only refer to
the relationship between a pair of spatial (raster) dimensions.
if dimension
is not specified, return the spatial raster type:
one of NA
(if the object does not have raster dimensions),
"curvilinear"
, "rectilinear"
, "affine"
, or "regular"
.
In case dimension(s) are specified, return one of "regular"
, "rectilinear"
(irregular but numeric), or "discrete"
(anything else).
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) st_raster_type(x) st_raster_type(x, 1:3)
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) st_raster_type(x) st_raster_type(x, 1:3)
rasterize simple feature geometries
st_rasterize( sf, template = guess_raster(sf, ...) %||% st_as_stars(st_bbox(sf), values = NA_real_, ...), file = tempfile(), driver = "GTiff", options = character(0), align = FALSE, proxy = FALSE, ... )
st_rasterize( sf, template = guess_raster(sf, ...) %||% st_as_stars(st_bbox(sf), values = NA_real_, ...), file = tempfile(), driver = "GTiff", options = character(0), align = FALSE, proxy = FALSE, ... )
sf |
object of class |
template |
optional; stars object with desired target geometry, or target geometry alignment if |
file |
temporary file name |
driver |
driver for temporary file |
options |
character; options vector for |
align |
logical; if |
proxy |
logical; should a proxy object be returned? |
... |
arguments passed on to st_as_stars |
if 'template' is a 'stars' object, non-NA cells that are not covered by 'sf' receive the value in 'template'; see also argument 'align'.
demo(nc, echo = FALSE, ask = FALSE) (x = st_rasterize(nc)) # default grid: plot(x, axes = TRUE) # a bit more customized grid: (x = st_rasterize(nc, st_as_stars(st_bbox(nc), nx = 100, ny = 50, values = NA_real_))) plot(x, axes = TRUE) (ls = st_sf(a = 1:2, st_sfc(st_linestring(rbind(c(0.1, 0), c(1.1, 1))), st_linestring(rbind(c(0, 0.05), c(1, 0.05)))))) (grd = st_as_stars(st_bbox(ls), nx = 10, ny = 10, xlim = c(0, 1.0), ylim = c(0, 1), values = NA_real_)) # Only the left-top corner is part of the grid cell: sf_extSoftVersion()["GDAL"] plot(st_rasterize(ls, grd), axes = TRUE, reset = FALSE) # ALL_TOUCHED=FALSE; plot(ls, add = TRUE, col = "red") plot(st_rasterize(ls, grd, options = "ALL_TOUCHED=TRUE"), axes = TRUE, reset = FALSE) plot(ls, add = TRUE, col = "red") # add lines to existing 0 values, summing values in case of multiple lines: (grd = st_as_stars(st_bbox(ls), nx = 10, ny = 10, xlim = c(0, 1.0), ylim = c(0, 1), values = 0)) r = st_rasterize(ls, grd, options = c("MERGE_ALG=ADD", "ALL_TOUCHED=TRUE")) plot(r, axes = TRUE, reset = FALSE) plot(ls, add = TRUE, col = "red")
demo(nc, echo = FALSE, ask = FALSE) (x = st_rasterize(nc)) # default grid: plot(x, axes = TRUE) # a bit more customized grid: (x = st_rasterize(nc, st_as_stars(st_bbox(nc), nx = 100, ny = 50, values = NA_real_))) plot(x, axes = TRUE) (ls = st_sf(a = 1:2, st_sfc(st_linestring(rbind(c(0.1, 0), c(1.1, 1))), st_linestring(rbind(c(0, 0.05), c(1, 0.05)))))) (grd = st_as_stars(st_bbox(ls), nx = 10, ny = 10, xlim = c(0, 1.0), ylim = c(0, 1), values = NA_real_)) # Only the left-top corner is part of the grid cell: sf_extSoftVersion()["GDAL"] plot(st_rasterize(ls, grd), axes = TRUE, reset = FALSE) # ALL_TOUCHED=FALSE; plot(ls, add = TRUE, col = "red") plot(st_rasterize(ls, grd, options = "ALL_TOUCHED=TRUE"), axes = TRUE, reset = FALSE) plot(ls, add = TRUE, col = "red") # add lines to existing 0 values, summing values in case of multiple lines: (grd = st_as_stars(st_bbox(ls), nx = 10, ny = 10, xlim = c(0, 1.0), ylim = c(0, 1), values = 0)) r = st_rasterize(ls, grd, options = c("MERGE_ALG=ADD", "ALL_TOUCHED=TRUE")) plot(r, axes = TRUE, reset = FALSE) plot(ls, add = TRUE, col = "red")
obtain resolution(s) of a stars object: by default only the (absolute) x/y raster dimensions, optionally all delta
dimension parameters
st_res(x, all = FALSE, absolute = !all)
st_res(x, all = FALSE, absolute = !all)
x |
an object of class |
all |
logical; if FALSE return a vector with the x/y raster resolution |
absolute |
logical; only works when |
if all = FALSE
a vector with x/y raster resolutions, otherwise a list with delta values
st_res(L7_ETMs) st_res(L7_ETMs, absolute = FALSE) st_res(L7_ETMs, all = TRUE) if (require(starsdata)) { paste0("netcdf/", c("avhrr-only-v2.19810901.nc", "avhrr-only-v2.19810902.nc", "avhrr-only-v2.19810903.nc", "avhrr-only-v2.19810904.nc")) |> system.file(package = "starsdata") |> read_stars(quiet = TRUE) -> x st_res(x) |> print() st_res(x, all = TRUE) |> print() }
st_res(L7_ETMs) st_res(L7_ETMs, absolute = FALSE) st_res(L7_ETMs, all = TRUE) if (require(starsdata)) { paste0("netcdf/", c("avhrr-only-v2.19810901.nc", "avhrr-only-v2.19810902.nc", "avhrr-only-v2.19810903.nc", "avhrr-only-v2.19810904.nc")) |> system.file(package = "starsdata") |> read_stars(quiet = TRUE) -> x st_res(x) |> print() st_res(x, all = TRUE) |> print() }
reduce dimension to rgb (alpha) hex values
st_rgb( x, dimension = 3, use_alpha = dim(x)[dimension] == 4, maxColorValue = 255L, probs = c(0, 1), stretch = NULL )
st_rgb( x, dimension = 3, use_alpha = dim(x)[dimension] == 4, maxColorValue = 255L, probs = c(0, 1), stretch = NULL )
x |
object of class |
dimension |
dimension name or number to reduce |
use_alpha |
logical; if TRUE, the fourth band will be used as alpha values |
maxColorValue |
integer; maximum value for colors |
probs |
probability values for quantiles used for stretching by "percent". |
stretch |
logical or character; if |
the dimension's bands are mapped to red, green, blue, alpha; if a different
ordering is wanted, use [.stars to reorder a dimension, see examples.
Alternatively, you can use plot.stars with the rgb
argument to create a three-band composition.
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) st_rgb(x[,,,3:1]) r = st_rgb(x[,,,c(6,5,4,3)], 3, use_alpha=TRUE) # now R=6,G=5,B=4,alpha=3 if (require(ggplot2)) { ggplot() + geom_stars(data = r) + scale_fill_identity() } r = st_rgb(x[,,,3:1], probs = c(0.01, 0.99), stretch = "percent") plot(r) r = st_rgb(x[,,,3:1], probs = c(0.01, 0.99), stretch = "histogram") plot(r)
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) st_rgb(x[,,,3:1]) r = st_rgb(x[,,,c(6,5,4,3)], 3, use_alpha=TRUE) # now R=6,G=5,B=4,alpha=3 if (require(ggplot2)) { ggplot() + geom_stars(data = r) + scale_fill_identity() } r = st_rgb(x[,,,3:1], probs = c(0.01, 0.99), stretch = "percent") plot(r) r = st_rgb(x[,,,3:1], probs = c(0.01, 0.99), stretch = "histogram") plot(r)
Transform rotated long/lat regular grid to unrotated curvilinear grid
## S3 method for class 'stars' st_rotate(.x, lon0, lat0, north = TRUE, ...) ## S3 method for class 'sfc' st_rotate(.x, lon0, lat0, north = TRUE, ...) ## S3 method for class 'sf' st_rotate(.x, lon0, lat0, north = TRUE, ...)
## S3 method for class 'stars' st_rotate(.x, lon0, lat0, north = TRUE, ...) ## S3 method for class 'sfc' st_rotate(.x, lon0, lat0, north = TRUE, ...) ## S3 method for class 'sf' st_rotate(.x, lon0, lat0, north = TRUE, ...)
.x |
object of class |
lon0 |
longitude of the rotated pole in degrees |
lat0 |
latitude of the rotated pole in degrees |
north |
logical; if |
... |
ignored |
curvilinear stars object with coordinates in regular long/lat (North pole at lat=90)
if (require("starsdata") && require("maps")) { # data downloaded from https://esgf-data.dkrz.de/search/cosmo-rea/ nc = "netcdf/ts_EUR-6km_ECMWF-ERAINT_REA6_r1i1p1f1_COSMO_v1_mon_201801-201812.nc" f = system.file(nc, package = "starsdata") m = read_mdim(f, "ts") print(m) # NOTE this function is obsolete when reading m as # m = read_mdim(f, "ts", curvilinear = c("longitude", "latitude")) if (require(RNetCDF)) { x = open.nc(f) lon = att.get.nc(x, "rotated_latitude_longitude", "grid_north_pole_longitude") lat = att.get.nc(x, "rotated_latitude_longitude", "grid_north_pole_latitude") close.nc(x) print(c(lon = lon, lat = lat)) } else { lon = -162 lat = 39.25 } m1 = st_rotate(m, lon, lat) print(m1) h = function() maps::map(add = TRUE) plot(m1, downsample = c(10, 10, 5), axes = TRUE, hook = h, mfrow = c(1, 2)) # curvilinear grid: downsample for plotting speed m2 = st_warp(m1, crs = st_crs("OGC:CRS84"), threshold = .1) plot(m2, hook = h, mfrow = c(3, 4)) # regular grid: plots fast }
if (require("starsdata") && require("maps")) { # data downloaded from https://esgf-data.dkrz.de/search/cosmo-rea/ nc = "netcdf/ts_EUR-6km_ECMWF-ERAINT_REA6_r1i1p1f1_COSMO_v1_mon_201801-201812.nc" f = system.file(nc, package = "starsdata") m = read_mdim(f, "ts") print(m) # NOTE this function is obsolete when reading m as # m = read_mdim(f, "ts", curvilinear = c("longitude", "latitude")) if (require(RNetCDF)) { x = open.nc(f) lon = att.get.nc(x, "rotated_latitude_longitude", "grid_north_pole_longitude") lat = att.get.nc(x, "rotated_latitude_longitude", "grid_north_pole_latitude") close.nc(x) print(c(lon = lon, lat = lat)) } else { lon = -162 lat = 39.25 } m1 = st_rotate(m, lon, lat) print(m1) h = function() maps::map(add = TRUE) plot(m1, downsample = c(10, 10, 5), axes = TRUE, hook = h, mfrow = c(1, 2)) # curvilinear grid: downsample for plotting speed m2 = st_warp(m1, crs = st_crs("OGC:CRS84"), threshold = .1) plot(m2, hook = h, mfrow = c(3, 4)) # regular grid: plots fast }
set bounding box parameters of regular grid
st_set_bbox(x, value, ...)
st_set_bbox(x, value, ...)
x |
object of class dimensions, stars or stars_proxy |
value |
object of class bbox |
... |
ignored |
replace POINT simple feature geometry list with an x y raster
st_sfc2xy(x, ...)
st_sfc2xy(x, ...)
x |
object of class |
... |
passed on to as.data.frame.stars |
object of class stars
with a POINT list replaced by x and y raster dimensions. This only works when the points are distributed over a regular or rectilinear grid.
Helper function for specifying the block parameters (nXOff
,
nYOff
, nXsize
, and nYSize
) required by RasterIO
argument in read_stars
st_tile(img_rows, img_cols, x_window, y_window, overlap = 0)
st_tile(img_rows, img_cols, x_window, y_window, overlap = 0)
img_rows |
number of input raster rows (integer) |
img_cols |
number of input raster columns (integer) |
x_window |
number of rows in block (integer) |
y_window |
number of columns in block (integer) |
overlap |
number of overlapping pixels (integer) |
matrix with specified nXOff
, nYOff
, nXsize
,
and nYSize
parameters for every block
## Not run: tif = system.file("tif/L7_ETMs.tif", package = "stars") r = read_stars(tif, proxy = TRUE) tiles = st_tile(nrow(r), ncol(r), 256, 256) for (i in seq_len(nrow(tiles))) { tile = read_stars(tif, proxy = FALSE, RasterIO = tiles[i, ]) # write tiles to separate files write_stars(tile, dsn = paste0(i, ".tif")) } ## End(Not run)
## Not run: tif = system.file("tif/L7_ETMs.tif", package = "stars") r = read_stars(tif, proxy = TRUE) tiles = st_tile(nrow(r), ncol(r), 256, 256) for (i in seq_len(nrow(tiles))) { tile = read_stars(tif, proxy = FALSE, RasterIO = tiles[i, ]) # write tiles to separate files write_stars(tile, dsn = paste0(i, ".tif")) } ## End(Not run)
transform geometries in stars objects to a new coordinate reference system, without warping
## S3 method for class 'stars' st_transform(x, crs, ...) st_transform_proj.stars(x, crs, ...)
## S3 method for class 'stars' st_transform(x, crs, ...) st_transform_proj.stars(x, crs, ...)
x |
object of class |
crs |
object of class |
... |
ignored |
For simple feature dimensions, st_transform is called, leading to lossless transformation. For gridded spatial data, a curvilinear grid with transformed grid cell (centers) is returned, which is also lossless. To convert this to a regular grid in the new CRS
, use st_warp (which is in general lossy).
If array values contain geometries and an array as a whole is of class 'sfc' and has a non-missing CRS, array geometries are also transformed.
geomatrix = system.file("tif/geomatrix.tif", package = "stars") (x = read_stars(geomatrix)) new = st_crs('OGC:CRS84') y = st_transform(x, new) plot(st_transform(st_as_sfc(st_bbox(x)), new), col = NA, border = 'red') plot(st_as_sfc(y, as_points=FALSE), col = NA, border = 'green', axes = TRUE, add = TRUE) image(y, col = heat.colors(12), add = TRUE) plot(st_as_sfc(y, as_points=TRUE), pch=3, cex=.5, col = 'blue', add = TRUE) plot(st_transform(st_as_sfc(x, as_points=FALSE), new), add = TRUE)
geomatrix = system.file("tif/geomatrix.tif", package = "stars") (x = read_stars(geomatrix)) new = st_crs('OGC:CRS84') y = st_transform(x, new) plot(st_transform(st_as_sfc(st_bbox(x)), new), col = NA, border = 'red') plot(st_as_sfc(y, as_points=FALSE), col = NA, border = 'green', axes = TRUE, add = TRUE) image(y, col = heat.colors(12), add = TRUE) plot(st_as_sfc(y, as_points=TRUE), pch=3, cex=.5, col = 'blue', add = TRUE) plot(st_transform(st_as_sfc(x, as_points=FALSE), new), add = TRUE)
Warp (resample) grids in stars objects to a new grid, possibly in an new coordinate reference system
st_warp( src, dest, ..., crs = NA_crs_, cellsize = NA_real_, segments = 100, use_gdal = FALSE, options = character(0), no_data_value = NA_real_, debug = FALSE, method = "near", threshold = NA_real_ )
st_warp( src, dest, ..., crs = NA_crs_, cellsize = NA_real_, segments = 100, use_gdal = FALSE, options = character(0), no_data_value = NA_real_, debug = FALSE, method = "near", threshold = NA_real_ )
src |
object of class |
dest |
object of class |
... |
ignored |
crs |
coordinate reference system for destination grid, only used when |
cellsize |
length 1 or 2 numeric; cellsize in target coordinate reference system units |
segments |
(total) number of segments for segmentizing the bounding box before transforming to the new crs |
use_gdal |
logical; if |
options |
character vector with options, passed on to gdalwarp |
no_data_value |
value used by gdalwarp for no_data (NA) when writing to temporary file;
not setting this when |
debug |
logical; if |
method |
character; see details for options; methods other than |
threshold |
numeric; distance threshold for warping curvilinear grids: new cells at distances larger than threshold are assigned NA values. |
method
should be one of near
, bilinear
, cubic
, cubicspline
, lanczos
, average
, mode
, max
, min
, med
, q1
or q3
; see https://github.com/r-spatial/stars/issues/109
For gridded spatial data (dimensions x
and y
), see figure; the existing grid is transformed into a regular grid defined by dest
, possibly in a new coordinate reference system. If dest
is not specified, but crs
is, the procedure used to choose a target grid is similar to that of projectRaster. This entails: (i) the envelope (bounding box polygon) is transformed into the new crs, possibly after segmentation (red box); (ii) a grid is formed in this new crs, touching the transformed envelope on its East and North side, with (if cellsize is not given) a cellsize similar to the cell size of src
, with an extent that at least covers x
; (iii) for each cell center of this new grid, the matching grid cell of x
is used; if there is no match, an NA
value is used.
geomatrix = system.file("tif/geomatrix.tif", package = "stars") (x = read_stars(geomatrix)) new_crs = st_crs('OGC:CRS84') y = st_warp(x, crs = new_crs) plot(st_transform(st_as_sfc(st_bbox(x)), new_crs), col = NA, border = 'red') plot(st_as_sfc(y, as_points=FALSE), col = NA, border = 'green', axes = TRUE, add = TRUE) image(y, add = TRUE, nbreaks = 6) plot(st_as_sfc(y, as_points=TRUE), pch=3, cex=.5, col = 'blue', add = TRUE) plot(st_transform(st_as_sfc(x, as_points=FALSE), new_crs), add = TRUE) # warp 0-360 raster to -180-180 raster: r = read_stars(system.file("nc/reduced.nc", package = "stars")) r %>% st_set_crs('OGC:CRS84') %>% st_warp(st_as_stars(st_bbox(), dx = 2)) -> s plot(r, axes = TRUE) # no CRS set, so no degree symbols in labels plot(s, axes = TRUE) # downsample raster (90 to 270 m) r = read_stars(system.file("tif/olinda_dem_utm25s.tif", package = "stars")) r270 = st_as_stars(st_bbox(r), dx = 270) r270 = st_warp(r, r270)
geomatrix = system.file("tif/geomatrix.tif", package = "stars") (x = read_stars(geomatrix)) new_crs = st_crs('OGC:CRS84') y = st_warp(x, crs = new_crs) plot(st_transform(st_as_sfc(st_bbox(x)), new_crs), col = NA, border = 'red') plot(st_as_sfc(y, as_points=FALSE), col = NA, border = 'green', axes = TRUE, add = TRUE) image(y, add = TRUE, nbreaks = 6) plot(st_as_sfc(y, as_points=TRUE), pch=3, cex=.5, col = 'blue', add = TRUE) plot(st_transform(st_as_sfc(x, as_points=FALSE), new_crs), add = TRUE) # warp 0-360 raster to -180-180 raster: r = read_stars(system.file("nc/reduced.nc", package = "stars")) r %>% st_set_crs('OGC:CRS84') %>% st_warp(st_as_stars(st_bbox(), dx = 2)) -> s plot(r, axes = TRUE) # no CRS set, so no degree symbols in labels plot(s, axes = TRUE) # downsample raster (90 to 270 m) r = read_stars(system.file("tif/olinda_dem_utm25s.tif", package = "stars")) r270 = st_as_stars(st_bbox(r), dx = 270) r270 = st_warp(r, r270)
replace x y raster dimensions with simple feature geometry list (points, or polygons = rasterize)
st_xy2sfc(x, as_points, ..., na.rm = TRUE)
st_xy2sfc(x, as_points, ..., na.rm = TRUE)
x |
object of class |
as_points |
logical; if |
... |
arguments passed on to |
na.rm |
logical; omit (remove) cells which are entirely missing valued (across other dimensions)? |
object of class stars
with x and y raster dimensions replaced by a single sfc geometry list column containing either points, or polygons. Adjacent cells with identical values are not merged; see st_rasterize
for this.
Sentinel-2 sample tile, downloaded from <https://scihub.copernicus.eu/> reads the four 10-m bands: B2 (490 nm), B3 (560 nm), B4 (665 nm) and B8 (842 nm)
stars_sentinel2
stars_sentinel2
An object of class stars_proxy
(inherits from stars
) of dimension 10980 x 10980 x 4.
subset stars objects
## S3 replacement method for class 'stars_proxy' x[i, downsample = 0] <- value ## S3 method for class 'stars' x[i = TRUE, ..., drop = FALSE, crop = !is_curvilinear(x)] ## S3 replacement method for class 'stars' x[i] <- value st_flip(x, which = 1)
## S3 replacement method for class 'stars_proxy' x[i, downsample = 0] <- value ## S3 method for class 'stars' x[i = TRUE, ..., drop = FALSE, crop = !is_curvilinear(x)] ## S3 replacement method for class 'stars' x[i] <- value st_flip(x, which = 1)
x |
object of class |
i |
first selector: integer, logical or character vector indicating attributes to select, or object of class |
downsample |
downsampling rate used in case |
value |
array of dimensions equal to those in |
... |
further (logical or integer vector) selectors, matched by order, to select on individual dimensions |
drop |
logical; if |
crop |
logical; if |
which |
character or integer; dimension(s) to be flipped |
If i
is an object of class sf
, sfc
or bbox
, the spatial subset covering this geometry is selected, possibly followed by cropping the extent. Array values for which the cell centre is not inside the geometry are assigned NA
. If i
is of class stars
, and attributes of i
are logical
, cells in x
corresponding to NA
or FALSE
cells in i
are assigned an NA
. Dimension ranges containing negative values or NA
may be partially supported.
in an assignment (or replacement form, [<-
), argument i
needs to be either (i) a stars
object with logical attribute(s) that has dimensions matching (possibly after recycling) those of x
, in which case the TRUE
cells will be replaced and i
and/or value
will be recycled to the dimensions of the arrays in x
, or (ii) a length-one integer or character vector indicating which array to replace, in which case value
may be stars object or a vector or array (that will be recycled).
st_flip
flips (reverts) the array values along the chosen dimension
without(s) changing the dimension properties
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) x[,,,1:3] # select bands x[,1:100,100:200,] # select x and y by range x["L7_ETMs.tif"] # select attribute xy = structure(list(x = c(293253.999046018, 296400.196497684), y = c(9113801.64775462, 9111328.49619133)), .Names = c("x", "y")) pts = st_as_sf(data.frame(do.call(cbind, xy)), coords = c("x", "y"), crs = st_crs(x)) image(x, axes = TRUE) plot(st_as_sfc(st_bbox(pts)), col = NA, add = TRUE) bb = st_bbox(pts) (xx = x[bb]) image(xx) plot(st_as_sfc(bb), add = TRUE, col = NA) image(x) pt = st_point(c(x = 290462.103109179, y = 9114202.32594085)) buf = st_buffer(st_sfc(pt, crs = st_crs(x)), 1500) plot(buf, add = TRUE) buf = st_sfc(st_polygon(list(st_buffer(pt, 1500)[[1]], st_buffer(pt, 1000)[[1]])), crs = st_crs(x)) image(x[buf]) plot(buf, add = TRUE, col = NA) image(x[buf, crop=FALSE]) plot(buf, add = TRUE, col = NA) # with i of class stars: x[x > 75] # generates lots of NA's; pattern for each band x[x[,,,1] > 75] # recycles a single band template for all bands x = read_stars(tif) # replace, using a logical stars selector: cuts all values above 90 to 90 x[x > 90] = 90 # replace a single attribute when there are more than one: s = split(x) names(s) = paste0("band", 1:6) # rescale only band 1: s[1] = s[1] * 0.75 # rescale only attribute named "band2": s["band2"] = s["band2"] * 0.85 # create a new attribute from a numeric vector: s["rnorm"] = rnorm(prod(dim(s))) s lc = read_stars(system.file("tif/lc.tif", package = "stars")) x = c(orig = lc, flip_x = st_flip(lc, "x"), flip_y = st_flip(lc, "y"), flip_xy = st_flip(lc, c("x", "y")), along = 3) plot(x)
tif = system.file("tif/L7_ETMs.tif", package = "stars") x = read_stars(tif) x[,,,1:3] # select bands x[,1:100,100:200,] # select x and y by range x["L7_ETMs.tif"] # select attribute xy = structure(list(x = c(293253.999046018, 296400.196497684), y = c(9113801.64775462, 9111328.49619133)), .Names = c("x", "y")) pts = st_as_sf(data.frame(do.call(cbind, xy)), coords = c("x", "y"), crs = st_crs(x)) image(x, axes = TRUE) plot(st_as_sfc(st_bbox(pts)), col = NA, add = TRUE) bb = st_bbox(pts) (xx = x[bb]) image(xx) plot(st_as_sfc(bb), add = TRUE, col = NA) image(x) pt = st_point(c(x = 290462.103109179, y = 9114202.32594085)) buf = st_buffer(st_sfc(pt, crs = st_crs(x)), 1500) plot(buf, add = TRUE) buf = st_sfc(st_polygon(list(st_buffer(pt, 1500)[[1]], st_buffer(pt, 1000)[[1]])), crs = st_crs(x)) image(x[buf]) plot(buf, add = TRUE, col = NA) image(x[buf, crop=FALSE]) plot(buf, add = TRUE, col = NA) # with i of class stars: x[x > 75] # generates lots of NA's; pattern for each band x[x[,,,1] > 75] # recycles a single band template for all bands x = read_stars(tif) # replace, using a logical stars selector: cuts all values above 90 to 90 x[x > 90] = 90 # replace a single attribute when there are more than one: s = split(x) names(s) = paste0("band", 1:6) # rescale only band 1: s[1] = s[1] * 0.75 # rescale only attribute named "band2": s["band2"] = s["band2"] * 0.85 # create a new attribute from a numeric vector: s["rnorm"] = rnorm(prod(dim(s))) s lc = read_stars(system.file("tif/lc.tif", package = "stars")) x = c(orig = lc, flip_x = st_flip(lc, "x"), flip_y = st_flip(lc, "y"), flip_xy = st_flip(lc, c("x", "y")), along = 3) plot(x)
write stars object to gdal dataset (typically: to file)
write_stars(obj, dsn, layer, ...) ## S3 method for class 'stars' write_stars( obj, dsn, layer = 1, ..., driver = detect.driver(dsn), options = character(0), type = if (is.factor(obj[[1]]) && length(levels(obj[[1]])) < 256) "Byte" else "Float32", NA_value = NA_real_, update = FALSE, normalize_path = TRUE, scale_offset = c(1, 0) ) ## S3 method for class 'stars_proxy' write_stars( obj, dsn, layer = 1, ..., driver = detect.driver(dsn), options = character(0), scale_offset = c(1, 0), type = "Float32", NA_value = NA_real_, chunk_size = c(dim(obj)[1], floor(2.5e+07/dim(obj)[1])), progress = TRUE ) detect.driver(filename)
write_stars(obj, dsn, layer, ...) ## S3 method for class 'stars' write_stars( obj, dsn, layer = 1, ..., driver = detect.driver(dsn), options = character(0), type = if (is.factor(obj[[1]]) && length(levels(obj[[1]])) < 256) "Byte" else "Float32", NA_value = NA_real_, update = FALSE, normalize_path = TRUE, scale_offset = c(1, 0) ) ## S3 method for class 'stars_proxy' write_stars( obj, dsn, layer = 1, ..., driver = detect.driver(dsn), options = character(0), scale_offset = c(1, 0), type = "Float32", NA_value = NA_real_, chunk_size = c(dim(obj)[1], floor(2.5e+07/dim(obj)[1])), progress = TRUE ) detect.driver(filename)
obj |
object of class |
dsn |
gdal dataset (file) name |
layer |
attribute name; if missing, the first attribute is written |
... |
passed on to gdal_write |
driver |
driver driver name; see st_drivers |
options |
character vector with dataset creation options, passed on to GDAL |
type |
character; output binary type, one of: |
NA_value |
non-NA value that should represent R's |
update |
logical; if |
normalize_path |
logical; see read_stars |
scale_offset |
length 2 numeric vector with scale, offset values: raw values computed by raw = (value - offset) / scale are written to dsn; scale and offset values are written to dsn or else a warning is raised |
chunk_size |
length two integer vector with the number of pixels (x, y) used in the read/write loop; see details. |
progress |
logical; if |
filename |
character; used for guessing driver short name based on file extension; see examples |
write_stars
first creates the target file, then updates it sequentially by writing blocks of chunk_size
.
in case obj
is a multi-file stars_proxy
object, all files are written as layers into the output file dsn
detect.driver("L7_ETMs.tif")
detect.driver("L7_ETMs.tif")