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

[Edzer Pebesma]

The stars part of this question has also been sent as an issue on
github, and that is where it has been answered (partially solved):

https://github.com/r-spatial/stars/issues/264

please feel free to send follow-up's here, or there.

On 3/23/20 6:45 PM, Frederico Faleiro wrote:
> 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]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: pEpkey.asc
Type: application/pgp-keys
Size: 3110 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20200324/c58f7c25/attachment.bin>