Retrieving Vertical Profiles

Read...

We've looked at NARR data in both space and time, now let's look at retrieving vertical data. We'll need to retrieve a levels file from the PSD Pressure Levels page. I grabbed the daily mean air temperature file for June 1988 and renamed the file "air.198806levels.nc". Notice that for the mono-level data, I was able to retrieve an entire year at a time, but for these data, I am only able to get a month's worth. To see why, run the following three lines of code from the command line (you can skip the library command if ncdf4 is already loaded).

> library(ncdf4)
> nc <- nc_open("air.198806levels.nc")
> nc

This, as you know, will dump the header for the file, giving you information about the dimensions of the file data. Look closely at the temperature variable "air". Can you see what's different? Instead of just 3 dimensions (lat, lon, time), this variable now has 4 dimensions (lat, lon, level, time). Furthermore, if you query the size of air with: nc$var$air$varsize, you will see that there are 29x30 grids in each file. That's why each level file covers only one month (each month's file is 2.5X larger than an entire year's worth of surface data).

Now, open up a new R-script, and let's look at how we might explore this data. My goal is to extract a daily mean temperature profile from the grid point nearest my lat/lon location and on a day of my choosing (in June). Could you modify the time-series code we used on the last page? Give it a try... or start with the code below... (you've seen much of this before).

# if for some reason we need to reinstall the
# package, then do so.
if (!require("ncdf4")) {
  install.packages("ncdf4") }

library(ncdf4)

filename<-"air.198806levels.nc"

nc <- nc_open(filename)

# choose a location/time
location<-c(40.85, -77.85)
target_datetime<-"1988-06-04"

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

#Get the properly formatted date strings/target time index
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)
time_index<-which(ncdates==as.POSIXct(target_datetime, tz='UTC'))

#Get the pressure levels
nclevels <- nc$dim$level$vals  #in mb

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)

# Find the closest point by minimizing the distance
distance<-(lats-location[1])**2+(lons-location[2])**2
closest_point<-which(distance==min(distance), arr.ind = TRUE)

In this code, I combine techniques of the spatial grid plotting script and the time-series script. Basically, I need to find the correct time and the lat/lon indices. Then, I can get the data and plot it. Notice that temp_out now has four dimensions, not three. Also, notice that I retrieve all of the times and levels for my given location. This is only a matrix of 870 values, so I'm not concerned with space. This way, I can plot not only a single profile, but a series of profiles or even a time-series at a particular level if I want.

temp_out = array(data = NA, dim = c(1,1,data_dims[3],data_dims[4]))

temp_out[,,,] = ncvar_get(nc, varid = 'air', 
                         start = c(closest_point[1],closest_point[2],1,1),
                         count = c(1,1,data_dims[3],data_dims[4]))
# change to degrees C
temp_out[1,1,,]<-(temp_out[1,1,,]-273.15)

# Do a line plot of the profile
plot(temp_out[1,1,,time_index],nclevels,type="n",
     xlab = "Temperature (C)", ylab = "Pressure",
     ylim=c(1000,200), log="y",
     main = paste0("NARR Mean Daily Temperature Profile\n University Park, PA\n",
                   ncdates[time_index])
)
lines(temp_out[1,1,,time_index],nclevels,lwd=2,lty=1,col="red")

Check out the profile of temperature plot. Notice that I changed a few lines in the plot routine. First, I reversed the order of the y-axis with ylim=c(1000,200). This puts the 1000mb tick on the bottom and lets pressure decrease upwards. Second, I made the y-axis log(P) with the parameter: log="y". If you wanted to think about it, you might play with how to plot temperature on a skewed-T axis. You might also want to check out the package RadioSonde for plotting proper skew-T/Log P diagrams (some data reformatting will be required). If someone wants to explore this package, I'll add your write-up as an Explore Further (and give you the by-line!).

How about one more way of looking at the data (this one's pretty interesting). Replace the plotting commands for the temperature profile with the following code:

# *** ADD THIS TO THE TOP!
library(fields)    # needed for image.plot

# *** REPLACE PLOT(...) AND LINES(...) WITH:
# Create a 2D matrix of values and reorder the columns of the matrix
z<-t(temp_out[1,1,,])
z<-z[,ncol(z):1]

# Create an image plot 
colormap<-rev(rainbow(100))
image.plot(ncdates, rev(nclevels), z, col = colormap,
           ylim=c(1000,500),log="y", zlim=c(-30,30),
           xlab = "Date", ylab = "Pressure (mb)",
           main = paste0("NARR Mean Daily Temperature Profiles\n University Park, PA\n June 1988")
)

This plots all of the profiles as a time-series and gives you a contour plot of sorts. In order to fit the restrictions of the image.plot(...) command, I had to pull out the temperature data into a 2-D array, transpose the rows and columns, and finally reverse the order of the levels (must be increasing values). Then I was able to plot the image (with the lowest 500mbs plotted, log y-axis, and limiting the colors to a -30C to 30C temperature range). Check out the resulting image... what do you observe?