Mapping temperature data from an .nc file in R

Question

I downloaded temperature data from [NARR] (https://www.esrl.noaa.gov/psd/data/gridded/data.narr.monolevel.html) specifically "Air temperature at 2m" -monthly mean

I opened the file using the package "ncdf4". The data has 4 dimensions- time, x, y, nbnds. y corresponds to lat and x corresponds to lon. There is a variable (not dimension) called air which I do not know how to use, although this is the temperature information.

My end goal is to map the temperature data on a map of North America, using averaged temperature data for each month for each year (12 maps, one for each month).

I am having trouble identifying how to use the data as all of the dimensions are just really long lists of numbers that don't seem to have meaning (eg. the x coordinates look like this: 6232896 6265359 6297822 6330285 6362748 6395211 6427674 6460137 6492600 6525063 6557526 6589989, and so do the y values and time).

Here is the code I am using to view the dimensions:

temp2m <- nc_open("air.2m.mon.mean.nc")
time <- temp2m$dim$time$vals
lat <- temp2m$dim$x$vals
lon <- temp2m$dim$y$vals
nbnds <- temp2m$dim$nbnds$vals

If someone could help me view the data as well as map temperature data onto North America that would be great.

Thank you!

Answer 1

You can use the raster package to read these into a stack:

> library(raster)
> air = stack("./air.2m.mon.mean.nc")

(Note, you may need a raster package compiled with netcdf drivers...)

You can then plot them by slice or by time-name:

> plot(air[[23]])
> plot(air[["X1979.10.01.01.01.15"]])

The stack prints like this:

> air
class       : RasterStack 
dimensions  : 277, 349, 96673, 450  (nrow, ncol, ncell, nlayers)
resolution  : 32462.99, 32463  (x, y)
extent      : -16231.49, 11313351, -16231.5, 8976020  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +x_0=5632642.22547 +y_0=4612545.65137 +lat_0=50 +lon_0=-107 +lat_1=50 +lat_2=50 +ellps=WGS84 
names       : X1979.01.01.00.01.15, X1979.02.01.00.01.15, X1979.03.01.00.01.15, X1979.04.01.01.01.15, X1979.05.01.01.01.15, X1979.06.01.01.01.15, X1979.07.01.01.01.15, X1979.08.01.01.01.15, X1979.09.01.01.01.15, X1979.10.01.01.01.15, X1979.11.01.00.01.15, X1979.12.01.00.01.15, X1980.01.01.00.01.15, X1980.02.01.00.01.15, X1980.03.01.00.01.15, ... 

and those coordinates are not really lat-long, but are in a transformed coordinate system described by that "coord. ref." string. If you want to put it on a lat-long map you need to warp it:

> air_ll = projectRaster(air[[1]],crs="+init=epsg:4326")
> plot(air_ll)

It might be better for you to transform any other data to this system, and keep the grid unprojected. Just look up how to deal with spatial data in R for more info on projections and transformations.