An HTML attachment was scrubbed... URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20201003/3a263faf/attachment.html>
Count of Non-NA pixels in stars object
3 messages · Vijay Lulla, Micha Silver
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 at 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 at 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]]
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 at gmail.com
<mailto:tsvibar at 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 at r-project.org <mailto:R-sig-Geo at 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