[R-sig-Geo] Count of Non-NA pixels in stars object

[Vijay Lulla]

Micha,
I am not sure why the function works individually but not on the stars
object. IMO, a much simpler function, at least for st_apply, of what you
are trying to do might be:

cnt_pixels1 <- function(s) { sum(!is.na(s)) }
st_apply(m_stars, MARGIN="DOY", FUN=cnt_pixels1)$cnt_pixels1

And if you wish to call this individually you will have to do
cnt_pixels1(as.array(m_stars[,,,1]$Value)) .

HTH,
Vijay.

On Sat, Oct 3, 2020 at 6:29 AM Micha Silver <tsvibar using gmail.com> wrote:

> 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-665918https://orcid.org/0000-0002-1128-1325
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>


-- 
Vijay Lulla, PhD
ORCID | <https://orcid.org/0000-0002-0823-2522> Homepage
<http://vlulla.github.io> | Google Scholar
<https://scholar.google.com/citations?user=VjhJWOgAAAAJ&hl=en> | Github
<https://github.com/vlulla>

	[[alternative HTML version deleted]]