Regridding High-resolution Data
Often time, we use observations to validate model simulations. In addition to data quality, a big challenge is that we need to regrid the observations (or model results) to match with its counterpart.
Here we are using the raster package to perform a bilinear
interpolation on the 5-min-by-5-min data of global maize production.
# ====== loading libraries =====
library(ncdf4) # for reading nc files
library(raster) # for resampling
## Loading required package: sp
library(reshape2) # for reshaping files
library(tidyverse) # for tidying up data
## ── Attaching packages ─────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.3 ✓ dplyr 1.0.1
## ✓ tidyr 1.1.1 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks raster::select()
The dataset we use here is in NetCDF format and can be downloaded here.
# ====== reading data =====
nc.name = "maize_AreaYieldProduction.nc"
# open the NC file
nc = nc_open(filename = nc.name)
# extract coordinate and value of interest
lon.old = ncvar_get(nc = nc, varid = "longitude")
lat.old = ncvar_get(nc = nc, varid = "latitude")
old.matrix = ncvar_get(nc = nc, varid = "maizeData")[,,6] # <--- picking the 6th layer of data = production of maize. check out the .nc file for details
nc_close(nc)
rm(nc)
# Different data could have different methods to assign values to the matrix (increasing vs decreasing lon)
# I found that the best way to avoid dealing with that is to transform the matrix into a named matrix and then a dataframe, i.e. a list of (lon, lat, and value)
dimnames(old.matrix) = list(lon = lon.old, lat = lat.old)
old.df = melt(old.matrix)
head(old.df)
## lon lat value
## 1 -179.9583 89.95833 NaN
## 2 -179.8750 89.95833 NaN
## 3 -179.7917 89.95833 NaN
## 4 -179.7083 89.95833 NaN
## 5 -179.6250 89.95833 NaN
## 6 -179.5417 89.95833 NaN
Now we have a dataframe storing the values. We first convert the dataframe into a raster object so that we can apply the resample function in the raster package. Then, we define the solution of the target grid size (usually matching with model).
# ===== regridding =====
# make raster for the old gridded data
r = rasterFromXYZ(old.df)
# define resolution of the new grid
dlat.new = 1.9 # new delta lat
nlat.new = 96 # number of new lat
dlon.new = 2.5 # new delta lon
nlon.new = 144 # number of new lon
# make raster for the new grid
s = raster(nrow = nlat.new, ncol = nlon.new)
# use resample to regridded data
t = resample(r, s, method = "bilinear")
Finally, we can convert the raster t back to a matrix or a
dataframe for further processing.
# converting the regridded raster back to a named matrix
new.matrix = as.matrix(t)
dimnames(new.matrix) = list(lat = seq(90, -90, length.out = nlat.new),
lon = seq(-177.5, 180, by = dlon.new)) # <--- this setting is a bit tricky as lon of -180º = 180º
# or, further, into a dataframe
new.df = melt(new.matrix)
head(new.df)
## lat lon value
## 1 90.00000 -177.5 NA
## 2 88.10526 -177.5 NA
## 3 86.21053 -177.5 NA
## 4 84.31579 -177.5 NA
## 5 82.42105 -177.5 NA
## 6 80.52632 -177.5 NA
Finally, we can compare the original data and the regridded data visually.
This is what the original data looks like:
# plot the original data to check
old.df %>% ggplot(aes(x = lon, y = lat, fill = value)) +
geom_raster(interpolate = FALSE) + # adding heat maps
scale_fill_viridis_b(na.value = NA ) + # change color
borders() + # adding country borders
coord_equal(expand = FALSE) # keeping a nice aspect ratio
## Warning: Removed 7123007 rows containing missing values (geom_raster).
# then the regrided data
new.df %>% ggplot(aes(x = lon, y = lat, fill = value)) +
geom_raster(interpolate = FALSE) + # adding heat maps
scale_fill_viridis_b(na.value = NA ) + # change color
borders() + # adding country borders
coord_equal(expand = FALSE) # keeping a nice aspect ratio
## Warning: Removed 9010 rows containing missing values (geom_raster).
Done!
