Changeset 271


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

bugfix for interpolation from ORCA2 without mask

Location:
trunk/SRC
Files:
1 deleted
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro

    r257 r271  
    4040; @restrictions 
    4141;  -  the input grid must be an "irregular 2D grid", defined as a grid made 
    42 ;  of quadrilateral cells which corners positions are defined with olonin and olat 
     42;  of quadrilateral cells 
    4343;  -  We supposed the data are located on a sphere, with a periodicity along 
    4444;  the longitude 
    … …  
    7171  IF n_elements(omsk) EQ 1 AND omsk[0] EQ -1 THEN BEGIN 
    7272    omsk = replicate(1b, jpio, jpjo) 
    73    IF jpio EQ 182 AND jpjo EQ 149 THEN BEGIN 
    74       lontmp = (olonin[1:180, *] + 3600) MOD 360 
    75       lattmp = olat[1:180, *] 
    76       bad = abs(abs(lontmp - shift(lontmp, 1, 0)) - 180) LT 176 $ 
    77             OR  abs(abs(lontmp - shift(lontmp, 0, 1)) - 180) LT 176 $ 
    78             OR  abs(abs(lontmp - shift(lontmp, 1, 1)) - 180) LT 176 $ 
    79             OR  abs(abs(lattmp - shift(lattmp, 1, 0)) - 180) LT 176 $ 
    80             OR  abs(abs(lattmp - shift(lattmp, 0, 1)) - 180) LT 176 $ 
    81             OR  abs(abs(lattmp - shift(lattmp, 1, 1)) - 180) LT 176 
    82       bad = bad AND lattmp LT 50 
    83       omsk[1:180, 1:148] = 1b - bad[*, 1:148] 
    84       omsk[0, *] = omsk[180, *] 
    85       omsk[181, *] = omsk[1, *] 
    86     ENDIF ELSE omsk = replicate(1b, jpio, jpjo) 
    87   ENDIF 
    88   IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja) 
    89   IF n_elements(omsk) NE jpio*jpjo THEN BEGIN 
    90     ras = report('input grid mask do not have the good size') 
    91     stop 
    92   ENDIF 
    93   IF n_elements(amsk) NE jpia*jpja THEN BEGIN 
    94     ras = report('output grid mask do not have the good size') 
    95     stop 
    96   ENDIF 
     73; if this is ORCA2 grid... 
     74    IF (jpio EQ 180 OR jpio EQ 182) AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ) THEN BEGIN 
     75; we look for ill defined cells. 
     76      IF jpio EQ 182 THEN lontmp = olonin[1:180, *] ELSE lontmp = olonin 
     77; longitudinal size of the cells... 
     78      a = (lontmp-shift(lontmp, 1, 0)) 
     79      d1 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) 
     80      a = (lontmp-shift(lontmp, 1, 1)) 
     81      d2 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) 
     82      a = (shift(lontmp, 0, 1)-shift(lontmp, 1, 0)) 
     83      d3 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) 
     84      a = (shift(lontmp, 0, 1)-shift(lontmp, 1, 1)) 
     85      d4 = 0.5 * max( [[[a]], [[360+a]]] MOD 360, dimension = 3) 
     86      md = [[[d1]], [[d2]], [[d3]], [[d4]]] 
     87      bad = max(md, dimension = 3) GE 1.5 
     88      bad = (bad + shift(bad, -1, -1) + shift(bad, 0, -1) + shift(bad, -1, 0)) < 1 
     89      IF jpio EQ 182 THEN BEGIN 
     90        omsk[1:180, 80:120] = 1b - bad[*, 80:120] 
     91        omsk[0, *] = omsk[180, *] 
     92        omsk[181, *] = omsk[1, *] 
     93      ENDIF ELSE omsk[*, 80:120] = 1b - bad[*, 80:120] 
     94    ENDIF 
     95  ENDIF ELSE BEGIN  
     96    IF n_elements(omsk) NE jpio*jpjo THEN BEGIN 
     97      ras = report('input grid mask do not have the good size') 
     98      stop 
     99    ENDIF 
     100  ENDELSE  
     101  IF n_elements(amsk) EQ 1 AND amsk[0] EQ -1 THEN amsk = replicate(1b, jpia, jpja) ELSE BEGIN  
     102    IF n_elements(amsk) NE jpia*jpja THEN BEGIN 
     103      ras = report('output grid mask do not have the good size') 
     104      stop 
     105    ENDIF 
     106  ENDELSE  
    97107; 
    98108; longitude, between 0 and 360 
    … …  
    175185; ORCA cases : orca grid is irregular only northward of 40N 
    176186        CASE 1 OF 
    177           jpio EQ 92  AND (jpjo EQ 76   OR jpjo EQ 75   OR jpjo EQ 74  ):onsphere = yy GT 40 
    178           jpio EQ 180 AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ):onsphere = yy GT 40 
    179           jpio EQ 720  AND (jpjo EQ 522  OR jpjo EQ 521  OR jpjo EQ 520 ):onsphere = yy GT 40 
    180           jpio EQ 1440 AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40 
     187          (jpio EQ 90 OR jpio EQ 92) AND (jpjo EQ 76   OR jpjo EQ 75   OR jpjo EQ 74  ):onsphere = yy GT 40 
     188          (jpio EQ 180 OR jpio EQ 182) AND (jpjo EQ 149  OR jpjo EQ 148  OR jpjo EQ 147 ):onsphere = yy GT 40 
     189          (jpio EQ 720 OR jpio EQ 722)  AND (jpjo EQ 522  OR jpjo EQ 521  OR jpjo EQ 520 ):onsphere = yy GT 40 
     190          (jpio EQ 1440 OR jpio EQ 1442) AND (jpjo EQ 1021 OR jpjo EQ 1020 OR jpjo EQ 1019):onsphere = yy GT 40 
    181191          ELSE:onsphere = 1 
    182192        ENDCASE 
    … …  
    233243          ELSE: BEGIN 
    234244; we keep its address 
    235         foraddr[n] = ind 
     245            foraddr[n] = ind 
    236246; keep the new coordinates 
    237247            forweight[n, 0] = xy[0] 
    … …  
    274284  a = reform(forweight[*, 0], 1, nawater) 
    275285  b = reform(forweight[*, 1], 1, nawater) 
    276   forweight =  -1 ; free memory 
     286  forweight =  -1               ; free memory 
    277287  newaweig = [(1-a)*(1-b), (1-b)*a, a*b, (1-a)*b] 
    278   a = -1 &  b = -1 ; free memory 
     288  a = -1 &  b = -1              ; free memory 
    279289; mask the weight to suppress the corner located on land 
    280290  newaweig = newaweig*((omsk)[newaaddr]) 

  • trunk/SRC/Interpolation/extrapolate.pro

    r262 r271  
    1616; put -1 if input data are not masked 
    1717; 
    18 ; @param nb_iteration {in}{optional}{type=integer scalar}{default=10.E20} 
     18; @param nb_iteration {in}{optional}{type=integer}{default=large enough to fill everything} 
    1919; Maximum number of iterations done in the extrapolation process. If there 
    2020; is no more masked values we exit extrapolate before reaching nb_iteration 
    21 ; (to be sure to fill everything, you can use a very large value) 
    2221; 
    2322; @keyword X_PERIODIC {type=scalar, 0 or 1}{default=0} 
    … …  
    3231; @keyword GE0 {type=scalar 0 or 1}{default=0} 
    3332; put 1 to force the extrapolated values to be larger than 0, same as using minval=0. 
     33; 
     34; @keyword  
     35; _EXTRA to be able to call extrapolate with _extra keyword 
    3436; 
    3537; @returns 
    … …  
    5456; 
    5557FUNCTION extrapolate, zinput, maskinput, nb_iteration, X_PERIODIC = x_periodic $ 
    56                       , MINVAL = minval, MAXVAL = maxval, GE0 = ge0 
     58                      , MINVAL = minval, MAXVAL = maxval, GE0 = ge0, _EXTRA = ex 
    5759; 
    5860  compile_opt idl2, strictarrsubs 
    5961; 
    6062; check the number of iteration used in the extrapolation. 
    61   IF n_elements(nb_iteration) EQ 0 THEN nb_iteration = 10.E20 
     63  szin = size(zinput) 
     64  IF szin[0] NE 2 THEN return, -1. ELSE szin = szin[1:2] 
     65  nx = szin[0] 
     66  ny = szin[1] 
     67  IF n_elements(nb_iteration) EQ 0 THEN nb_iteration = max(szin) 
    6268  IF nb_iteration EQ 0 THEN return, zinput 
    63   nx = (size(zinput))[1] 
    64   ny = (size(zinput))[2] 
    6569; take care of the boundary conditions... 
    6670; 

  • trunk/SRC/Interpolation/extrapsmooth.pro

    r262 r271  
    2727; put 1 to force the extrapolated values to be larger than 0, same as using minval=0. 
    2828; 
     29; @keyword  
     30; _EXTRA to be able to call extrapsmooth with _extra keyword 
     31; 
    2932; @returns 
    3033; the extrapolated array with no more masked values 
    … …  
    4750;- 
    4851; 
    49 FUNCTION extrapsmooth, in, mskin, X_PERIODIC = x_periodic, MINVAL = minval, MAXVAL = maxval, GE0 = ge0 
     52FUNCTION extrapsmooth, in, mskin, X_PERIODIC = x_periodic, MINVAL = minval, MAXVAL = maxval, GE0 = ge0, _EXTRA = ex 
    5053; 
    5154  compile_opt strictarr, strictarrsubs 

  • trunk/SRC/Interpolation/fromirr.pro

    r238 r271  
    4949; lonout, latout are used only to know the output domain size 
    5050; 
     51; @keyword  
     52; _EXTRA to be able to call fromirr with _extra keyword 
     53; 
    5154; @returns 
    5255; 2D array the interpolated data 
    … …  
    8790; 
    8891FUNCTION fromirr, method, datain, lonin, latin, mskin, lonout, latout, mskout $ 
    89                   , WEIG = weig, ADDR = addr 
     92                  , WEIG = weig, ADDR = addr, _EXTRA = ex 
    9093; 
    9194  compile_opt strictarr, strictarrsubs 

  • trunk/SRC/Interpolation/fromreg.pro

    r238 r271  
    5353; of the input data when performing the interpolation. 
    5454; 
     55; @keyword  
     56; _EXTRA to be able to call fromreg with _extra keyword 
     57; 
    5558; @returns 
    5659; 2D array the interpolated data 
    … …  
    9093                  , WEIG = weig, ADDR = addr $ 
    9194                  , NONORTHERNLINE = nonorthernline $ 
    92                   , NOSOUTHERNLINE = nosouthernline 
     95                  , NOSOUTHERNLINE = nosouthernline, _EXTRA = ex 
    9396; 
    9497  compile_opt idl2, strictarrsubs 

  • trunk/SRC/Interpolation/get_gridparams.pro

    r242 r271  
    22; 
    33; @file_comments 
    4 ; 1) extract from a NetCDF file the longitude, latitude, and their dimensions 
     4; Case 1: extract from a NetCDF file longitude and latitude arrays, their dimensions 
    55; and make sure it is 1D or 2D arrays 
    66; 
    7 ; or 
    8 ; 2) given longitude and latitude arrays, get their dimensions and make 
     7; Case 2: given longitude and latitude arrays, get their dimensions and make 
    98; sure they are 1D or 2D arrays 
    109; 
    … …  
    1413; @examples 
    1514; 
    16 ; 1) 
    17 ; IDL> get_gridparams, file, lonname, latname, lon, lat, jpi, jpj, n_dimensions 
    18 ; 
    19 ; or 
    20 ; 
    21 ; 2) 
     15; Case 1: 
     16; IDL> get_gridparams, file name/id, lonname, latname, lon, lat, jpi, jpj, n_dimensions 
     17; 
     18; Case 2: 
    2219; IDL> get_gridparams, lon, lat, jpi, jpj, n_dimensions 
    2320; 
    2421; @param in1 {in}{required} 
    25 ; Case 1: the name of the netcdf file 
    26 ; Case 2: 1d or 2d arrays defining longitudes and latitudes. 
    27 ; Out: the variable that will contain the longitudes 
     22; Case 1: name or the id (returned by ncdf_open) of the netcdf file 
     23; Case 2: 1D or 2D arrays defining longitudes, will be forced to have 
     24;         n_dimensions dimensions when returned 
    2825; 
    2926; @param in2 {in}{required} 
    30 ; Case 1: the name of the variable that contains the longitude in the NetCDF file 
    31 ; Case 2: 1d or 2d arrays defining longitudes and latitudes. 
    32 ;         Note that these arrays are also outputs and can therefore be modified. 
    33 ; Out: the variable that will contain the latitudes 
    34 ; 
    35 ; @param in3 {in}{required} 
    36 ; Case 1: the name of the variable that contains the latitude in the NetCDF file 
    37 ; Case 2: the number of points in the longitudinal direction. 
     27; Case 1: name of the variable containing the longitude in the NetCDF file 
     28; Case 2: 1D or 2D arrays defining latitudes, will be forced to have 
     29;         n_dimensions dimensions when returned 
     30; 
     31; @param in3  
     32; Case 1: name of the variable containing the latitude in the NetCDF file 
     33; Case 2: returned number of points in longitudinal direction. 
    3834; 
    3935; @param in4 {out} 
    40 ; Case 1: the number of points in the longitudinal direction 
    41 ; Case 2: the number of points in the latitudinal direction 
    42 ; 
    43 ; @param in5 {out} 
    44 ; Case 1: the number of points in the latitudinal direction 
    45 ; Case 2: 1 or 2 to specify if lon and lat should be 1D (jpi or jpj) 
    46 ; arrays or 2D arrays (jpi,jpj). Note that of  n_dimensions = 1, then the 
    47 ; grid must be regular (each longitude must be the same for all latitudes 
    48 ; and each latitude should be the same for all longitudes). 
     36; Case 1: returned longitudes array, with n_dimensions dimensions 
     37; Case 2: returned number of points in latitudinal direction 
     38; 
     39; @param in5 
     40; Case 1: returned latitudes array, with n_dimensions dimensions  
     41; Case 2: input scalar (1 or 2) to specify if lon and lat should be returned  
     42;         as 1D or 2D arrays. Note that if n_dimensions = 1, the grid must be 
     43;         regular (longitude and latitudes can be described as 1D arrays). 
    4944; 
    5045; @param in6 {out} 
    51 ; the variable that will contain the longitudes 
     46; Case 1: returned number of points in longitudinal direction. 
    5247; 
    5348; @param in7 {out} 
    54 ; the variable that will contain the latitudes 
    55 ; 
    56 ; @param in8 {out} 
    57 ; 1 or 2 to specify if lon and lat should be 1D (jpi or jpj) 
    58 ; 
    59 ; @keyword DOUBLE 
    60 ; use double precision to perform the computation 
     49; Case 1: returned number of points in latitudinal direction 
     50; 
     51; @param in8 {in} 
     52; Case 1: input scalar (1 or 2) to specify if lon and lat should be returned  
     53;         as 1D or 2D arrays. Note that if n_dimensions = 1, the grid must be 
     54;         regular (longitude and latitudes can be described as 1D arrays). 
     55; 
     56; @keyword DOUBLE {default=0}{type=scalar: 0 or 1} 
     57; activate to force double precision for lon/lat arrays 
    6158; 
    6259; @examples 
    … …  
    8481    8:BEGIN 
    8582; get longitude and latitude 
    86       IF file_test(in1) EQ 0 THEN BEGIN 
    87         ras = report('file ' + in1 + ' does not exist') 
    88         stop 
    89       ENDIF 
    90       cdfido = ncdf_open(in1) 
     83      IF size(in1, /type) EQ 7 THEN BEGIN 
     84        IF file_test(in1) EQ 0 THEN BEGIN 
     85          ras = report('file ' + in1 + ' does not exist') 
     86          stop 
     87        ENDIF 
     88        cdfido = ncdf_open(in1) 
     89      ENDIF ELSE cdfido = in1 
    9190      ncdf_varget, cdfido, in2, lon 
    9291      ncdf_varget, cdfido, in3, lat 
    93       ncdf_close, cdfido 
     92      IF size(in1, /type) EQ 7 THEN ncdf_close, cdfido 
    9493      n_dimensions = in8 
    9594    END 

  • trunk/SRC/ReadWrite/ncdf_getaxis.pro

    r232 r271  
    55; 
    66; @categories 
    7 ; Reading 
     7; Read NetCDF file 
    88; 
    99; @param cdfid {in}{required}{type=scalar} 
    … …  
    2929; 
    3030; @keyword XDIMNAME {default='longitude', 'lon', 'x', 'longitude*', 'lon*', 'x*', '*longitude*', '*lon*' or '*x*'}{type=scalar string} 
    31 ; A string giving the name of the x dimension 
     31; A string giving the name of the x dimension or/and a named variable 
     32; in which x dimension name is returned. 
    3233; 
    3334; @keyword YDIMNAME {default='latitude', 'lat', 'y', 'latitude*', 'lat*', 'y*', 'eta_*', '*latitude*', '*lat*', '*y*'}{type=scalar string} 
    34 ; A string giving the name of the y dimension 
     35; A string giving the name of the y dimension or/and a named variable 
     36; in which y dimension name is returned. 
    3537; 
    3638; @keyword XAXISNAME {default='x', 'longitude', 'nav_lon', 'lon', 'lon_rho' or 'NbLongitudes'}{type=scalar string} 
    3739; A string giving the name of the variable in the file 
    38 ; that contains the [xyz]axis. 
     40; that contains the x axis or/and a named variable 
     41; in which this variable name is returned. 
    3942; 
    4043; @keyword YAXISNAME {default='y', 'latitude', 'nav_lat','lat', 'lat_rho' or 'NbLatitudes'}{type=scalar string} 
    4144; A string giving the name of the variable in the file 
    42 ; that contains the [xyz]axis. 
     45; that contains the y axis or/and a named variable 
     46; in which this variable name is returned. 
    4347; 
    4448; @keyword XYINDEX {default=0}{type=scalar: 0 or 1} 
    … …  
    4751; x/yaxis = keyword_set(start1) + findgen(jpi/jpj) 
    4852; 
    49 ; @keyword XYINDEX {default=0}{type=scalar: 0 or 1} 
    50 ; 
    5153; @keyword _EXTRA 
    5254; Used to be able to call ncdf_getaxis with _extra 
    … …  
    6062;- 
    6163; 
    62 PRO ncdf_getaxis, cdfid, dimidx, dimidy, xaxis, yaxis $ 
     64PRO ncdf_getaxis, fileid, dimidx, dimidy, xaxis, yaxis $ 
    6365                  , XAXISNAME = xaxisname, YAXISNAME = yaxisname $ 
    6466                  , XDIMNAME = xdimname, YDIMNAME = ydimname $ 
    … …  
    6769  compile_opt idl2, strictarrsubs 
    6870; 
    69  
     71; should we open the file? 
     72  IF size(fileid, /type) EQ 7 THEN cdfid = ncdf_open(fileid) ELSE cdfid = fileid 
    7073; what is inside the file 
    7174  inside = ncdf_inquire(cdfid) 
    … …  
    97100                    , '''longitude'', ''nav_lon'', ''lon'', ''lon_rho'', ''nblongitudes''' $ 
    98101                    , ' we use a fake xaxis based on x dimension size (or use XAXISNAME keyword)'], /simple) 
     102    xaxisname = 'Not Found' 
    99103; try to get the dimension corresponding to x 
    100104; roms file? 
    … …  
    125129        ENDELSE 
    126130      ENDCASE 
    127     ENDIF 
    128     romsgrid = 0b 
     131      romsgrid = 0b 
     132    ENDIF ELSE romsgrid = 1b 
    129133  ENDIF ELSE BEGIN 
    130134    romsgrid = strmid(namevar[xvarid], 0, 4) EQ 'lon_' 
    131135    xinq = ncdf_varinq(cdfid, xvarid) 
     136    xaxisname = xinq.name 
    132137    dimidx = xinq.dim[0] 
    133138    IF xinq.ndims GE 2 THEN ncdf_diminq, cdfid, xinq.dim[1], blabla, jpjfromx 
    134139  ENDELSE 
     140  IF arg_present(xdimname) THEN ncdf_diminq, cdfid, dimidx,  xdimname, jpifromx 
    135141; 
    136142  IF arg_present(xaxis) THEN BEGIN 
    … …  
    162168                    , '''latitude'', ''nav_lat'', ''lat'', ''lat_rho'', ''nblatitudes''' $ 
    163169                    , ' we use a fake yaxis based on y dimension size (or use YAXISNAME keyword)'], /simple) 
     170    yaxisname = 'Not Found' 
    164171; try to get the dimension corresponding to y 
    165172; roms file? 
    … …  
    193200  ENDIF ELSE BEGIN 
    194201    yinq = ncdf_varinq(cdfid, yvarid) 
     202    yaxisname = yinq.name 
    195203    IF yinq.ndims GE 2 THEN BEGIN 
    196204      ncdf_diminq, cdfid, yinq.dim[0], blabla, jpifromy 
    … …  
    198206    ENDIF ELSE dimidy = yinq.dim[0] 
    199207  ENDELSE 
     208  IF arg_present(ydimname) THEN ncdf_diminq, cdfid, dimidy,  ydimname, jpjfromy 
    200209; 
    201210  IF arg_present(yaxis) THEN BEGIN 
    … …  
    224233  ENDIF 
    225234 
     235  IF size(fileid, /type) EQ 7 THEN ncdf_close, cdfid 
     236 
    226237  return 
    227238END 

  • 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.