We have one of these files on hand:
## dp is the root to our local data repository
f <- file.path(dp, "data", "ftp.cdc.noaa.gov/Datasets/ncep.reanalysis2.derived/pressure/air.mon.mean.nc")
library(raster)
b <- brick(f, level = 1)
b
On my system that uses the "ncdf4" package, but it could also use "ncdf".
The global Pacific-view extent looks close, but I would reset it after checking carefully in the file:
extent : -1.25, 358.75, -91.25, 91.25 (xmin, xmax, ymin, ymax)
library(ncdf4)
nc <- nc_open(f)
## print(nc) to see details
lon <- ncvar_get(nc, "lon")
range(lon)
##[1] 0.0 357.5
So, the longitudes that start at 0 and end somewhere near 0 need to be fixed, and so do the latitudes. I think it's pretty uncontroversial to assume they mean to fill the whole world (raster makes the reasonable assumption that each coordinate belongs to the centre, but clearly that's not what these are. It's typical for these incomplete specifications to exist, you are expected to fix it yourself. NetCDF has no standard concept of an affine transformation afaics so it's stuck with these redundant and often broken 1D coordinate arrays for regular grids.
extent(b) <- extent(0, 360, -90, 90)
The function rotate() will put the object into Atlantic view.
You can read other levels ("air" is a 4D variable with 432 time steps atm and 17 height levels) by setting that in the brick() argument.
I tend to always work this way:
use raster to see what happens
read the metadata with some netcdf package (ncdf, ncdf4, RNetCDF, or possibly GDAL or ncdump on the command line, occasionally Manifold GIS)
use those investigations to hone the use of raster.
You can't really use raster for non-regular coordinates bound to arrays, unless you really know what you're doing (but this air temp data is fine for XYT). See the available levels with
ncvar_get(nc, "level")
#[1] 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
10 is level = 1 for raster, etc.
Using raster you can get to individual months like this:
b[[10]]
subset(b, 10:13)
but now we get into "go learn the raster package" territory.
I use the "CountriesLow" data set in "rworldmap" to get mean country temperatures:
asub <- which(as.Date(getZ(b)) >= as.Date("2001-01-01"))
library(rworldmap)
data(countriesLow)
## extract the first year's data
imonth <- 1
d <- data.frame(Name = countriesLow$SOVEREIGNT, meanairtemp = extract(rotate(b[[asub[imonth]]]), countriesLow, fun = mean, na.rm = TRUE),
date = as.Date(getZ(b)[asub[imonth]]), stringsAsFactors = FALSE)
I've assumed you would build a long-form table and just append each year to this, but you might want a column per month? One of the values comes out as NA but I think that can be controlled by small/sampling options with extract (your polygon layer may differ).
Don't forget to rotate() to get the Pacific view to match the Atlantic polygons, and that the data is in Kelvin.
Use rgdal::readOGR to get your shapefile or whatever in. extract() will automatically reproject from different projections, but that won't always work or apply (certainly not for the Pacific/Atlantic conventions).
