Error in st_warp convert NetCDF from rectilinear to regular grid

[FVFaleiro]

I am trying convert a NetCDF from rectilinear to regular grid with st_warp without success (see bellow). What am I doing wrong?

The data are available in:
GDrive: https://drive.google.com/file/d/1HKxFpAwvY9k_cFbasagKHAwi1luzLv90/view?usp=sharing
or
CMIP5 portal (you must make a login): http://aims3.llnl.gov/thredds/fileServer/cmip5_css01_data/cmip5/output1/BCC/bcc-csm1-1/rcp45/mon/atmos/Amon/r1i1p1/v20120705/pr/pr_Amon_bcc-csm1-1_rcp45_r1i1p1_200601-209912.nc

library(stars)
library(raster)
bc <- read_ncdf("pr_Amon_bcc-csm1-1_rcp45_r1i1p1_200601-209912.nc")
# create a new standard grid
new_res <- c(2.8125, 2.8125)
wgs84 <- "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"
grid_r <- raster(res = new_res, crs = wgs84)
grid_st <- st_as_stars(grid_r)
# change from rectilinear to regular grid
bc_reg <- bc %>% st_warp(grid_st, method = "average")
# Error in 1:prod(dims[dxy]) : NA/NaN argument

[edzer]

Cross post from https://stat.ethz.ch/pipermail/r-sig-geo/2020-March/027997.html

[edzer]

This makes the following work:

library(stars)
library(raster)

# prepare data
bc <- read_ncdf("pr_Amon_bcc-csm1-1_rcp45_r1i1p1_200601-209912.nc")

# create a new standard grid, with lat 0...360
grid_st <- st_as_stars(st_bbox(c(xmin=0,xmax=360,ymin=-90,ymax=90))) %>% st_set_crs(4326)

# change from rectilinear to regular grid
# stars package: use st_warp to change from rectilinear to regular raster with different methods
(bc_reg <- bc %>% st_warp(grid_st, method = "near"))
plot(bc_reg[,,,1:12])

Note that the new raster also needs to cover latitude 0-360, as the old one, that is something that also needs fixing.

As the documentation says, the other method options are for use_gdal=TRUE, which means do this through GDAL, which means no rectilinear grids supported to start with. One alternative would be to do this using polygons (st_as_sf), and the using sf::st_interpolate_aw.

[FVFaleiro]

Hi Edzer, I answer at r-sig-geo list (http://r-sig-geo.2731867.n2.nabble.com/Convert-rectilinear-to-regular-grid-in-R-stars-and-raster-td7593462.html#a7593477), but I think you did not see it.

Now the code is working fine, but the structure of the output was changed. When I try extract the values converting to data.frame all values was replaced by NA. However the plot is OK as you have tested.

1. How can I create a new grid with the new resolution in stars? Am I doing it correctly (see bellow)?
2. How can I get the data.frame with the new values?
3. I have tried the function sf::st_interpolate_aw (see bellow), but there is some warnings that I do not know the consequences for the calculation.

# update packages
library(devtools)
install_github("r-spatial/sf")
install_github("r-spatial/stars")

library(stars)
# read file
bc <- read_ncdf("pr_Amon_bcc-csm1-1_rcp45_r1i1p1_200601-209912.nc")
# create the new grid
grid_st <- st_as_stars(st_bbox(c(xmin=0,xmax=360,ymin=-90,ymax=90))) %>% st_set_crs(4326)
# change resolution. Is it right?
attr(grid_st, "dimensions")[[1]]$delta <- 2.8125
attr(grid_st, "dimensions")[[2]]$delta <- - 2.8125

# change from rectilinear to regular grid with stars
bc_reg <- bc %>% st_warp(grid_st, method = "near")
# plot
plot(bc_reg[, , , 1:12])

# convert to data.frame
bc_reg_df <- as.data.frame(bc_reg_st)
# there is only NA in the variable
str(bc_reg_df)
range(bc_reg_df$pr)

# change from rectilinear to regular grid with sf
grid_sf <- st_as_sf(grid_st)
plot(grid_sf)
bc_sf <- st_as_sf(bc)
plot(bc_sf)
bc_sf_reg <- st_interpolate_aw(bc_sf, grid_sf, F)
#although coordinates are longitude/latitude, st_intersection assumes that they are planar
#Warning message:
#In st_interpolate_aw.sf(bc_sf, grid_sf, F) :
#  st_interpolate_aw assumes attributes are constant over areas of x

Best regards.

[edzer]

This is caused by changing both delta values without changing the number of rows/cols; now you have a grid that extends beyond the Earth (meaningless); read into st_as_stars.bbox to see how you can create rasters properly.

[FVFaleiro]

@edzer I did it with other grids, but I have the same problem (i.e. the variable have only NA values).

library(stars)
# read file
bc <- read_ncdf("pr_Amon_bcc-csm1-1_rcp45_r1i1p1_200601-209912.nc")
# create the new grid
grid_st <- st_as_stars(st_bbox(c(xmin = 0, xmax = 360, ymin = -90, ymax = 90)), dx = 2.8125, dy = 2.8125) %>% st_set_crs(4326)
# change from rectilinear to regular grid with stars
bc_st_reg <- bc %>% st_warp(grid_st, method = "near")
# convert to data.frame
bc_st_reg_df <- as.data.frame(bc_reg_st)
# the variable have only NA before and after transform to data.frame, despite I can plot it
summary(bc_st_reg$pr)
summary(bc_st_reg_df$pr)
plot(bc_st_reg[, , , 1:12])

Best regards.

[edzer]

Your script has some typos, but the plot looks good to me. Which versions of stars and sf are you using?

[FVFaleiro]

Sorry about the typos, now it is correct. The version of stars is 0.4-1 installed from github.

library(stars)
# read file
bc <- read_ncdf("pr_Amon_bcc-csm1-1_rcp45_r1i1p1_200601-209912.nc")
# create the new grid
grid_st <- st_as_stars(st_bbox(c(xmin = 0, xmax = 360, ymin = -90, ymax = 90)), dx = 2.8125, dy = 2.8125) %>% st_set_crs(4326)
# change from rectilinear to regular grid with stars
bc_st_reg <- bc %>% st_warp(grid_st, method = "near")
# convert to data.frame
bc_st_reg_df <- as.data.frame(bc_st_reg)
# the variable have only NA before and after transform to data.frame, despite I can plot it
summary(bc_st_reg$pr)
summary(bc_st_reg_df$pr)
plot(bc_st_reg[, , , 1:12])