[R-sig-Geo] Convert rectilinear to regular grid in R (stars and raster)

[Frederico Faleiro]

Dear all,

I am trying convert netCDF data from rectilinear to regular grid in R, but
I am not sure what is the best approach for it.
I would like to convert the grid modifying the original resolution only
when necessary. In some cases I must change resolution to keep regular cell
size and cover all the world extent.
I have tried without success st_warp from stars and I find a way with
rasterize from raster package (see script bellow).

I have the following questions:
1. What am I doing wrong in stars package (see script bellow)?
2. Is there a way to rasterize points to raster using metrics that consider
distance (e.g. bilinear, nearest neighbor, etc)? Note that I need
use points in raster package because it do not handle with rectilinear grid.
3. Do you recommend any other approach?

The data are available in:
GDrive:
https://drive.google.com/file/d/1HKxFpAwvY9k_cFbasagKHAwi1luzLv90/view?usp=sharing
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

# script
library(stars)
library(raster)

# prepare data
bc <- read_ncdf("pr_Amon_bcc-csm1-1_rcp45_r1i1p1_200601-209912.nc")
# get coordinates and data from the first time period
bc1 <- as.data.frame(bc[, , , 1])
# convert longitude to variate between -180 and 180
bc1$lon <- ifelse(bc1$lon > 180, bc1$lon - 360, bc1$lon)

# original resolution
bc1$lon %>% diff %>% table
bc1$lat %>% diff %>% table
# new resolution
new_res <- c(2.8125, 2.8125)
# Is it multiple of 180 or 360?
180 %% new_res

# create a new standard grid
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
# stars package: use st_warp to change from rectilinear to regular raster
with different methods
bc_reg <- bc %>% st_warp(grid_st, method = "average")
bc_reg <- bc %>% st_warp(grid_st, method = "near")
bc_reg <- bc %>% st_warp(grid_st, method = "bilinear")
# In all cases I have this error: Error in 1:prod(dims[dxy]) : NA/NaN
argument

# raster package: rasterize the points values to the new regular grid
bc1_sp <- bc1
coordinates(bc1_sp) <- ~ lon + lat
bc1_reg <- rasterize(bc1_sp, grid_r, field = "pr", fun = mean)

# plot the first time period
dev.new(title = "rectilinear")
plot(bc[, , , 1])
dev.new(title = "regular with raster package")
plot(bc1_mn_reg)

Best regards,

Frederico Faleiro
Postdoctoral Researcher in the INCT-EECBio (https://www.eecbio.ufg.br/)
Federal University of Goiás | Brazil
RG: https://www.researchgate.net/profile/Frederico_Faleiro

	[[alternative HTML version deleted]]