How to read and plot netCDF files with R

Reading and plotting netCDF files with R packages ncdf and maps. With longitude range from 0 to 360 and  from -180 to 180. The netCDF data used in this tutorial is the 4-times Daily Air Temperature of the NCEP/NCAR Reanalysis-1 Pressure data provided by NOAA (National Oceanic and Atmospheric Administration).

Download netCDF

Selecte variable and year to download the netCDF file in the working directory.
URL: ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/pressure/air.2015.nc

# Select: variable and year
var='air'
year=2015

# Download netCDF file
url=paste('ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/pressure/',var,'.',year,'.nc', sep='')
dfile=paste(var,'.',year,'.nc', sep='')
download.file(url, destfile=dfile)

Read netCDF and get selection

To open and read netCDF files the R package ncdf is used.

# Library
library(ncdf)

# Open netCDF file
nc=open.ncdf(dfile)

The list nc has 4 dimensions: Longitude, Latitude, Level, and Time. With the following commands you can look at each dimension:

# 4 dimensions: lon,lat,level,time
nc
nc$dim$lon$vals
nc$dim$lat$vals
nc$dim$level$vals
nc$dim$time$vals
nc$dim$time$units

The times are encoded as the number of hours since 1800-01-01 00:00:0.0 (NOAA: How to Read PSD’s netCDF Files)

In the next step the level 1000 and a short time period is selected with the function get.var.ncdf(). This function gets the data from the starting point and reads along the count of values. To use this function the indices of the selected level and timestamps must be identified. The value “-1” indicates that all entries along that dimension should be read. This is used for longitude and latitude. After getting the selected array the nc file can be closed.

# Level
lev=1000
z1=which(nc$dim$level$val==lev)
dz=1

# Time
time.s=as.POSIXct('2015-01-01 12:00',tz='UTC')
time.e=as.POSIXct('2015-01-02 12:00',tz='UTC')
tseq=seq(time.s, time.e, by='6 hours')
times=as.POSIXct(nc$dim$time$vals*60*60, origin='1800-01-01 00:00:0.0', tz='UTC')
t1=which(times==time.s)
t2=which(times==time.e)
dt=t2-t1+1

# Get array from netCDF: x,y,z,t
nc.a=get.var.ncdf(nc, var, start=c(1,1,z1,t1), count=c(-1,-1,dz,dt))

# Close netCDF file
close.ncdf(nc)

Plot netCDF

To display maps of the word in the plot the R package maps is used.

library(maps)

There are two different ways to plot the air temperature data:

Plot longitude from 0 to 360

In order to plot the data, the latitude has to be set from [90 -90] to [-90 90] and the array has to be arrange in the same way. As the longitudes goes from [0 360] the map ‘word2’ should be used in the plot. The unit of the temperature is Kelvin. To change the unit to °C subtract 273.5 from the array. The vector int contains 30 equal intervals of the temperature values.

## Set lat from [90 -90] to [-90 90]
lat=nc$dim$lat$vals
py=c(length(lat):1)
lat=lat[ py ]
nc.a=nc.a[,py,]

## Set lon [0 357.5]
lon=nc$dim$lon$vals

## Plot settings
rgb.palette=colorRampPalette(c('darkblue','palegreen1','yellow','red2'),interpolate='spline')
nc.a=nc.a-273.15
int=seq(min(nc.a),max(nc.a),length.out=30)

## Plot first selected timestamp with map 'world2'
i=1
filled.contour(lon, lat, nc.a[,,i], color.palette=rgb.palette, levels=int,
    plot.title=title(main=tseq[i], xlab='Longitude [°]', ylab='Latitude [°]'),
    plot.axes={axis(1); axis(2);map('world2', add=TRUE);grid()})

lon_0-360Plot longitude from -180 to 180

You can also plot the data with the longitude from -180 to 180. Therefore you have to rearrange the array. In this case the map ‘word’ should be used in the plot:

# Set lon [-180 177.5]
lon=nc$dim$lon$vals
k=which(lon==180)
px=c(k:length(lon),1:(k-1))
lon2=lon[px]
nc.a2=nc.a[px,,]
lon2[1:(which(lon2==0.0))-1]=lon2[1:(which(lon2==0.0)-1)]-360

# Plot with map 'world'
filled.contour(lon2, lat, nc.a2[,,i], color.palette=rgb.palette, levels=int,
    plot.title=title(main=tseq[i], xlab='Longitude [°]', ylab='Latitude [°]'),
    plot.axes={axis(1); axis(2);map('world', add=TRUE);grid()})

 

lon_-180-180