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

[Frederico Faleiro]

Hi Edzer,

thank you for your reply and modifications in the function.
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,
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


Em ter., 24 de mar. de 2020 às 11:59, Edzer Pebesma <
edzer.pebesma using uni-muenster.de> escreveu:

> 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
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]