Recipe: IDL program to read data from downloaded file
Back to G5NR Data Access Guide.
Problem
We want to read a downloaded data file using IDL.
Solution
For the purpose of this example, we assume that we have already downloaded the file c1440_NR.inst01hr_3d_T_Cv.20060918_0900z.nc4
from the ftp/http server. For more information about file naming conventions, and how to download a file from the ftp server, please follow the links in the #See Also section.
Code
The code below reads the global temperature data from file, computes the maximum and minimum temperatures and plots (using the matplotlib
package) the aire temperature at the surface (level=71). It has been tested with IDL v 8.3.
; file
file = 'c1440_NR.inst01hr_3d_T_Cv.20060918_0900z.nc4'
; read data from file
ncid = ncdf_open(file)
ncdf_varget, ncid, 'T', T
ncdf_close, ncid
; max/min of T
print, 'max(T):', max(T)
print, 'min(T):', min(T)
; plot T at the surface
lons = findgen(720)*0.5-180.
lats = findgen(361)*0.5-90.
basemap = map('Cylindrical Equal Area', limit=[-90., -180., 90., 180.])
cntrplot = contour( $
T[*,*,71], $ ;; level=71 => surface
lons, $
lats, $
grid_units=2, $ ;; degrees
n_levels=20, $
rgb_table=34, $
;; /fill, $
/overplot, $
title='surface air temperature' $
)
cntnt = mapcontinents(/countries)
END
Output
Running this IDL script
IDL> .run g5nr_reader.pro
we get the text output
max(T): 315.651 min(T): 180.367
and the plot
Discussions
Modifications to read a subset of the data
The above IDL script can be easily modified to read a subset of the data instead. The modified code, to read temperature data inside the box bounded by latitudes 25oN, 50oN and longitudes -130oW, -65oW, is
; file
file = 'c1440_NR.inst01hr_3d_T_Cv.20060918_0900z.nc4'
; read data from file
ncid = ncdf_open(file)
; bounding box
; lons = -130:0.5:-65
; lats = 25:0.5:50
imin = round((-130. + 180.)/0.5) - 1
imax = round(( -65. + 180.)/0.5) - 1
jmin = round(( 25 + 90.)/0.5) - 1
jmax = round(( 50 + 90.)/0.5) - 1
; corresponding array sizes
im = imax-imin+1
jm = jmax-jmin+1
lm = 72 ; read all 72 levels
; start and count vectors
offset = [imin, jmin, 0, 0]
count = [im, jm, lm, 1]
ncdf_varget, ncid, 'T', T, count=count, offset=offset
ncdf_close, ncid
; max/min of T
print, 'max(T):', max(T)
print, 'min(T):', min(T)
; plot T at the surface
lons = findgen(im)*0.5-130.
lats = findgen(jm)*0.5+25.
basemap = map('Cylindrical Equal Area', limit=[25., -130., 50., -65.])
cntrplot = contour( $
T[*,*,71], $ ;; level=71 => surface
lons, $
lats, $
grid_units=2, $ ;; degrees
n_levels=20, $
rgb_table=34, $
;; /fill, $
/overplot, $
title='surface air temperature' $
)
m1 = mapcontinents(/countries)
END
Running the modified script, produces the text output
max(T): 305.605 min(T): 191.696
and the plot