|
Hello:
I'm trying to get the number of pixels with non-NA values in the
'Value' attrib. of a stars object.
My first try was to read a list of geotiff files into raster
objects, then use this function to get the number of pixels:
c = raster::extract(r, site,fun=function(x, ...)
length(na.omit(x)))
This worked OK within sapply(), but was rather slow.
Instead I read the files into a stars object and I'm trying to
switch to st_apply() to gain efficiency. Here's what I've tried
(my reprex):
library(stars)
# An example stars structure
m_stars = structure(list(Value = structure(c(NA, 13458, 13315,
13381, NA,
13483, 13400, 13251,
13251, 13282, NA, NA, 13663, 13663, 13174,
NA, 13783, 13754, NA,
NA, 13664, 13643, 13800, 13800, 13797,
NA, NA, 13796, 13796,
NA, NA, 13515, 13574, 13539, NA, 13383,
13407, 13407, 13407,
13541, NA, NA, 13399, 13399, 13399, NA,
NA, NA, NA, NA,
12402, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
13761, 13662, 13662,
NA, 13651, 13627, 13627, 13627, 13647, NA,
NA, 13607, 13607,
NA), .Dim = c(x = 5L, y = 3L, DOY = 5L))), dimensions =
structure(list(
x =
structure(list(from = 1L, to = 5L, offset = -3.8610742258177,
delta = 0.01151941823202, refsys = structure(list(input =
"GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS
84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]",
wkt = "GEOGCS[\"WGS 84\",\n DATUM[\"WGS_1984\",\n
SPHEROID[\"WGS 84\",6378137,298.257223563,\n
AUTHORITY[\"EPSG\",\"7030\"]],\n AUTHORITY[\"EPSG\",\"6326\"]],\n
PRIMEM[\"Greenwich\",0],\n UNIT[\"degree\",0.0174532925199433],\n
AUTHORITY[\"EPSG\",\"4326\"]]"), class = "crs"),
point = FALSE, values = NULL), class = "dimension"),
y =
structure(list(from = 1L, to = 3L, offset = 57.118130712839,
delta = -0.0141360917583313, refsys = structure(list(
input = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS
84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]]",
wkt = "GEOGCS[\"WGS 84\",\n DATUM[\"WGS_1984\",\n
SPHEROID[\"WGS 84\",6378137,298.257223563,\n
AUTHORITY[\"EPSG\",\"7030\"]],\n AUTHORITY[\"EPSG\",\"6326\"]],\n
PRIMEM[\"Greenwich\",0],\n UNIT[\"degree\",0.0174532925199433],\n
AUTHORITY[\"EPSG\",\"4326\"]]"), class = "crs"),
point = FALSE, values = NULL), class = "dimension"),
DOY =
structure(list(from = 1L, to = 5L, offset = NA_real_,
delta = NA_real_, refsys = NA_character_, point = NA,
values = c("DOY_2002_001", "DOY_2002_009", "DOY_2002_017",
"DOY_2002_025", "DOY_2002_033")), class = "dimension")), raster =
structure(list(
affine = c(0, 0), dimensions = c("x", "y"), curvilinear = FALSE),
class = "stars_raster"), class = "dimensions"), class = "stars")
# My count function
cnt_pixels = function(s, na.rm = TRUE, ...) {
# Count number of non NA values in $Value attrib
sdf = as.data.frame(s)
if (na.rm) {
sdf = sdf[complete.cases(sdf),]
}
return(length(sdf$Value))
}
# Using st_apply
cnt = st_apply(m_stars,
MARGIN = "DOY",
FUN = cnt_pixels,
na.rm = TRUE)$cnt_pixels
# The function alone seems to work...
(cnt_pixels(m_stars[,,,1]))
(cnt_pixels(m_stars[,,,4]))
# But run thru st_apply...
(cnt) # shows all zeros, ??
Why am I getting all zeros in the resulting "cnt" list?
Am I going about this the right way? Is there some built in method
that I missed?
Thanks, Micha
--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
https://orcid.org/0000-0002-1128-1325
|
