source: trunk/step1_diff.pro @ 48

;+
;
; @file_comments
; compute delta between cond_sed*.pro Netcdf files with same dimension
;
; @param file1 {in}{required} {type=string}
; name of the first file to be read
;
; @param file2 {in}{required} {type=string}
; name of the second file to be read
;
; @param file3 {in}{required} {type=string}
; file where differences between variables are written
;
; @examples
; test (for header and delta)
; to compute difference between cond_sed in cond_sed_ORCA2.nc and
; itself and write it in ginette.nc :
; $ cd ${GEOMAG}
; $ idl
; IDL> file1=getenv('GEOMAG') + '/cond_sed_ORCA2.nc'
; IDL> file2=getenv('GEOMAG_OD') + '/cond_sed_ORCA2.nc'
; IDL> step1_diff, file1, file2, 'ginette.nc'
; values of delta must be 0 everywhere
; $ rm ginette.nc
; idem with BR file
; IDL> file1=getenv('GEOMAG') + '/Br_ORCA2.nc'
; IDL> file2=getenv('GEOMAG_OD') + '/Br_ORCA2.nc'
; IDL> step1_diff, file1, file2, 'ginette.nc'
;
; real life
; to compute difference between cond_sed in cond_sed_ORCA2.nc and
; /home/rech/eee/reee522/geomag_d/cond_sed_ORCA2.nc and write it in
; cond_sed_ORCA2_diff.nc :
; $ cd ${GEOMAG}
; $ idl
; IDL> file1=getenv('GEOMAG') + '/cond_sed_ORCA2.nc'
; IDL> file2=getenv('GEOMAG_OD') + '/cond_sed_ORCA2.nc'
; IDL> step1_diff, file1, file2, 'cond_sed_ORCA2_diff.nc'
; IDL> file1=getenv('GEOMAG') + '/Br_ORCA2.nc'
; IDL> file2=getenv('GEOMAG_OD') + '/Br_ORCA2.nc'
; IDL> step1_diff, file1, file2, 'Br_ORCA2_diff.nc'
;
; to see the difference, for example
; IDL> xxx,'cond_sed_ORCA2_diff.nc'
; IDL> xxx,'Br_ORCA2_diff.nc'
; select xy on plt wigdet
;
; @history
; reee522 2006-12-18T10:38:49Z rhodes (IRIX64)
; thanks to example in SAXO
; SRC/ReadWrite/idl-NetCDF/ncdf_quickread/README_ncdf_quickread.txt
;
; @todo
; check difference between headers
;
; @version
; $Id$
;
PRO step1_diff, file1, file2, file3
;
ncverbose=1
;
; check if file1 exists
fullfile1 = isafile(file1, NEW=0,/MUST_EXIST)
IF fullfile1[0] EQ '' THEN BEGIN
   msg = 'eee : the file ' + fullfile1 + ' was not found.'
   ras = report(msg)
   RETURN
ENDIF
;
; protection
IF (FILE_TEST(fullfile1[0], /READ) EQ 0) THEN BEGIN
   msg = 'eee : the file ' + fullfile1[0] + ' is not readable.'
   ras = report(msg)
   RETURN
ENDIF
;
; check if file2 exists
fullfile2 = isafile(file2, NEW=0,/MUST_EXIST)
IF fullfile2[0] EQ '' THEN BEGIN
   msg = 'eee : the file ' + fullfile2 + ' was not found.'
   ras = report(msg)
   RETURN
ENDIF
;
; protection
IF (FILE_TEST(fullfile2[0], /READ) EQ 0) THEN BEGIN
   msg = 'eee : the file ' + fullfile2[0] + ' is not readable.'
   ras = report(msg)
   RETURN
ENDIF
;
; read the first file
ncdf_read,  fullfile1[0], file1info, file1dinfo, file1vinfo, $
            file1gatts, file1vatts, file1data
;
; read the second file
ncdf_read,  fullfile2[0], file2info, file2dinfo, file2vinfo, $
            file2gatts, file2vatts, file2data
;
; compute delta
delta = *file1data[3] - *file2data[3]
minarr = min(delta, max = maxarr)
;
; check for number of non zero in delta
checknonzero=where(delta NE 0.,count)
IF count EQ 0 THEN BEGIN
   msg = 'iii : delta is zero everywhere'
   ras = report(msg)
ENDIF ELSE BEGIN
   msg = 'iii : delta is not zero ' + STRING(count)  + ' times'
   ras = report(msg)
ENDELSE
;
; write output
; in order to avoid unexpected overwritten
IF (FILE_TEST(file3) EQ 1) THEN BEGIN
   msg = 'eee : the file ' + file3  + ' already exists.'
   ras = report(msg)
   RETURN
ENDIF
;
;---------------------------
; Creation of the NetCdf output file
;---------------------------
;
  netcdf_id = NCDF_CREATE(file3, /clobber)
  NCDF_CONTROL, netcdf_id, /NOFILL
;
; dimension
  dimidx = NCDF_DIMDEF(netcdf_id, file1dinfo[0].name ,  file1dinfo[0].size)
  dimidy = NCDF_DIMDEF(netcdf_id, file1dinfo[1].name ,  file1dinfo[1].size)
  dimidt = NCDF_DIMDEF(netcdf_id, file1dinfo[2].name, /UNLIMITED)
;
; global attributes
  NCDF_ATTPUT, netcdf_id, 'Conventions', 'GDT 1.2', /GLOBAL ; ++ conformite
  NCDF_ATTPUT, netcdf_id, 'file_name'  , file3, /GLOBAL
  NCDF_ATTPUT, netcdf_id, 'Title'      , $
               'delta between '+ file1 + ' and ' + file2, /GLOBAL
;
; declaration of variables
; 4 common variables for the two files to produce
  varid = lonarr(4)
;
  varid[0] = NCDF_VARDEF(netcdf_id, file1vinfo[0].name  , [dimidx, dimidy], /FLOAT)
  NCDF_ATTPUT, netcdf_id, varid[0], $
               (*file1vatts[0])[0].name , $
               STRING(*(*file1vatts[0])[0].values)
  NCDF_ATTPUT, netcdf_id, varid[0], $
               (*file1vatts[0])[1].name, $
               *(*file1vatts[0])[1].values
  NCDF_ATTPUT, netcdf_id, varid[0], $
               (*file1vatts[0])[2].name, $
               *(*file1vatts[0])[2].values
  NCDF_ATTPUT, netcdf_id, varid[0], $
               (*file1vatts[0])[3].name, $
               STRING(*(*file1vatts[0])[3].values)
;
  varid[1] = NCDF_VARDEF(netcdf_id, file1vinfo[1].name  , [dimidx, dimidy], /FLOAT)
  NCDF_ATTPUT, netcdf_id, varid[1], $
               (*file1vatts[1])[0].name, $
               STRING(*(*file1vatts[1])[0].values)
  NCDF_ATTPUT, netcdf_id, varid[1], $
               (*file1vatts[1])[1].name, $
               *(*file1vatts[1])[1].values
  NCDF_ATTPUT, netcdf_id, varid[1], $
               (*file1vatts[1])[2].name, $
               *(*file1vatts[1])[2].values
  NCDF_ATTPUT, netcdf_id, varid[1], $
               (*file1vatts[1])[3].name, $
               STRING(*(*file1vatts[1])[3].values)
;
  varid[2] = NCDF_VARDEF(netcdf_id, file1vinfo[2].name, [dimidt], /FLOAT)
;
  varid[3] =  NCDF_VARDEF(netcdf_id, 'delta', [dimidx, dimidy, dimidt], /FLOAT)
  NCDF_ATTPUT, netcdf_id, varid[3], 'long_name', 'delta'
  NCDF_ATTPUT, netcdf_id, varid[3], $
               (*file1vatts[3])[0].name, $
               STRING(*(*file1vatts[3])[0].values)
;
  NCDF_ATTPUT, netcdf_id, varid[3], 'valid_min', minarr ,/FLOAT
  NCDF_ATTPUT, netcdf_id, varid[3], 'valid_max', maxarr ,/FLOAT
;
;---------------------------
; end of header definition, writing of the NetCdf files
;---------------------------
;
  NCDF_CONTROL, netcdf_id, /ENDEF
;
  NCDF_VARPUT, netcdf_id, file1vinfo[0].name, (*file1data[0])
  NCDF_VARPUT, netcdf_id, file1vinfo[1].name, (*file1data[1])
  NCDF_VARPUT, netcdf_id, file1vinfo[2].name, (*file1data[2])
;
; write the data
  NCDF_VARPUT, netcdf_id, 'delta', delta
;
;---------------------------
; close the netcdf files
  NCDF_CLOSE, netcdf_id
;
  msg = 'iii : ' + file3 + ' created'
  ras = report(msg)
;
END