Convert rectilinear to regular grid in R (stars and raster)

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

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]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Convert rectilinear to regular grid in R (stars and raster)

edzer
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
> [hidden email]
> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

pEpkey.asc (3K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: 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 <
[hidden email]> 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
> > [hidden email]
> > 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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo