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.