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

[Micha Silver]

On 03/10/2020 16:32, Vijay Lulla wrote:
> 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 <http://is.na>(s)) }
> st_apply(m_stars, MARGIN="DOY", FUN=cnt_pixels1)$cnt_pixels1


Excellent, thanks!



>
> 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 
> <mailto: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-665918
>     https://orcid.org/0000-0002-1128-1325
>
>     _______________________________________________
>     R-sig-Geo mailing list
>     R-sig-Geo using r-project.org <mailto: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>
>
-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918
https://orcid.org/0000-0002-1128-1325