@time_lib
;calculate the field size based on the size of the dimensions:
function calc_field_size, dim
field_size=1L
ndim=n_elements(dim)
for i=0, ndim-1 do begin
if dim[i] ne 0 then field_size=field_size*dim[i]
endfor
return, field_size
end
;writes the file headers in a customized binary format:
;-dimensions are the dimensions of the fields
;-dim_names are the names of the dimensions
;-vars are the names of the fields contained in the file
;-err returns an error code
pro write_headers, unit, dimensions, dim_names, vars, err=err
on_ioerror, fail
magic=-49L
ndim=4L
; dimensions2=lonarr(ndim)
; dimensions2[0]=dimensions
dimensions2=long(dimensions)
dim_names2=strarr(ndim)
dim_names2[0]=dim_names
dim_len=strlen(dim_names2)
var_len=strlen(vars)
nvars=n_elements(vars)
offset=long(24+ndim*8+total(dim_len)+nvars*9+total(var_len)+total(dimensions[0:2])*4)
fieldsize=4*dimensions2[0]*dimensions2[1]*dimensions2[2]
; openw, unit, filename, /get_lun
dim_type=bytarr(ndim)+4b ;reserve some space to hold the types of each of the dimensions
var_type=bytarr(nvars)+4b ; " " " " " " vars.
point_lun, unit, 0
writeu, unit, magic, ndim, nvars, fieldsize, offset, dimensions2, dim_type, dim_len, dim_names2
writeu, unit, var_type, var_len, vars, lonarr(nvars)
;(the last parameter in the write statement is reserved for variable descriptors... )
point_lun, -unit, cur
; print, "Offset; current:", cur, " predicted:", offset-total(dimensions[0:2])*4
; stop
; free_lun, unit
err=0
return
fail:
err=-1
end
;reads the file headers in a customized binary format:
;-dimensions are the dimensions of the fields
;-dim_names are the names of the dimensions
;-vars are the names of the fields contained in the file
;-err returns an error code
;-offset returns the starting offset of the first record
;-fieldsize returns the size in bytes of each record
pro read_headers, unit, dims, dim_names, var_names, offset=offset, $
err=err, fieldsize=fieldsize
; openr, unit, filename, /get_lun ;open the file
on_ioerror, fail
point_lun, unit, 4 ;skip the "magic" number
ndim=0L ;number of dimensions (should be 4)
nvar=0L ;number of variables
fieldsize=0L ;size in bytes of one field
offset=0L ;offset after reading headers and dimensions
readu, unit, ndim, nvar, fieldsize, offset
dims=lonarr(ndim) ;read the dimension sizes and name lengths
dim_len=lonarr(ndim)
dim_type=bytarr(ndim)
readu, unit, dims, dim_type, dim_len
dim_names=string(strarr(ndim), format='(a'+strjoin(string(dim_len), ',/,a')+')')
readu, unit, dim_names
;read the variable names and descriptors:
var_type=bytarr(nvar)
var_len=lonarr(nvar)
desc_len=lonarr(nvar)
readu, unit, var_type, var_len
var_names=string(strarr(nvar), format='(a'+strjoin(string(var_len), ',/,a')+')')
readu, unit, var_names
readu, unit, desc_len
desc=string(strarr(nvar), format='(a'+strjoin(desc_len, ',/,a')+')')
readu, unit, desc
; stop
; free_lun, unit
err=0
return
fail:
err=-1
end
;finds the location of a record starting from the current
;position using a binary search:
;-unit is the file unit
;-offset is the starting position in bytes of the first record
;-rec_size is the size of each record
;-nrec is the number of records
;-time is the date to search on
;-found is 1 if the record was found, 0 if not
function find_field, unit, offset, rec_size, nrec, time, found
on_ioerror, fail
first=0
last=nrec-1
first_time={time_str} & last_time={time_str} & mid_time={time_str}
;read the dates of the first and last records:
point_lun, unit, offset+first*rec_size
readu, unit, first_time
point_lun, unit, offset+last*rec_size
readu, unit, last_time
;check to see if the date is less than or equal to that of the first record:
comp=time_compare(time, first_time)
if comp eq 0 then begin
found=1
return, first
endif else if comp lt 0 then begin
found=0
return, first
endif
;check if the date is greater than or equal to that of the last record:
comp=time_compare(time, last_time)
if comp eq 0 then begin
found=1
return, last
endif else if comp gt 0 then begin
found=0
return, last+1
endif
while last-first gt 1 do begin
mid=(first+last)/2
point_lun, unit, offset+mid*rec_size
readu, unit, mid_time
comp=time_compare(time, mid_time)
if comp eq 0 then begin
found=1
return, mid
endif else if comp lt 0 then begin
last=mid
endif else begin
first=mid
endelse
endwhile
found=0
return, last
fail:
return, -1L
end
;inserts a field into an existing file:
;offset is where the data actually begins
;nvar is the number of types of data
;nrec is the number of dates in the file
;field is the data itself
;time is the date
;returns 0 if the record existed, 1 if it didn't,
;-1 if an I/O error occurred
function insert_field, unit, offset, nvar, var_ind, nrec, field, time, $
overwrite=overwrite
on_ioerror, fail
field_size=n_elements(field)*4 ;size in bytes of one global field for one date
rec_size=nvar*field_size+14 ;size of all global fields for one date
;first search for the field:
ind=find_field(unit, offset, rec_size, nrec, time, found)
if ind lt 0 then return, ind
;if the field was found only overwrite if specified:
if found then begin
if keyword_set(overwrite) then begin
point_lun, unit, offset+ind*rec_size+var_ind*field_size+14
writeu, unit, float(field)
endif else begin
message, "Field exists, ("+time_string(time)+") was not overwritten", /continue
endelse
new=0
;if the field was not found, then insert it:
endif else begin
;shift all subsequent records forward:
record={time:time, data:make_array(/float, dimension=[n_elements(field), nvar])}
data=assoc(unit, record, offset, /packed)
for i=nrec-1, ind, -1 do begin
data[i+1]=data[i]
endfor
;write the new record:
record.data[*,var_ind]=field
data[ind]=record
new=1
endelse
; stop
return, new
fail:
return, -1
end
pro write_field, filename, name, field, time, lon, lat, lev, overwrite=overwrite, $
clobber=clobber, nocheck=nocheck, err=err, continue_if_err=continue
;+***********************************************************************************
; WRITE_FIELD
;************************************************************************************
;
; usage: write_field, filename, fieldname, data, time, lon, lat, lev
;
; purpose: Experimental routine to store one or more 3-d atmospheric fields in a
; customized binary format which is optimized for this type of data.
;
; Features of the format:
; -stores 3-d fields at multiple time grids
; -unlimited time dimension
; -fields can be inserted at arbitrary time grids
; -multiple fields can be stored in one file
; -data can be retrieved all at once or by searching on a single time grid
; -fields can be read or written to using a single subroutine call
; -very fast
; -portable across different platforms
;
; parameters: filename The name of the file in which to store the field.
;
; fieldname The name of the field(s) to store. If you are saving to an
; old file, it must contain a field of this name. If you are
; saving to a new file, this should contain the names of all
; the fields you would like stored in the file. Data is written
; to the first field in the list.
;
; data The data to store in the file. If saving to a new file,
; should have four dimensions: the first corresponds to
; longitude, the second to latitude, the third to vertical
; level and the fourth to time. If saving to an existing file
; should have only three dimensions. (Time is omitted.)
;
; time The time grid. If saving to an existing file, can store
; only one time grid at a time, but data can be added at any
; location. Time data must be in the form of a structure, a
; string or a longword integer. (see "time_lib.pro" and
; "convert_date")
;
; lon The longitude grid. If saving to an existing file, should
; agree with the data already in there.
;
; lat The latidutde grid. If saving to an existing file, should
; agree with the data already in there.
;
; lev The vertical levels.. If saving to an existing file, should
; agree with the data already in there.
;
; keywords: /overwrite Tells the routine to overwrite existing fields.
;
; /clobber Tells the routine to always create a new file and clobber
; any existing file of the same name.
;
; /nocheck If writing to an existing file, normally the longitude, latitude
; and vertical grids are checked against those already saved in
; the file. This tells the routine to skip that, thereby making
; those arguments optional.
;
; err An error code. Zero indicates success.
;
; /continue_if_err If an error occurs, tells the routine to return and
; continue in the calling routine. Default behaviour is
; to stop.
;
; dependencies: All of the routines in this file except "read_field", below.
; Also: time_lib.pro, convert_date
;
; history: First formally documented 2002-6-24.
; Written by Peter Mills (peter.mills@nrl.navy.mil)
; 2003-10-18 PM: detail changes to improve efficiency and logic
;
;-*****************************************************************************************
on_ioerror, fail
dims=[n_elements(lon), n_elements(lat), n_elements(lev), n_elements(time)]
time1=convert_date(time)
if size(field, /type) ne 4 then field=float(field)
;simplest case: file does not exist, write only the one field
if file_test(filename, /read, /write) eq 0 or keyword_set(clobber) then begin
;open the file:
openw, unit, filename, /get_lun, err=err
if err ne 0 then begin
message, "Could not open file '"+filename+"'", /ioerror
endif
;write out all the headers:
nvar=n_elements(name)
write_headers, unit, dims, ["lon", "lat", "lev", "time"], name, err=err
if err ne 0 then begin
message, "An I/O error occurred while writing the headers", /ioerror
endif
if calc_field_size(dims) ne n_elements(field) then begin
err=3
message, "Size of field does not match expected", /ioerror
endif
;write out all the dimensions:
if dims[0] ne 0 then writeu, unit, float(lon)
if dims[1] ne 0 then writeu, unit, float(lat)
if dims[2] ne 0 then writeu, unit, float(lev)
;make sure the field has four dimensions:
field=reform(field, dims[0], dims[1], dims[2], dims[3], /overwrite)
;if the user has specified only one variable, write is more simple:
if nvar eq 1 then begin
for i=0, dims[3]-1 do begin
writeu, unit, time1[i], field[*, *, *, i]
endfor
endif else begin
;if the user has specified multiple variables, must skip those not included:
skip=dims[0]*dims[1]*dims[2]*4*(nvar-1)
for i=0, dims[3]-1 do begin
writeu, unit, time1[i], field[*, *, *, i]
point_lun, -unit, loc ;skip locations for missing variables
point_lun, unit, loc+skip
endfor
endelse
free_lun, unit
endif else begin
;file exists, must insert field at the correct date locations, one by one:
openu, unit, filename, /get_lun, err=err
if err ne 0 then begin
message, "Could not open file '"+filename+"'", /ioerror
endif
; on_ioerror, fail
magic=0L
readu, unit, magic
if magic ne -49 then begin
if swap_endian(magic) eq -49 then begin
free_lun, unit
openu, unit, filename, /get_lun, /swap_endian
endif else begin
err=1
message, "Not a recognized file format", /ioerror
endelse
endif
;read the headers and compare them with the parameters passed to the subroutine:
read_headers, unit, read_dim, dim_names, var_names, offset=offset, fieldsize=field_size, err=err
if err ne 0 then begin
message, "An I/O error occurred while reading the headers", /ioerror
endif
ndim=n_elements(read_dim) ;number of dimensions in file
nvar=n_elements(var_names) ;number of variables in file
;check the number of dimensions (only 4 is allowed):
if ndim ne 4 then begin
message, "File has wrong number of dimensions", /ioerror
endif
;check the stored field size against that passed to the function:
; field_size=calc_field_size(read_dim[0:2])
if dims[3]*field_size ne 4*n_elements(field) then begin
message, "Field size does not match file headers", /ioerror
endif
;see if variable is actually stored in the file:
var_ind=(where(var_names eq name))[0]
if var_ind eq -1 then begin
message, "Field '"+name+"' is not stored in file", /ioerror
endif
;check the saved dimensions against those passed to the function:
if not keyword_set(no_check) then begin
if read_dim[0] ne 0 then begin
read_lon=fltarr(read_dim[0])
readu, unit, read_lon
endif
if read_dim[1] ne 0 then begin
read_lat=fltarr(read_dim[1])
readu, unit, read_lat
endif
if read_dim[2] ne 0 then begin
read_lev=fltarr(read_dim[2])
readu, unit, read_lev
endif
if dims[0] ne 0 then begin
if not array_equal(lon, read_lon) then begin
message, "Stored and supplied longitude data do not agree", /ioerror
endif
endif
if dims[1] ne 0 then begin
if not array_equal(lat, read_lat) then begin
message, "Stored and supplied latitude data do not agree", /ioerror
endif
endif
if dims[2] ne 0 then begin
if not array_equal(lev, read_lev) then begin
message, "Stored and supplied vertical level data do not agree", /ioerror
endif
endif
endif
;make sure the field has four dimensions:
field=reform(field, read_dim[0], read_dim[1], read_dim[2], dims[3], /overwrite)
;insert the field at the correct locations in the file:
nt=read_dim[3]
for i=0, dims[3]-1 do begin
added_new=insert_field(unit, offset, nvar, var_ind, nt, $
field[*, *, *, i], time1[i], overwrite=overwrite)
if added_new lt 0 then begin
message, "An I/O error occurred in insert_field", /ioerror
endif
nt=nt+added_new
endfor
;update the number of dates:
point_lun, unit, 32
writeu, unit, nt
free_lun, unit
endelse
;remove degenerate dimension that might have been introduced:
field=reform(field, /overwrite)
err=0
return
fail:
message, !error_state.msg, continue=continue
err=!error_state.code
if n_elements(unit) ne 0 then free_lun, unit
end
pro read_field, filename, name, data, time, lon, lat, lev, all=all, err=err, $
scan=scan, continue_if_err=continue
;+***********************************************************************************
; READ_FIELD
;************************************************************************************
;
; usage: read_field, filename, fieldname, data, time, lon, lat, lev
;
; purpose: Experimental routine to load one or more 3-d atmospheric fields from
; a binary file in a customized format which is optimized for this type of data.
;
; Features of the format:
; -stores 3-d fields at multiple time grids
; -unlimited time dimension
; -fields can be inserted at arbitrary time grids
; -multiple fields can be stored in one file
; -data can be retrieved all at once or by searching on a single time grid
; -fields can be read or written to using a single subroutine call
; -very fast
; -portable across different platforms
;
; parameters: filename The name of the file in which to store the field.
;
; fieldname The name of the field to retrieve. If "/scan" is set,
; returns all of the fields contained in the file.
;
; data The retrieved data. If the "/all" keyword is set, will
; have four dimensions: the first corresponds to longitude,
; the second to latitude, the third to vertical level and
; the fourth to time. Otherwise will have only three dimensions.
; (Time is omitted.)
;
; time If the "/all" or "/scan" keyword is set, will return all of
; the time grids in the file as an array of structures.
; Otherwise should contain the desired time grid at which to
; retrieve the data as a structure, string or longword integer.
; For information on the format of the time values
; consult the file "time_lib.pro" and the routine "convert_date."
;
; lon The retrieved longitude grid.
;
; lat The retrieved latidutde grid.
;
; lev The retrieved vertical levels.
;
; dependencies: All of the routines in this file except "read_field", below.
; keywords: /all Retrieves all of the time grids for the specified field.
;
; /scan Does not return any data, but simply scans the file to retrieve
; the names of all of the fields contained in the file, as well as
; the longitude, latitude, vertical and time grids.
;
; err An error code. Zero indicates success.
;
; /continue_if_err If an error occurs, tells the routine to return and
; continue in the calling routine. This is the default.
;
; dependencies: All of the routines in this file except "read_field", below.
; Also: time_lib.pro, convert_date
;
; history: First formally documented 2002-6-24.
; Written by Peter Mills (peter.mills@nrl.navy.mil)
; 2003-10-18 PM: changed documentation to match an earlier revision
;
;-*****************************************************************************************
if n_elements(continue) ne 1 then continue=1
on_ioerror, fail
openr, unit, filename, /get_lun
magic=0L
readu, unit, magic
if magic ne -49 then begin
if swap_endian(magic) eq -49 then begin
free_lun, unit
openr, unit, filename, /get_lun, /swap_endian
endif else begin
message, "Not a recognized file format", /ioerror
endelse
endif
read_headers, unit, dim, dim_names, var_names, offset=offset, fieldsize=field_size, err=err
if err ne 0 then begin
message, "An I/O error occurred while reading the headers", /ioerror
endif
; print, "Dimensions: ", dim
ndim=n_elements(dim)
if ndim ne 4 then begin
message, "File has wrong number of dimensions", /ioerror
endif
nvar=n_elements(var_names)
; time1=systime(1)
if dim[0] ne 0 then begin
lon=fltarr(dim[0])
readu, unit, lon
endif
if dim[1] ne 0 then begin
lat=fltarr(dim[1])
readu, unit, lat
endif
if dim[2] ne 0 then begin
lev=fltarr(dim[2])
readu, unit, lev
endif
; time2=systime(1)
; readu, unit, lon, lat, lev
; print, "Reading dimensions:", time2-time1
rec_size=nvar*field_size+14
if keyword_set(all) then begin
var_ind=(where(var_names eq name))[0]
if var_ind eq -1 then begin
message, "Field '"+name+"' is not stored in file", /ioerror
endif
; time3=systime(1)
;get all the dates:
slice=fltarr(dim[0], dim[1], dim[2])
data=make_array(/float, dimension=dim)
time_read={time_str}
time=replicate({time_str}, dim[3])
; time4=systime(1)
; print, "Initializing arrays:", time4-time3
; stop
for i=0, dim[3]-1 do begin
; time5=systime(1)
point_lun, unit, offset+rec_size*i
readu, unit, time_read
time[i]=time_read
point_lun, unit, offset+rec_size*i+var_ind*field_size+14
readu, unit, slice
; time6=systime(1)
data[*, *, *, i]=slice
; time7=systime(1)
; print, "Field,", i, time6-time5, time7-time6
endfor
endif else if keyword_set(scan) then begin
; time2=systime(1)
name=var_names
data=offset+rec_size*lindgen(dim[3])
time_read={time_str}
time=replicate({time_str}, dim[3])
for i=0, dim[3]-1 do begin
point_lun, unit, offset+rec_size*i
readu, unit, time_read
time[i]=time_read
endfor
; time3=systime(1)
; print, "Scan,", time3-time2
endif else begin
var_ind=(where(var_names eq name))[0]
if var_ind eq -1 then begin
message, "Field '"+name+"' is not stored in file", /ioerror
endif
time1=convert_date(time)
ind=find_field(unit, offset, rec_size, dim[3], time1, found)
if ind lt 0 then begin
message, "An I/O error occurred in find_field", /ioerror
endif
if found ne 1 then begin
message, "Date doesn't exist in file ("+time_string(time)+")", /ioerror
endif
data=fltarr(dim[0], dim[1], dim[2])
point_lun, unit, offset+rec_size*ind+var_ind*field_size+14
readu, unit, data
endelse
free_lun, unit
err=0
return
fail:
message, !error_state.msg, continue=continue
err=!error_state.code
if n_elements(unit) ne 0 then free_lun, unit
end