[R-sig-Geo] Issue with st_warp, gdal and NAs

[Edzer Pebesma]

Dear Julian, you can do that by passing a non-NA value as no_data_value, 
e.g. by

ag3 <- st_warp(src = s1,
               dest = st_as_stars(st_bbox(s1), dx = 1, dy = 1),
               use_gdal = TRUE,
               method = "average",
               no_data_value = -9999)


On 16/07/2021 12:17, Julian M. Burgos wrote:
> Dear list,
> 
> My apologies for posting this again.  I am not sure if it went through the first time.
> 
> I am having a bit of trouble using st_warp() to aggregate star objects with the option "use_gdal = TRUE".  I am using the "average" method, and I want cells with NA values to be ignored when computing the averages.  To demonstrate, I will make a small raster for testing, and will include cells with NA values:
> 
> #-------------------------------------------------
>    
> library(raster)
> 
> # Create raster for testing
> myfile <- "~/myraster.tif"
> myraster <- raster(xmn = 0, ymn = 0, xmx = 10, ymx = 10, resolution = 0.1,
>                     crs = 4326, vals = sample(1:20, size = 10000, replace = TRUE))
> myraster[sample(1:10000, 200)] <- NA
> writeRaster(myraster, filename = myfile, overwrite = TRUE)
> 
> #-------------------------------------------------
> 
> With the raster package, I can aggregate the cells ignoring the NA values like this:
> 
> 
> r1 <- raster(myfile)
> ag1 <- raster::aggregate(x = raster(myfile), fact = 10,
>                           fun = mean, na.rm = TRUE)
> 
> ...which is what I want.  I can also aggregate the cells including the NAs, like this:
>                                                                                     
> ag2 <- raster::aggregate(x = raster(myfile), fact = 10,
>                           fun = mean, na.rm = FALSE)
> 
> .... which produces a raster with a lots of gaps.  Now, to do a similar aggregation with the stars package I can do this:
> 
> s1 <-  read_stars(myfile)
> 
> ag3 <- st_warp(src = s1,
>                dest = st_as_stars(st_bbox(s1), dx = 1, dy = 1),
>                use_gdal = TRUE,
>                method = "average")
> 
> The resulting raster (ag3) has gaps like ag2, which means that the NAs are not being ignored when computing the averages.  With use_gdal_TRUE and method="average" we are using gdalwarp with the average resampling method, which according to the documentation "computes the weighted average of all non-NODATA contributing pixels.".  So I guess the problem is that the NAs are not being recognized as NODATA pixels.
> 
> Does anybody know how to solve this?
> 
> Many thanks,
> 
> Julian
> 

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081