METEO 810
Weather and Climate Data Sets

Time Series NARR data

Prioritize...

You will learn how to extract different "slices" through a NetCDF dataset by retrieving variables in various dimensions. 

Read...

Now that you've performed a basic plot from a NetCDF file, let's try something more interesting. Let's plot the daily mean temperature for a point closest to a chosen latitude and longitude. For this task, we will want to go back to the PSD mono-level website and get a daily means file for temperature. I'm sticking with 1988, however, notice that the file names for the means are the same as the 3-hourly data (so I renamed my new file: "air.sfc.1988mean.nc"). Here's the code:

# if for some reason we need to reinstall the
# package, then do so.

if (!require("ncdf4")) {
  install.packages("ncdf4") }

# include the library
library(ncdf4)

# open the CDF file
nc <- nc_open("air.sfc.1988mean.nc")

# instead of a specific date, we need a location
location<-c(40.85, -77.85)

# How big is the dataset?
data_dims<-nc$var$air$varsize

#Get the properly formatted date strings
ncdates = nc$dim$time$vals  #hours since 1800-1-1
origin = as.POSIXct('1800-01-01 00:00:00', tz = 'UTC')
ncdates = origin + ncdates*(60*60)

# obtain the lat/lon grids
lats=ncvar_get(nc, varid = 'lat', 
               start = c(1,1),
               count = c(data_dims[1],data_dims[2]))

lons=ncvar_get(nc, varid = 'lon', 
               start = c(1,1),
               count = c(data_dims[1],data_dims[2]))
lons<-ifelse(lons>0,lons-360,lons)

# Here's where we find the grid index of closest point
distance<-(lats-location[1])**2+(lons-location[2])**2
closest_point<-which(distance==min(distance), arr.ind = TRUE)

# create the temperature array... Notice the different dim parameter
temp_out = array(data = NA, dim = c(1,1,data_dims[3]))

# get the temperature array and convert it
temp_out[,,] = ncvar_get(nc, varid = 'air', 
                         start = c(closest_point[1],closest_point[2],1),
                         count = c(1,1,data_dims[3]))
temp_out[1,1,]<-(temp_out[1,1,]-273.15)*9/5+32

# Create a standard line plot but with no x-axis lables
plot(ncdates,temp_out[1,1,],type="n",
     xlab = "Date", ylab = "Temperature (F)",
     xaxt="n",
     xlim=c(as.POSIXct('1988-01-01 00:00:00', tz = 'UTC'),
            as.POSIXct('1988-12-31 00:00:00', tz = 'UTC')),
     main = paste0("NARR Mean Daily Surface Temperature\n University Park, PA\n",
                   head(ncdates,1)," to ",tail(ncdates,1))
)
lines(ncdates,temp_out[1,1,],lwd=2,lty=1,col="red")

# create some custom labels
labs<-format(ncdates, "%b-%d")
spacing<-seq(1,length(ncdates),30)
axis(side = 1, at = ncdates[spacing], labels = labs[spacing], 
     tcl = -0.7, cex.axis = 1, las = 1)

# close the file
nc_close(nc)

Here is the graph of mean daily surface temperature produced by the code above. You can narrow the view by changing the dates in the xlim plot parameter. How might these daily means compare with the observed daily means? Hmm.... that's a good question to try and answer. But before we get to that, let's look at some vertical data from the NARR.

Read on.