Changeset 271 for trunk/SRC/ToBeReviewed


Ignore:
Timestamp:
08/30/07 14:44:23 (17 years ago)
Author:
smasson
Message:

bugfix for interpolation from ORCA2 without mask

Location:
trunk/SRC/ToBeReviewed
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/SRC/ToBeReviewed/HOPE/read_hope.pro

    r232 r271  
    452452  timename = strlowcase((tag_names(dimlist))[wathinside.recdim]) 
    453453; read the time axis in julina days 
    454   time = ncdf_timeget(cdfid, timename) 
     454  time = ncdf_gettime(filename, cdfid) 
    455455; update the dimlist structure 
    456456  dimlist.(wathinside.recdim) = time 

  • trunk/SRC/ToBeReviewed/INIT/initncdf.pro

    r238 r271  
    1111; A string giving the name of the NetCdf file 
    1212; 
    13 ; @keyword INVMASK {default=0}{type=scalar: 0 or 1} 
    14 ; Inverse the land/sea mask (that should have 0/1 values for land/sea): mask = 1-mask 
    15 ; 
    16 ; @keyword MASKNAME {type=string} 
    17 ; A string giving the name of the variable in the file 
    18 ; that contains the land/sea mask 
    19 ; 
    20 ; @keyword MISSING_VALUE {type=scalar} 
    21 ; To define (or redefine if the attribute is 
    22 ; already existing) the missing values used with USEASMASK 
    23 ; keyword 
    24 ; 
    2513; @keyword START1 {default=0}{type=scalar: 0 or 1} 
    2614; Index the axis from 1 instead of 0 when using 
    2715; /xyindex and/or /zindex 
    28 ; 
    29 ; @keyword USEASMASK {type=scalar string} 
    30 ; A string giving the name of the variable in the file 
    31 ; that will be used to build the land/sea mask. In this case the 
    32 ; mask is based on the first record (if record dimension 
    33 ; exists). The mask is build according to : 
    34 ;    1 the keyword missing_value if existing 
    35 ;    2 the attribute 'missing_value' if existing 
    36 ;    3 NaN values if existing 
    3716; 
    3817; @keyword ZAXISNAME {default='z', 'level', 'lev', 'depth...'}{type=scalar string} 
    … …  
    5231; 
    5332; @keyword _EXTRA 
    54 ; Used to pass keywords to <pro>computegrid</pro> and  
    55 ; <pro>ncdf_getaxis</pro> 
     33; Used to pass keywords to <pro>computegrid</pro>,  
     34; <pro>ncdf_getaxis</pro>, <pro>ncdf_getmask</pro> and <pro>isafile</pro> 
    5635; 
    5736; @uses 
    … …  
    7756; 
    7857PRO initncdf, ncfilein $ 
    79               , ZAXISNAME = zaxisname, MASKNAME = maskname $ 
    80               , INVMASK = invmask, USEASMASK = useasmask $ 
    81               , MISSING_VALUE = missing_value, START1 = start1 $ 
     58              , ZAXISNAME = zaxisname, START1 = start1 $ 
    8259              , XYINDEX = xyindex, ZINDEX = zindex $ 
    8360              , _EXTRA = ex 
    … …  
    10077; what is inside the file 
    10178  inside = ncdf_inquire(cdfid) 
    102 ;------------------------------------------------------------ 
    103 ; name of all dimensions 
    104   namedim = strarr(inside.ndims) 
    105   for dimiq = 0, inside.ndims-1 do begin 
    106     ncdf_diminq, cdfid, dimiq, tmpname, value 
    107     namedim[dimiq] = strlowcase(tmpname) 
    108   ENDFOR 
    10979;---------------------------------------------------------- 
    11080; name of the variables 
    … …  
    149119;---------------------------------------------------------- 
    150120; mask 
    151   IF NOT (keyword_set(maskname) OR keyword_set(useasmask)) AND keyword_set(romsgrid) THEN maskname = 'mask_rho' 
    152   CASE 1 OF 
    153     keyword_set(maskname):BEGIN 
    154       mskid = (where(namevar EQ strlowcase(maskname)))[0] 
    155       if mskid NE -1 THEN BEGIN 
    156         mskinq = ncdf_varinq(cdfid, mskid) 
    157 ; is the mask variable containing the record dimension? 
    158         withrcd = (where(mskinq.dim EQ inside.recdim))[0] 
    159         IF withrcd NE -1 THEN BEGIN 
    160 ; in order to read only the first record 
    161 ; we need to get the size of each dimension 
    162           count = replicate(1L, mskinq.ndims) 
    163           FOR d = 0, mskinq.ndims -1 DO BEGIN 
    164             IF d NE withrcd THEN BEGIN 
    165               ncdf_diminq, cdfid, mskinq.dim[d], name, size 
    166               count[d] = size 
    167             ENDIF 
    168           ENDFOR 
    169 ; read the variable for the first record 
    170           ncdf_varget, cdfid, mskid, tmask, count = count 
    171         ENDIF ELSE ncdf_varget, cdfid, mskid, tmask 
    172 ; check if we need to applay add_offset and scale factor 
    173         FOR a = 0, mskinq.natts-1 DO BEGIN 
    174           attname = ncdf_attname(cdfid, mskid, a) 
    175           CASE strlowcase(attname) OF 
    176             'add_offset':ncdf_attget, cdfid, mskid, attname, add_offset 
    177             'scale_factor':ncdf_attget, cdfid, mskid, attname, scale_factor 
    178             ELSE: 
    179           ENDCASE 
    180         ENDFOR 
    181         IF n_elements(scale_factor) NE 0 THEN tmask = tmask*scale_factor 
    182         IF n_elements(add_offset) NE 0 THEN tmask = tmask+add_offset 
    183         if keyword_set(invmask) then tmask = 1-tmask 
    184         tmask = byte(round(tmask)) 
    185       ENDIF ELSE tmask = -1 
    186     END 
    187 ;.................. 
    188     keyword_set(useasmask):BEGIN 
    189       mskid = (where(namevar EQ strlowcase(useasmask)))[0] 
    190       if mskid NE -1 THEN BEGIN 
    191         mskinq = ncdf_varinq(cdfid, mskid) 
    192 ; is the mask variable containing the record dimension? 
    193         withrcd = (where(mskinq.dim EQ inside.recdim))[0] 
    194         IF withrcd NE -1 THEN BEGIN 
    195 ; in order to read only the first record 
    196 ; we need to get the size of each dimension 
    197           count = replicate(1L, mskinq.ndims) 
    198           FOR d = 0, mskinq.ndims -1 DO BEGIN 
    199             IF d NE withrcd THEN BEGIN 
    200               ncdf_diminq, cdfid, mskinq.dim[d], name, size 
    201               count[d] = size 
    202             ENDIF 
    203           ENDFOR 
    204 ; read the variable for the first record 
    205           ncdf_varget, cdfid, mskid, tmask, count = count 
    206         ENDIF ELSE ncdf_varget, cdfid, mskid, tmask 
    207 ; check if we need to applay add_offset and scale factor 
    208         FOR a = 0, mskinq.natts-1 DO BEGIN 
    209           attname = ncdf_attname(cdfid, mskid, a) 
    210           CASE strlowcase(attname) OF 
    211             'add_offset':ncdf_attget, cdfid, mskid, attname, add_offset 
    212             'scale_factor':ncdf_attget, cdfid, mskid, attname, scale_factor 
    213             'missing_value':IF n_elements(missing_value) EQ 0 THEN $ 
    214                ncdf_attget, cdfid, mskid, attname, missing_value 
    215             ELSE: 
    216           ENDCASE 
    217         ENDFOR 
    218         IF n_elements(scale_factor) NE 0 THEN tmask = tmask*scale_factor 
    219         IF n_elements(add_offset) NE 0 THEN tmask = tmask+add_offset 
    220         IF n_elements(missing_value) NE 0 THEN BEGIN 
    221 ; we have to take care of the float accuracy 
    222           CASE 1 OF 
    223             missing_value GE 1.e6:tmask = tmask LT (missing_value-10) 
    224             missing_value LE -1.e6:tmask = tmask GT (missing_value-10) 
    225             abs(missing_value) LE 1.e-6:tmask = abs(tmask) GT 1.e-6 
    226             ELSE:tmask = tmask NE missing_value 
    227           ENDCASE 
    228           if keyword_set(invmask) then tmask = 1-tmask 
    229         ENDIF ELSE BEGIN 
    230           tmask = finite(tmask) 
    231           IF min(tmask) EQ 1 THEN BEGIN 
    232             ras = report( 'missing or nan values not found...') 
    233             tmask = -1 
    234           ENDIF 
    235         ENDELSE 
    236       ENDIF ELSE tmask = -1 
    237     END 
    238 ;.................. 
    239     ELSE:tmask  = -1 
    240   ENDCASE 
     121  tmask = ncdf_getmask(cdfid, _extra = ex) 
     122;---------------------------------------------------------- 
    241123; 
    242124  ncdf_close, cdfid 
    243125; 
    244 ; compute the grid 
     126;---------------------------------------------------------- 
     127; call compute the grid 
    245128  if NOT keyword_set(zaxis) then BEGIN 
    246129    computegrid, xaxis = xaxis, yaxis = yaxis $ 

  • trunk/SRC/ToBeReviewed/LECTURE/read_ncdf.pro

    r240 r271  
    3131; Useless, defined for compatibility 
    3232; 
     33; @keyword ADDSCL_BEFORE {default=0}{type=scalar: 0 or 1} 
     34; put 1 to apply add_offset ad scale factor on data before looking for 
     35; missing values 
     36; 
    3337; @keyword BOXZOOM 
    3438; Contain the boxzoom on which we have to do the reading 
    … …  
    4145; 
    4246; @keyword INIT {default=0}{type=scalar: 0 or 1} 
    43 ; To call automatically initncdf, filename and thus 
     47; To call automatically initncdf with filename as input argument and thus 
    4448; redefine all the grid parameters 
    4549; 
    … …  
    6266; but only the array referring to the field. 
    6367; 
    64 ; @keyword TIMEVAR {type=string} 
    65 ; It define the name of the variable that 
    66 ; contains the time axis. This keyword can be useful if there 
    67 ; is no unlimited dimension or if the time axis selected by default 
    68 ; (the first 1D array with unlimited dimension) is not the good one. 
    69 ; 
    7068; @keyword ZETAFILENAME {default=FILENAME}{type=string} 
    7169; For ROMS outputs. The filename of the file where zeta variable should be read 
    … …  
    7573; 
    7674; @keyword _EXTRA 
    77 ; Used to pass keywords 
     75; Used to pass keywords to <pro>isafile</pro>, <pro>initncdf</pro>, 
     76; <pro>ncdf_gettime</pro> and <pro>domdef</pro> 
    7877; 
    7978; @returns 
    … …  
    9594; 
    9695FUNCTION read_ncdf, name, beginning, ending, compatibility, BOXZOOM = boxzoom, FILENAME = filename $ 
    97                     , PARENTIN = parentin, TIMESTEP = timestep, TIMEVAR = timevar $ 
     96                    , PARENTIN = parentin, TIMESTEP = timestep, ADDSCL_BEFORE = addscl_before $ 
    9897                    , TOUT = tout, NOSTRUCT = nostruct, CONT_NOFILL = CONT_NOFILL, INIT = init $ 
    9998                    , GRID = grid, CALLITSELF = callitself $ 
    … …  
    115114;  print,filename 
    116115; is parent a valid widget ? 
    117   if keyword_set(parentin) then BEGIN 
     116  IF keyword_set(parentin) THEN BEGIN 
    118117    parent = long(parentin) 
    119118    parent = parent*widget_info(parent, /managed) 
    … …  
    123122; Opening of the name file 
    124123;------------------------------------------------------------ 
    125   if size(filename, /type) NE 7 then $ 
     124  IF size(filename, /type) NE 7 THEN $ 
    126125    return, report('read_ncdf cancelled') 
    127126  IF !version.OS_FAMILY EQ 'unix' THEN spawn, '\file '+filename+' > /dev/null' 
    128127  cdfid = ncdf_open(filename) 
    129   contient = ncdf_inquire(cdfid) 
     128  inq = ncdf_inquire(cdfid) 
    130129;------------------------------------------------------------ 
    131130; we check if the variable name exists in the file. 
    132131;------------------------------------------------------------ 
    133   if ncdf_varid(cdfid, name) EQ -1 then BEGIN 
     132  IF ncdf_varid(cdfid, name) EQ -1 THEN BEGIN 
    134133    ncdf_close, cdfid 
    135134    return, report('variable '+name+' !C not found in the file '+filename) 
    136135  ENDIF 
    137   varcontient = ncdf_varinq(cdfid, name) 
    138   IF varcontient.ndims LT 2 THEN return, report('read_ncdf cannot read scalar or 1D data') 
     136  varinq = ncdf_varinq(cdfid, name) 
     137  IF varinq.ndims LT 2 THEN return, report('read_ncdf cannot read scalar or 1D data') 
    139138; look for the dimension names 
    140   dimnames = strarr(varcontient.ndims) 
    141   FOR i = 0, varcontient.ndims-1 DO BEGIN 
    142     ncdf_diminq, cdfid, varcontient.dim[i], tmp, dimsize 
     139  dimnames = strarr(varinq.ndims) 
     140  FOR i = 0, varinq.ndims-1 DO BEGIN 
     141    ncdf_diminq, cdfid, varinq.dim[i], tmp, dimsize 
    143142    dimnames[i] = tmp 
    144143  ENDFOR 
    … …  
    146145; shall we redefine the grid parameters 
    147146;------------------------------------------------------------ 
    148   if keyword_set(init) THEN initncdf, filename, _extra = ex 
     147  IF keyword_set(init) THEN initncdf, filename, _extra = ex 
    149148;------------------------------------------------------------ 
    150149; check the time axis and the debut and ending dates 
    151150;------------------------------------------------------------ 
    152   if n_elements(beginning) EQ 0 then begin 
     151  IF n_elements(beginning) EQ 0 THEN BEGIN 
    153152    beginning = 0 
    154153    timestep = 1 
    155   endif 
    156   if keyword_set(timestep) then begin 
    157     firsttps = beginning[0] 
    158     if n_elements(ending) NE 0 then lasttps = ending[0] ELSE lasttps = firsttps 
    159     jpt = lasttps-firsttps+1 
    160     IF NOT keyword_set(callitself) then time = julday(1, 1, 1) + lindgen(jpt) 
    161   ENDIF ELSE BEGIN 
    162     if keyword_set(parent) then BEGIN 
     154  ENDIF 
     155; define time and jpt 
     156  CASE 1 OF 
     157    keyword_set(timestep):BEGIN  
     158      firsttps = beginning[0] 
     159      IF n_elements(ending) NE 0 THEN lasttps = ending[0] ELSE lasttps = firsttps 
     160      jpt = lasttps-firsttps+1 
     161      IF NOT keyword_set(callitself) then time = julday(1, 1, 1) + lindgen(jpt) 
     162    END 
     163    keyword_set(parent):BEGIN  
    163164      widget_control, parent, get_uvalue = top_uvalue 
    164165      filelist = extractatt(top_uvalue, 'filelist') 
    … …  
    168169      date1 = date2jul(beginning[0]) 
    169170      if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1 
    170       firsttps = where(time EQ date1) & firsttps = firsttps[0] 
    171       lasttps = where(time EQ date2) & lasttps = lasttps[0] 
    172     ENDIF ELSE BEGIN 
    173       IF keyword_set(timevar) THEN BEGIN 
    174         timeid = ncdf_varid(cdfid, timevar) 
    175         IF timeid EQ -1 THEN BEGIN 
    176           ncdf_close, cdfid 
    177           return, report('the file '+filename+' as no variable ' + timevar $ 
    178                          + '. !C Use the TIMESTEP keyword') 
    179         endif 
    180         timecontient = ncdf_varinq(cdfid, timeid) 
    181         contient.recdim = timecontient.dim[0] 
    182       ENDIF ELSE BEGIN 
    183 ; we find the infinite dimension 
    184         timedim = contient.recdim 
    185         if timedim EQ -1 then BEGIN 
    186           ncdf_close, cdfid 
    187           return, report('the file '+filename+' as no infinite dimension. !C Use TIMESTEP or TIMEVAR keyword') 
    188         endif 
    189 ; we find the FIRST time axis 
    190         timeid = 0 
    191         repeat BEGIN       ; As long as we have not find a variable having only one dimension: the infinite one 
    192           timecontient = ncdf_varinq(cdfid, timeid) ; that the variable contain. 
    193           timeid = timeid+1 
    194         endrep until (n_elements(timecontient.dim) EQ 1 $ 
    195                       AND timecontient.dim[0] EQ contient.recdim) $ 
    196           OR timeid EQ contient.nvars+1 
    197 ; 
    198         if timeid EQ contient.nvars+1 then BEGIN 
    199           ncdf_close, cdfid 
    200           return, report('the file '+filename+' as no time axis variable. !C Use the TIMESTEP keyword') 
    201         endif 
    202         timeid = timeid-1 
    203       ENDELSE 
    204 ; we must found the time origin of the julian calendar used in the 
    205 ; time axis. 
    206 ; does the attribut units an dcalendar exist for the variable time axis? 
    207       if timecontient.natts EQ 0 then BEGIN 
     171      firsttps = (where(abs(time - date1) LT 0.9d/86400.d))[0]   ; beware of rounding errors... 
     172      lasttps = (where(abs(time - date2) LT 0.9d/86400.d))[0] 
     173      jpt = lasttps-firsttps+1 
     174    END 
     175    ELSE:BEGIN  
     176      time = ncdf_gettime(filename, cdfid, caller = 'read_ncdf', err = err, _extra = ex) 
     177      IF n_elements(err) NE 0 THEN BEGIN 
     178        dummy = report(err) 
    208179        ncdf_close, cdfid 
    209         return, report('the variable '+timecontient.name+' has no attribut.!C Use the TIMESTEP keyword or add the attribut ''units'' to the variable') 
    210       endif 
    211       attnames = strarr(timecontient.natts) 
    212       for attiq = 0, timecontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, timeid, attiq) 
    213       if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN 
    214         ncdf_close, cdfid 
    215         return, report('Attribut ''units'' not found for the variable '+timecontient.name+' !C Use the TIMESTEP keyword') 
     180        return, -1 
    216181      ENDIF 
    217 ; 
    218 ; now we try to find the attribut called calendar... 
    219 ; the attribute "calendar" exists? 
    220 ; If no, we suppose that the calendar is gregorian calendar 
    221 ; 
    222       if (where(attnames EQ 'calendar'))[0] NE -1 then BEGIN 
    223         ncdf_attget, cdfid, timeid, 'calendar', value 
    224         value = string(value) 
    225         CASE value OF 
    226           'noleap':key_caltype = 'noleap' 
    227           '360d':key_caltype = '360d' 
    228           'greg':IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg' 
    229           ELSE:BEGIN 
    230 ;            notused = report('Unknown calendar: '+value+', we use greg calendar.') 
    231             key_caltype = 'greg' 
    232           END 
    233         ENDCASE 
    234       ENDIF ELSE BEGIN 
    235 ;        notused = report('Unknown calendar, we use '+key_caltype+' calendar.') 
    236         IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg' 
    237       ENDELSE 
    238 ; 
    239 ; now we take acre of units attribut 
    240       ncdf_attget, cdfid, timeid, 'units', value 
    241 ; 
    242 ; time_counter:units = "seconds since 0001-01-01 00:00:00" ; 
    243 ; time_counter:units = "hours since 0001-01-01 00:00:00" ; 
    244 ; time_counter:units = "days since 1979-01-01 00:00:00" ; 
    245 ; time_counter:units = "months since 1979-01-01 00:00:00" ; 
    246 ; time_counter:units = "years since 1979-01-01 00:00:00" ; 
    247 ; 
    248 ; we decript the "units" attribut to find the time origin 
    249       value = strtrim(strcompress(string(value)), 2) 
    250       mots = str_sep(value, ' ') 
    251       unite = mots[0] 
    252       unite = strlowcase(unite) 
    253       IF strpos(unite, 's', strlen(unite)-1) NE -1 THEN unite = strmid(unite, 0, strlen(unite)-1) 
    254       IF strpos(unite, 'julian_') NE -1 THEN unite = strmid(unite, 7) 
    255       IF unite NE 'second' AND unite NE 'hour' AND unite NE 'day' $ 
    256          AND unite NE 'month' AND unite NE 'year' THEN BEGIN 
    257         ncdf_close, cdfid 
    258         return, report('time units does not start with seconds/hours/days/months/years') 
    259       ENDIF 
    260       IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN 
    261         ncdf_close, cdfid 
    262         return, report('attribut units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*') 
    263       ENDIF 
    264       depart = str_sep(mots[2], '-') 
    265       ncdf_varget, cdfid, timeid, time 
    266       time = double(time) 
    267       case unite of 
    268         'second':time = julday(depart[1], depart[2], depart[0])+time/86400.d 
    269         'hour':time = julday(depart[1], depart[2], depart[0])+time/24.d 
    270         'day':time = julday(depart[1], depart[2], depart[0])+time 
    271         'month':BEGIN 
    272           if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m 
    273             time = julday(depart[1], depart[2], depart[0])+round(time*30) $ 
    274             ELSE for t = 0, n_elements(time)-1 DO $ 
    275             time[t] = julday(depart[1]+time[t], depart[2], depart[0]) 
    276         END 
    277         'year':BEGIN 
    278           if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y 
    279             time = julday(depart[1], depart[2], depart[0])+round(time*365) $ 
    280             ELSE for t = 0, n_elements(time)-1 do $ 
    281             time[t] = julday(depart[1], depart[2], depart[0]+time[t]) 
    282         END 
    283         ELSE:BEGIN 
    284           ncdf_close, cdfid 
    285           return, report('The "units" attribute of the time axis must be something like: !C "seconds since 0001-01-01 ..." !C "days since 1979-01-01 ..." !C "months since 1979-01-01 ..." !C "years since 1979-01-01 ..."') 
    286         end 
    287       ENDCASE 
     182; date1 
    288183      date1 = date2jul(beginning[0]) 
     184; date2 
    289185      if n_elements(ending) NE 0 then date2 = date2jul(ending[0]) ELSE date2 = date1 
    290       time = double(time) 
    291       firsttps = where(time GE date1) & firsttps = firsttps[0] 
     186; firsttps 
     187      firsttps = where(time GE (date1 - 0.9d/86400.d)) & firsttps = firsttps[0] 
    292188      if firsttps EQ -1 THEN BEGIN 
    293189        ncdf_close, cdfid 
    294190        return, report('date 1: '+strtrim(jul2date(date1), 1)+' is not found in the time axis.') 
    295191      ENDIF 
    296       lasttps = where(time LE date2) 
    297       if lasttps[0] EQ -1 THEN BEGIN 
     192; lasttps 
     193      lasttps = where(time LE (date2 + 0.9d/86400.d)) & lasttps = lasttps[n_elements(lasttps)-1] 
     194      if lasttps EQ -1 THEN BEGIN 
    298195        ncdf_close, cdfid 
    299196        return, report('the time axis has no date before date 2: '+strtrim(jul2date(date2), 1)) 
    300197      endif 
    301       lasttps = lasttps[n_elements(lasttps)-1] 
    302198      if lasttps LT firsttps then BEGIN 
    303199        ncdf_close, cdfid 
    304200        return, report('the time axis has no dates between date1 and  date 2: '+strtrim(jul2date(date1), 1)+' '+strtrim(jul2date(date2), 1)) 
    305201      endif 
    306     ENDELSE 
    307     time = time[firsttps:lasttps] 
    308     jpt = lasttps-firsttps+1 
    309   ENDELSE 
     202      time = time[firsttps:lasttps] 
     203      jpt = lasttps-firsttps+1 
     204    END 
     205  ENDCASE 
    310206;------------------------------------------------------------ 
    311207; Name of the grid on which the field refer to. 
    … …  
    316212; are we in one of the case corresponding to ROMS conventions? 
    317213      CASE 1 OF 
    318         dimnames[2 <(varcontient.ndims-1)] EQ 's_w':vargrid = 'W' 
     214        dimnames[2 <(varinq.ndims-1)] EQ 's_w':vargrid = 'W' 
    319215        dimnames[0] EQ 'xi_rho' AND dimnames[1] EQ 'eta_rho':vargrid = 'T' 
    320216        dimnames[0] EQ 'xi_u'   AND dimnames[1] EQ 'eta_u'  :vargrid = 'U' 
    … …  
    342238  ENDELSE 
    343239;--------------------------------------------------------------- 
    344 ; call the init function ? 
    345 ;--------------------------------------------------------------- 
    346  
    347 ;--------------------------------------------------------------- 
    348240; redefinition of the  domain 
    349241;--------------------------------------------------------------- 
    350   if keyword_set(tout) then begin 
     242  IF keyword_set(tout) THEN BEGIN 
    351243    nx = jpi 
    352244    ny = jpj 
    … …  
    413305;--------------------------------------------------------------------- 
    414306;--------------------------------------------------------------------- 
    415 ; We define global variable joined with the variable. 
     307; We define common (cm_4data) variables associated with the variable. 
    416308;--------------------------------------------------------------------- 
    417309; varname 
    418310  IF NOT keyword_set(callitself) THEN varname = name 
    419 ; varunit 
    420   if varcontient.natts NE 0 then begin 
    421     attnames = strarr(varcontient.natts) 
    422     for attiq = 0, varcontient.natts-1 do attnames[attiq] = ncdf_attname(cdfid, name, attiq) 
    423     lowattnames = strlowcase(attnames) 
    424 ; 
    425     found = (where(lowattnames EQ 'units'))[0] 
    426     IF found NE -1 then ncdf_attget, cdfid, name, attnames[found], value ELSE value = '' 
    427     IF NOT keyword_set(callitself) THEN varunit = strtrim(string(value), 2) 
    428 ; 
    429     found = (where(lowattnames EQ 'add_offset'))[0] 
    430     if found NE -1 then ncdf_attget, cdfid, name, attnames[found], add_offset ELSE add_offset = 0. 
    431 ; 
    432     found = (where(lowattnames EQ 'scale_factor'))[0] 
    433     if found NE -1 then ncdf_attget, cdfid, name, attnames[found], scale_factor ELSE scale_factor = 1. 
    434 ; 
    435     missing_value = 'no' 
    436     found = (where(lowattnames EQ '_fillvalue'))[0] 
    437     if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value 
    438     found = (where(lowattnames EQ 'missing_value'))[0] 
    439     if found NE -1 then ncdf_attget, cdfid, name, attnames[found], missing_value 
    440 ; 
    441   ENDIF ELSE BEGIN 
    442     IF NOT keyword_set(callitself) THEN varunit = '' 
    443     add_offset = 0. 
    444     scale_factor = 1. 
    445     missing_value = 'no' 
    446   ENDELSE 
     311; varunit and others... 
     312  ncdf_getatt, cdfid, name, add_offset = add_offset, scale_factor = scale_factor, missing_value = missing_value, units = units 
     313  IF NOT keyword_set(callitself) THEN varunit = units 
    447314; vardate 
    448315; We make a legible date in function of the specified language. 
    … …  
    455322; We apply the value valmask on land points. 
    456323  if NOT keyword_set(cont_nofill) then begin 
    457     valmask = 1e20 
     324    valmask = 1.e20 
    458325    case 1 of 
    459       varcontient.ndims eq 2:BEGIN ;xy array 
    460         mask = mask[*, *, 0] 
     326      varinq.ndims eq 2:BEGIN ;xy array 
     327        earth = where(mask[*, *, 0] EQ 0) 
     328      END 
     329      varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] EQ -1:BEGIN ;xyz array 
    461330        earth = where(mask EQ 0) 
    462331      END 
    463       varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] EQ -1:BEGIN ;xyz array 
    464         earth = where(mask EQ 0) 
    465       END 
    466       varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] NE -1:BEGIN ;xyt array 
    467         mask = mask[*, *, 0] 
    468         earth = where(mask EQ 0) 
     332      varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1:BEGIN ;xyt array     
     333        earth = where(mask[*, *, 0] EQ 0) 
    469334        if earth[0] NE -1 then BEGIN 
    470335          earth = earth#replicate(1, jpt)+replicate(nx*ny, n_elements(earth))#lindgen(jpt) 
    471336        END 
    472337      END 
    473       varcontient.ndims eq 4:BEGIN ;xyzt array 
     338      varinq.ndims eq 4:BEGIN ;xyzt array 
    474339        earth = where(mask EQ 0) 
    475340        if earth[0] NE -1 then BEGIN 
    … …  
    477342        END 
    478343      END 
    479     endcase 
     344    ENDCASE 
    480345  ENDIF ELSE earth = -1 
    481346; we look for  missing_value 
     347; we apply add_offset, scale_factor and missing_value 
     348  IF keyword_set(addscl_before) THEN BEGIN 
     349    if scale_factor NE 1 then res = temporary(res)*scale_factor 
     350    if add_offset NE 0 then res = temporary(res)+add_offset 
     351  ENDIF 
    482352  IF size(missing_value, /type) NE 7 then BEGIN 
    483     IF size(missing_value, /type) EQ 1 THEN BEGIN 
    484       missing_value = strlowcase(string(missing_value)) 
    485       IF strmid(missing_value, 0, 1, /reverse_offset) EQ 'f' THEN $ 
    486         missing_value = strmid(missing_value, 0, strlen(missing_value)-1) 
    487       IF isnumber(string(missing_value), tmp) EQ 1 THEN missing_value = tmp ELSE BEGIN 
    488         ras = report('Warning: missing value is not a number: ' + missing_value) 
    489         missing_value = - 1 
    490       ENDELSE 
    491     ENDIF 
    492 ;    if missing_value NE valmask then begin 
    493     if abs(missing_value) LT 1e6 then missing = where(res EQ missing_value) $ 
    494     ELSE missing = where(abs(res) gt abs(missing_value)/10.) 
    495 ;    ENDIF ELSE missing = -1 
     353    IF finite(missing_value) EQ 1 THEN BEGIN 
     354        CASE 1 OF 
     355          missing_value GT 1.e6:missing = where(res GT missing_value-10.) 
     356          missing_value LT -1.e6:missing = where(res LT missing_value+10.) 
     357          abs(missing_value) LT 1.e-6:missing = where(res LT 1.e-6) 
     358          ELSE:missing = where(res EQ missing_value) 
     359        ENDCASE 
     360    ENDIF ELSE missing = where(finite(res) EQ 0) 
    496361  ENDIF ELSE missing = -1 
    497 ; we apply add_offset, scale_factor and missing_value 
    498   if scale_factor NE 1 then res = temporary(res)*scale_factor 
    499   if add_offset NE 0 then res = temporary(res)+add_offset 
     362  IF NOT keyword_set(addscl_before) THEN BEGIN 
     363    if scale_factor NE 1 then res = temporary(res)*scale_factor 
     364    if add_offset NE 0 then res = temporary(res)+add_offset 
     365  ENDIF 
    500366  if missing[0] NE -1 then res[temporary(missing)] = !values.f_nan 
    501367  if earth[0] NE -1 then res[temporary(earth)] = 1.e20 
    … …  
    508374      ncdf_attget, cdfid, 'hc', hc, /global 
    509375; look for all variables names 
    510       allvarnames = strarr(contient.nvars) 
    511       FOR i = 0, contient.nvars-1 DO BEGIN 
     376      allvarnames = strarr(inq.nvars) 
     377      FOR i = 0, inq.nvars-1 DO BEGIN 
    512378        tmp = ncdf_varinq( cdfid, i) 
    513379        allvarnames[i] = tmp.name 

  • trunk/SRC/ToBeReviewed/LECTURE/read_ncdf_varget.pro

    r231 r271  
    9797;...................................................................... 
    9898;...................................................................... 
    99    varcontient.ndims eq 2:BEGIN ;xy array 
     99   varinq.ndims eq 2:BEGIN ;xy array 
    100100;...................................................................... 
    101101;...................................................................... 
    … …  
    185185;...................................................................... 
    186186;...................................................................... 
    187    varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] EQ -1:BEGIN ;xyz array 
     187   varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] EQ -1:BEGIN ;xyz array 
    188188;...................................................................... 
    189189;...................................................................... 
    … …  
    276276;...................................................................... 
    277277;...................................................................... 
    278    varcontient.ndims eq 3 AND (where(varcontient.dim EQ contient.recdim))[0] NE -1:BEGIN ;xyt array 
     278   varinq.ndims eq 3 AND (where(varinq.dim EQ inq.recdim))[0] NE -1:BEGIN ;xyt array 
    279279;...................................................................... 
    280280;...................................................................... 
    … …  
    376376;...................................................................... 
    377377;...................................................................... 
    378    varcontient.ndims eq 4:BEGIN ;xyzt array 
     378   varinq.ndims eq 4:BEGIN ;xyzt array 
    379379;...................................................................... 
    380380;...................................................................... 
    … …  
    481481; we apply reverse 
    482482  IF keyword_set(key_yreverse) AND ny NE 1 THEN BEGIN 
    483     IF varcontient.ndims - ((where(varcontient.dim EQ contient.recdim))[0] NE -1) EQ 2 THEN $ 
     483    IF varinq.ndims - ((where(varinq.dim EQ inq.recdim))[0] NE -1) EQ 2 THEN $ 
    484484       res = reverse(reform(res, nx, ny, jpt, /overwrite),  2) $ 
    485485    ELSE res = reverse(reform(res, nx, ny, nz, jpt, /overwrite),  2) 
    486486  ENDIF 
    487487  if keyword_set(key_zreverse) AND nz NE 1 $ 
    488      AND varcontient.ndims - ((where(varcontient.dim EQ contient.recdim))[0] NE -1) EQ 3 THEN $ 
     488     AND varinq.ndims - ((where(varinq.dim EQ inq.recdim))[0] NE -1) EQ 3 THEN $ 
    489489        res = reverse(reform(res, nx, ny, nz, jpt, /overwrite),  3) 
    490490; 

  • trunk/SRC/ToBeReviewed/WIDGET/AUTOUR_de_XXX/scanfile.pro

    r238 r271  
    1111; 
    1212; @keyword _EXTRA 
    13 ; Used to pass keywords to <pro>isafile</pro>  
    14 ; and <pro>ncdf_getaxis</pro> 
     13; Used to pass keywords to <pro>isafile</pro>,  
     14; <pro>ncdf_getaxis</pro> and <pro>ncdf_gettime</pro> 
    1515; 
    1616; @returns 
    … …  
    8181; What contains the file? 
    8282;------------------------------------------------------------ 
    83   inside = ncdf_inquire(cdfid)  ; 
     83  inq = ncdf_inquire(cdfid)  ; 
    8484;------------------------------------------------------------ 
    8585; name of all dimensions 
    8686;------------------------------------------------------------ 
    87   namedim = strarr(inside.ndims) 
    88   for dimiq = 0, inside.ndims-1 do begin 
     87  namedim = strarr(inq.ndims) 
     88  for dimiq = 0, inq.ndims-1 do begin 
    8989    ncdf_diminq, cdfid, dimiq, tmpname, value 
    9090    namedim[dimiq] = strlowcase(tmpname) 
    … …  
    9898;------------------------------------------------------------ 
    9999; we keep only the variables containing at least x, y and time dimension (if existing...) 
    100   namevar = strarr(inside.nvars) 
    101   for varid = 0, inside.nvars-1 do begin 
    102     invar = ncdf_varinq(cdfid, varid) ; what contains the variable? 
    103     if (inter(invar.dim, dimidx))[0] NE -1 AND $ 
    104        (inter(invar.dim, dimidy))[0] NE -1 AND $ 
    105        ((where(invar.dim EQ inside.recdim))[0] NE -1 OR inside.recdim EQ -1) $ 
    106     THEN namevar[varid] = invar.name 
     100  namevar = strarr(inq.nvars) 
     101  for varid = 0, inq.nvars-1 do begin 
     102    varinq = ncdf_varinq(cdfid, varid) ; what contains the variable? 
     103    if (inter(varinq.dim, dimidx))[0] NE -1 AND $ 
     104       (inter(varinq.dim, dimidy))[0] NE -1 AND $ 
     105       ((where(varinq.dim EQ inq.recdim))[0] NE -1 OR inq.recdim EQ -1) $ 
     106    THEN namevar[varid] = varinq.name 
    107107  ENDFOR 
    108108  namevar = namevar[where(namevar NE '')] 
    … …  
    117117; for each variable, look if we in one of the case corresponding to ROMS conventions? 
    118118    FOR i = 0, n_elements(namevar)-1 do begin 
    119       invar = ncdf_varinq(cdfid, namevar[i]) 
    120       tmpnm = namedim[invar.dim] 
     119      varinq = ncdf_varinq(cdfid, namevar[i]) 
     120      tmpnm = namedim[varinq.dim] 
    121121; are we in one of the case corresponding to ROMS conventions? 
    122122      CASE 1 OF 
    123         tmpnm[2 <(invar.ndims-1)] EQ 's_w':vargrid = 'W' 
     123        tmpnm[2 < (varinq.ndims-1)] EQ 's_w':vargrid = 'W' 
    124124        tmpnm[0] EQ 'xi_rho' AND tmpnm[1] EQ 'eta_rho':listgrid[i] = 'T' 
    125125        tmpnm[0] EQ 'xi_u'   AND tmpnm[1] EQ 'eta_u'  :listgrid[i] = 'U' 
    … …  
    152152; time axis 
    153153;------------------------------------------------------------ 
     154  time = ncdf_gettime(fullname, cdfid, caller = 'scanfile', err = err, _extra = ex) 
     155  IF n_elements(err) NE 0 THEN BEGIN 
     156    dummy = report(err) 
     157    jpt = abs(time) 
     158    fakecal = 1 
     159  ENDIF ELSE jpt = n_elements(time)  
     160; high frequency calendar: more than one element per day 
     161  IF max(histogram([long(time-time[0])])) GT 1 THEN fakecal = 1 ELSE fakecal = 0 
    154162  date0fk = date2jul(19000101) 
    155   IF inside.recdim EQ -1 THEN BEGIN 
    156     jpt = 1 
    157     time = date0fk 
    158     fakecal = 1 
    159   ENDIF ELSE BEGIN 
    160     ncdf_diminq, cdfid, inside.recdim, timedimname, jpt 
    161 ; we look for the variable containing the time axis 
    162 ; we look for the first variable having for only dimension inside.recdim 
    163     varid = 0 
    164     repeat BEGIN 
    165       IF varid LT inside.nvars THEN BEGIN 
    166         invar = ncdf_varinq(cdfid, varid) 
    167         varid = varid+1 
    168       ENDIF ELSE varid = 0 
    169     endrep until (n_elements(invar.dim) EQ 1 AND invar.dim[0] EQ inside.recdim) OR (varid EQ 0) 
    170     varid = varid-1 
    171 ; 
    172     CASE 1 OF 
    173       varid EQ -1:BEGIN 
    174         dummy = report('the file '+fullname+' has no time axis.!C we create a fake calendar ...') 
    175         fakecal = 1 
    176         time = date0fk + lindgen(1>jpt) 
    177       END 
    178       invar.natts EQ 0:BEGIN 
    179         dummy = report('the variable '+invar.name+' has no attribut.!C we create a fake calendar ...') 
    180         fakecal = 1 
    181         time = date0fk + lindgen(1>jpt) 
    182       END 
    183       ELSE:BEGIN 
    184 ; 
    185 ; we want to know which attributes are attached to the time variable... 
    186 ; 
    187         attnames = strarr(invar.natts) 
    188         for attiq = 0, invar.natts-1 do attnames[attiq] = ncdf_attname(cdfid, varid, attiq) 
    189         if (where(attnames EQ 'units'))[0] EQ -1 then BEGIN 
    190           dummy = report('Attribut ''units'' not found for the variable '+invar.name+'!C we create a fake calendar ...') 
    191           fakecal = 1 
    192           time = date0fk + lindgen(1>jpt) 
    193         ENDIF ELSE BEGIN 
    194 ; we read the time axis 
    195           ncdf_varget, cdfid, varid, time 
    196           time = double(time) 
    197           ncdf_attget, cdfid, varid, 'units', value 
    198 ; time_counter:units = "seconds since 0001-01-01 00:00:00" ; 
    199 ; time_counter:units = "hours since 0001-01-01 00:00:00" ; 
    200 ; time_counter:units = "days since 1979-01-01 00:00:00" ; 
    201 ; time_counter:units = "months since 1979-01-01 00:00:00" ; 
    202 ; time_counter:units = "years since 1979-01-01 00:00:00" ; 
    203           value = strtrim(strcompress(string(value)), 2) 
    204           mots = str_sep(value, ' ') 
    205           unite = mots[0] 
    206           unite = strlowcase(unite) 
    207           IF strpos(unite, 's', strlen(unite)-1) NE -1 THEN unite = strmid(unite, 0, strlen(unite)-1) 
    208           IF strpos(unite, 'julian_') NE -1 THEN unite = strmid(unite, 7) 
    209           err = 0 
    210           IF unite NE 'second' AND unite NE 'hour' AND unite NE 'day' $ 
    211              AND unite NE 'month' AND unite NE 'year' THEN BEGIN 
    212             dummy = report('time units does not start with seconds/hours/days/months/years') 
    213             err = 1 
    214           ENDIF 
    215           IF stregex(value, '[^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*', /boolean) EQ 0 THEN BEGIN 
    216             dummy = report('attribut units of time has not the good format: [^ ]* since ([0-9]){1,4}-([0-9]){1,2}-([0-9]){1,2}.*') 
    217             err = 1 
    218           ENDIF 
    219           IF err GT 0 THEN BEGIN 
    220             fakecal = 1 
    221             time = date0fk + lindgen(1>jpt) 
    222           ENDIF ELSE BEGIN 
    223             debut = str_sep(mots[2], '-') 
    224 ; 
    225 ; now we try to find the attribut called calendar... 
    226 ; the attribute "calendar" exists? 
    227 ; If no, we suppose that the calendar is gregorian calendar 
    228 ; 
    229             if (where(attnames EQ 'calendar'))[0] NE -1 then BEGIN 
    230               ncdf_attget, cdfid, varid, 'calendar', value 
    231               value = string(value) 
    232               CASE value OF 
    233                 'noleap':key_caltype = 'noleap' 
    234                 '360d':key_caltype = '360d' 
    235                 'greg':IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg' 
    236                 ELSE:BEGIN 
    237 ;            notused = report('Unknown calendar: '+value+', we use greg calendar.') 
    238                   key_caltype = 'greg' 
    239                 END 
    240               ENDCASE 
    241             ENDIF ELSE BEGIN 
    242 ;        notused = report('Unknown calendar, we use '+key_caltype+' calendar.') 
    243               IF n_elements(key_caltype) EQ 0 THEN key_caltype = 'greg' 
    244             ENDELSE 
    245 ; 
    246 ; BEWARE we have to get back the calendar attribute and ajust time by consequence... 
    247 ; 
    248 ; 
    249 ; We pass time in IDL julian days 
    250 ; 
    251             case unite of 
    252               'second':time = julday(debut[1], debut[2], debut[0])+time/86400.d 
    253               'hour':time = julday(debut[1], debut[2], debut[0])+time/24.d 
    254               'day':time = julday(debut[1], debut[2], debut[0])+time 
    255               'month':BEGIN 
    256                 if total(fix(time) NE time) NE 0 then $ ; we switch to days with 30d/m 
    257                    time = julday(debut[1], debut[2], debut[0])+round(time*30) $ 
    258                 ELSE for t = 0, n_elements(time)-1 DO $ 
    259                    time[t] = julday(debut[1]+time[t], debut[2], debut[0]) 
    260               END 
    261               'year':BEGIN 
    262                 if total(fix(time) NE time) NE 0 then $ ; we switch to days with 365d/y 
    263                    time = julday(debut[1], debut[2], debut[0])+round(time*365) $ 
    264                 ELSE for t = 0, n_elements(time)-1 do $ 
    265                    time[t] = julday(debut[1], debut[2], debut[0]+time[t]) 
    266               END 
    267             ENDCASE 
    268 ; 
    269 ; high frequency calendar: more than one element per day 
    270             IF max(histogram([long(time-time[0])])) GT 1 THEN fakecal = 1 ELSE fakecal = 0 
    271             date0fk = date2jul(19000101) 
    272             IF keyword_set(fakecal) THEN time = date0fk+lindgen(1>jpt) $ 
    273             ELSE time = long(time) 
    274 ; 
    275           ENDELSE 
    276         ENDELSE 
    277       END 
    278     ENDCASE 
    279   ENDELSE 
     163  IF keyword_set(fakecal) THEN time = date0fk+lindgen(1 > jpt) 
    280164;------------------------------------------------------------ 
    281165  ncdf_close, cdfid 

  • trunk/SRC/ToBeReviewed/WIDGET/COMPOUND_WIDGET/cw_calendar.pro

    r262 r271  
    6868; 
    6969@cm_4cal 
     70 
    7071; get back the calendar and its related informations 
    7172  winfo_id = widget_info(id, find_by_uname = 'infocal') 
    … …  
    9293    day = value MOD 100L 
    9394; check that the date exists in the calendar 
    94   if (where(infowid.calendar EQ julday(month, day, year)))[0] EQ - 1 then return 
     95  if (where(abs(infowid.calendar - julday(month, day, year)) LT 1./86400.))[0] EQ - 1 then return 
    9596; update the value of infocal 
    9697    infowid.date = value 
Note: See TracChangeset for help on using the changeset viewer.