Regridding of irregular model grids
by Mark R Payne
DTU-Aqua, Charlottenlund, Denmark
http://www.staff.dtu.dk/mpay
Wed Jun 10 02:45:27 2015
Many (most) oceanographic and climate models work on grids that are not based on standard lat-long coordinates. While it is common for these models to provide remapping and gridded outputs to make the user’s life a bit easier, it is not by any means unusual to have to do the regridding yourself. This demo shows how it can be done using the standard tools in R. We use outputs from the Max Planck Ocean Model (MPI-OM) to illustrate how this is done.
Setup
The MPI-OM model provides its outputs as NetCDF files, so its easier to work directly with them using the raster package. The NetCDF files store data in 4D (time, depth, x, y) on the native model grid. We can get a good overview of the structure of this file using the netcdf4 package - alternatively, we can also use the “ncdump” tool if we’ere working in linux.
# Get NetCDF meta information using the netcdf4 package
library(ncdf4)
dat.fname <- "data//MPI_temp.nc"
ncid <- nc_open(dat.fname)
ncid## File data//MPI_temp.nc (NC_FORMAT_CLASSIC):
##
## 5 variables (excluding dimension variables):
## float lat[x,y]
## standard_name: latitude
## long_name: latitude
## units: degrees
## _CoordinateAxisType: Lat
## bounds: lat_bnds
## float lat_bnds[nv4,x,y]
## float lon[x,y]
## standard_name: longitude
## long_name: longitude
## units: degrees
## _CoordinateAxisType: Lon
## bounds: lon_bnds
## float lon_bnds[nv4,x,y]
## float var2[x,y,lev,time]
## coordinates: lon lat
## _FillValue: -8.99999987309029e+33
## missing_value: -8.99999987309029e+33
##
## 5 dimensions:
## y Size:220
## x Size:254
## nv4 Size:4
## lev Size:40
## long_name: generic
## units: level
## axis: Z
## time Size:1 *** is unlimited ***
## standard_name: time
## units: day as %Y%m%d.%f
## calendar: proleptic_gregorian
##
## 5 global attributes:
## CDI: Climate Data Interface version 1.6.2 (http://code.zmaw.de/projects/cdi)
## Conventions: CF-1.4
## history: Wed Dec 2 15:50:36 2015: ncks -d time,1 ncep_07c_tho_all_1948-2010.mon.nc MPI_temp.nc
## Thu May 01 09:39:41 2014: cdo -a -f nc setgrid,../GR15s.nc -selindexbox,2,255,1,220 -setgrid,r256x220 ncep_07c_tho_all_1948-2010.mon.ext4 ncep_07c_tho_all_1948-2010.mon.nc
## CDO: Climate Data Operators version 1.6.2 (http://code.zmaw.de/projects/cdo)
## NCO: 4.4.2Here, “var2” is the variable of interest (in this case temperature). However, the variables “lon” and “lat” are also the key, as they give the longitude and latitude at the locations on the model grid. We will use these to develop a remapping between the model grid, (x,y) and geographical space (lon, lat).
We now setup the appropriate raster objects as follows. Note that it is good practice to specify the exact variable. Bricks are also limited to 3D, so we also need to specify the depth level that we are interested in - in this case, level 1, which corresponds to SST.
library(raster)
b <- brick(dat.fname,varname="var2",level=1)
b## class : RasterBrick
## dimensions : 220, 254, 55880, 1 (nrow, ncol, ncell, nlayers)
## resolution : 1, 1 (x, y)
## extent : 0.5, 254.5, 0.5, 220.5 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : /home/mpayne/dtu/Lab/2015_06_10_Remapping_of_model_outputs/data/MPI_temp.nc
## names : X19480229
## z-value : 19480229
## varname : var2
## level : 1Lets visualise how the first value looks
It kind of looks like the world, but there is also some serious distortion going on - look at the size of Iceland for example versus Africa and South America.
Now we create the lookup table by extracting the lon-lat coordinates at each model grid location.
lon <- raster(dat.fname,varname="lon")
lat <- raster(dat.fname,varname="lat")
#We now convert these to a data.frame
library(reshape)
remap.tbl <- data.frame(coordinates(lon),
lon=as.vector(lon),lat=as.vector(lat))
tail(remap.tbl)## x y lon lat
## 55875 249 1 315.0766 76.58982
## 55876 250 1 314.5945 76.55020
## 55877 251 1 314.1198 76.50690
## 55878 252 1 313.6530 76.45997
## 55879 253 1 313.1947 76.40948
## 55880 254 1 312.7453 76.35550Note the we have longitudes of 315 here! It is common in climate models for langitude to run from 0 -360, instead of -180-180, which is more common. We can correct for this as follows
remap.tbl$lon <- ifelse(remap.tbl$lon>180,
remap.tbl$lon-360,
remap.tbl$lon)
range(remap.tbl$lon)## [1] -179.9865 179.9932We can get a good idea of how the grid actually looks like by then plotting the grid points on a geographical mop
library(maps)
map("world",fill=TRUE,mar=c(0,0,0,0))
points(lat ~lon,remap.tbl,pch=".",col="red")
Note that the North pole here has been shifted so that it lies over Greenland. This has the advantage of giving high resolution in the North Atlantic, which is one of the strengths of this model.
Extract values from model
The most common task when using such model data is to extract the values at given points. However, this always raises a question - the points will never fall prefectly on the grid and interpolation will be required, so where do you do this - in model space or geographical space? The answer will depend on the exact situation. Some variables are highly sensitive to this choice e.g. when calculating volumetric fluxs through a section. Fortunately, scalar variables such as SST and salinity are generally not so sensitive. From a computational point of view, it is almost always more efficient to move the extraction points to the model grid, rather than vice-versa.
To illustrate this, lets say we want to take a transect of temperatures across 65 N, from Greenland to Norway. e.g.
library(mapdata)
map("worldHires",xlim=c(-50,15),ylim=c(50,80),fill=TRUE,mar=c(0,0,0,0))
lon.pts <- seq(-50,15,by=2)
lat.pts <- rep(65,length(lon.pts))
points(lon.pts,lat.pts,pch=4,col="red")
We can then use the akima package to do the spatial interpolation from geographic space (lon-lat) to model space, using the remap.tbl that we created earlier. Akima has a number of useful tools for doing this type of analysis. To make the script run a little quicker, we restrict the remap table to the North Atlantic first.
library(akima)
remap.NA <- subset(remap.tbl,lon>-60 & lon<30 & lat > 30)
x.pts <- interpp(remap.NA$lon,remap.NA$lat,remap.NA$x,
xo=lon.pts,yo=lat.pts)
y.pts <- interpp(remap.NA$lon,remap.NA$lat,remap.NA$y,
xo=lon.pts,yo=lat.pts)This can be understood as follows. We have a set of points characterised by their lon-lat (not necessarily on a regular grid) for each of them we have associated a value of some variable - in this case it is one of x or y, the coordinates in the model grid. We then use interpp() to interpolate this on this irregularly spaced set of points to find the value of x (or y) at each of the lookup points. By doing so we have therefore transformed the lookup points from geographic coordinates to model grid coordinates.
We can then visualise our line in model space
pts <- data.frame(lon=lon.pts,lat=lat.pts,x=x.pts$z,y=y.pts$z)
plot(b[[1]])
points(y~x,pts,col="red",pch=4)
Note the distortion in what is a straight line in geographical coordinates!
Finally, we can lookup the individual values in the standard manner with raster.
pts.sp <- pts
coordinates(pts.sp) <- ~ x + y
pts$temp <- extract(b[[1]],pts.sp,method="bilinear")
plot(temp ~ lon,pts,type="b")
rug(pts$lon)
Note the gap in the middle associated with Iceland.
Remap the entire field
Alternatively, there may be instances (e.g. publishing an article), where you need to remap the entire field. This can be done in a similar manner, again using the remap table that we defined. However, akima has some useful tools again that make it go a bit quicker.
First, lets define the raster that we want to remap to.
geo.r <- raster(extent(-50,0,50,75))
res(geo.r) <- c(0.25,0.25)Then get the corresponding x and y coordinates using interp(). We could alternatively use interpp() but interp() should be quicker.
r.coords.x <- interp(remap.NA$lon,remap.NA$lat,remap.NA$x,
xo=xFromCol(geo.r),yo=yFromRow(geo.r))
r.coords.y <- interp(remap.NA$lon,remap.NA$lat,remap.NA$y,
xo=xFromCol(geo.r),yo=yFromRow(geo.r))
r.coords <- expand.grid(lon=xFromCol(geo.r),
lat=yFromRow(geo.r))
r.coords$x <- as.vector(r.coords.x$z)
r.coords$y <- as.vector(r.coords.y$z)Then extract as previously, form a raster, and plot!
r.coords.sp <- r.coords
coordinates(r.coords.sp) <- ~x +y
r.coords$temp <- extract(b[[1]],r.coords.sp,method="bilinear")
geo.r <- setValues(geo.r,r.coords$temp)
plot(geo.r)
map("worldHires",add=TRUE)
Looking pretty good!
Bonus info: In New Zealand, when I was growing up, Santa Claus lived at the North Pole. However, when I moved to Denmark, I found that he lived in Greenland. I never understood the difference - however, this model solves that problem by moving the North Pole to Greenland! Viola!!!
| This work by Mark R Payne is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License. For details, see http://creativecommons.org/licenses/by-nc-sa/3.0/deed.en_US Basically, this means that you are free to “share” and “remix” for non-commerical purposes as you see fit, so long as you “attribute” me for my contribution. Derivatives can be distributed under the same or similar license. |
| This work comes with ABSOLUTELY NO WARRANTY or support. |
| This work should also be considered as BEER-WARE. For details, see http://en.wikipedia.org/wiki/Beerware |