npp/01.read_ascii_933.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Thu, 26 Mar 2009 14:02:21 -0400
changeset 1 4be95183fbcd
permissions -rw-r--r--
Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB.
forrest@0
     1
;----------------------------------------------------------------------
forrest@0
     2
; This example reads an ASCII file that is formatted a specific way, and
forrest@0
     3
; writes out the results to a netCDF file.
forrest@0
     4
;
forrest@0
     5
; The first line in the ASCII file must be a header, with each field
forrest@0
     6
; separated by a single character delimiter (like a ","). The rest of
forrest@0
     7
; the file must be such that each row contains all fields, each
forrest@0
     8
; separated by the designated delimiter.
forrest@0
     9
;
forrest@0
    10
; The fields can be integer, float, double, character, or string.
forrest@0
    11
; String fields cannot be written to a netCDF file. They have to
forrest@0
    12
; be read in as character arrays and written out that way.
forrest@0
    13
;----------------------------------------------------------------------
forrest@0
    14
forrest@0
    15
;----------------------------------------------------------------------
forrest@0
    16
; This function returns the index locations of the given delimiter
forrest@0
    17
; in a row or several rows of strings.
forrest@0
    18
;----------------------------------------------------------------------
forrest@0
    19
function delim_indices(strings,nfields,delimiter)
forrest@0
    20
local cstrings, cdelim
forrest@0
    21
begin
forrest@0
    22
  nrows = dimsizes(strings)
forrest@0
    23
;
forrest@0
    24
; Handle special case if we only have one string. Make sure it
forrest@0
    25
; is put into a 2D array.
forrest@0
    26
;
forrest@0
    27
  if(nrows.eq.1) then
forrest@0
    28
    cstrings = new((/1,strlen(strings)+1/),character)
forrest@0
    29
  end if
forrest@0
    30
forrest@0
    31
  cstrings = stringtochar(strings)        ; Convert to characters.
forrest@0
    32
  cdelim   = stringtochar(delimiter)      ; Convert delimiter to character.
forrest@0
    33
;
forrest@0
    34
; Som error checking here. Make sure delimiter is one character.
forrest@0
    35
;
forrest@0
    36
  nc   = dimsizes(cdelim)
forrest@0
    37
  rank = dimsizes(nc)
forrest@0
    38
  if(rank.ne.1.or.(rank.eq.1.and.nc.ne.2)) then
forrest@0
    39
    print("delim_indices: fatal: the delimiter you've selected")
forrest@0
    40
    print("must be a single character. Can't continue.")
forrest@0
    41
    exit
forrest@0
    42
  end if
forrest@0
    43
forrest@0
    44
;
forrest@0
    45
; Create array to hold indices of delimiter locations, and then loop
forrest@0
    46
; through each row and find all the delimiters. Make sure each row has
forrest@0
    47
; the correct number of delimiters.
forrest@0
    48
;
forrest@0
    49
  ndelims  = nfields-1
forrest@0
    50
  cindices = new((/nrows,ndelims/),integer)
forrest@0
    51
  do i = 0, nrows-1
forrest@0
    52
    ii = ind(cstrings(i,:).eq.cdelim(0))
forrest@0
    53
;
forrest@0
    54
; Make sure there were delimiters on this row. If not, we just quit.
forrest@0
    55
; This could probably be modified to do this more gracefully.
forrest@0
    56
;
forrest@0
    57
    if(any(ismissing(ii))) then
forrest@0
    58
      print("delim_indices: fatal: I didn't find any delimiters")
forrest@0
    59
      print("('" + delimiter + "') on row " + i + ". Can't continue.")
forrest@0
    60
      exit
forrest@0
    61
    end if
forrest@0
    62
    if(dimsizes(ii).ne.ndelims) then
forrest@0
    63
      print("delim_indices: fatal: I expected to find " + ndelims)
forrest@0
    64
      print("delimiters on row " + i + ". Instead, I found " + dimsizes(ii) + ".")
forrest@0
    65
      print("Can't continue.")
forrest@0
    66
      exit
forrest@0
    67
    end if
forrest@0
    68
forrest@0
    69
    cindices(i,:) = ii
forrest@0
    70
forrest@0
    71
    delete(ii)            ; For next time through loop
forrest@0
    72
  end do
forrest@0
    73
forrest@0
    74
  return(cindices)
forrest@0
    75
end
forrest@0
    76
forrest@0
    77
;----------------------------------------------------------------------
forrest@0
    78
; This function reads in a particular field from a string array,
forrest@0
    79
; given the field number to read (fields start at #1 and go to #nfield),
forrest@0
    80
; and the indices of the delimiters.
forrest@0
    81
;
forrest@0
    82
; It returns either an integer, float, double, character, or a string,
forrest@0
    83
; depending on the input flag "return_type".
forrest@0
    84
;----------------------------------------------------------------------
forrest@0
    85
function read_field(strings,ifield,indices,return_type)
forrest@0
    86
local nstring, cstrings, nf, tmp_str
forrest@0
    87
begin
forrest@0
    88
  nrows = dimsizes(strings)
forrest@0
    89
;
forrest@0
    90
; Handle special case if we only have one string. Make sure it
forrest@0
    91
; is put into a 2D array.
forrest@0
    92
;
forrest@0
    93
  if(nrows.eq.1) then
forrest@0
    94
    cstrings = new((/1,strlen(strings)+1/),character)
forrest@0
    95
  end if
forrest@0
    96
forrest@0
    97
  cstrings = stringtochar(strings)
forrest@0
    98
  nf       = dimsizes(indices(0,:))+1     ; indices is nrows x (nfields-1)
forrest@0
    99
forrest@0
   100
;
forrest@0
   101
; Error checking. Make sure user has entered a valid field.
forrest@0
   102
;
forrest@0
   103
  if(ifield.le.0.or.ifield.gt.nf) then
forrest@0
   104
    print("read_field: fatal: you've selected a field that is")
forrest@0
   105
    print("out-of-range of the number of fields that you have (" + nf + ").")
forrest@0
   106
    exit
forrest@0
   107
  end if
forrest@0
   108
forrest@0
   109
;
forrest@0
   110
; Set up array to return. For string, int, float, or double arrays,
forrest@0
   111
; we don't have to do anything special. For character arrays,
forrest@0
   112
; however, we do.
forrest@0
   113
;
forrest@0
   114
  if(return_type.ne."character") then
forrest@0
   115
    return_array = new(nrows,return_type)
forrest@0
   116
  else
forrest@0
   117
;
forrest@0
   118
; We don't know what the biggest character array is at this point, so
forrest@0
   119
; make it bigger than necessary, and then resize later as necessary.
forrest@0
   120
;
forrest@0
   121
    tmp_return_array = new((/nrows,dimsizes(cstrings(0,:))/),"character")
forrest@0
   122
forrest@0
   123
    max_len = 0     ; Use to keep track of max lengths of strings.
forrest@0
   124
  end if
forrest@0
   125
forrest@0
   126
  do i = 0,nrows-1
forrest@0
   127
;
forrest@0
   128
; Special case of first field in row.
forrest@0
   129
;
forrest@0
   130
    if(ifield.eq.1) then
forrest@0
   131
      ibeg = 0
forrest@0
   132
      iend = indices(i,ifield-1)-1
forrest@0
   133
    else
forrest@0
   134
;
forrest@0
   135
; Special case of first field in row.
forrest@0
   136
;
forrest@0
   137
      if(ifield.eq.nf) then
forrest@0
   138
        ibeg = indices(i,ifield-2)+1
forrest@0
   139
        iend = dimsizes(cstrings(i,:))-1
forrest@0
   140
;
forrest@0
   141
; Any field between first and last field.
forrest@0
   142
;
forrest@0
   143
      else
forrest@0
   144
        ibeg = indices(i,ifield-2)+1
forrest@0
   145
        iend = indices(i,ifield-1)-1
forrest@0
   146
      end if  
forrest@0
   147
    end if
forrest@0
   148
;
forrest@0
   149
; Here's the code that pulls off the correct string, and converts it
forrest@0
   150
; to float if desired.
forrest@0
   151
;
forrest@0
   152
    if(return_type.eq."integer") then
forrest@0
   153
      return_array(i) = stringtointeger(chartostring(cstrings(i,ibeg:iend)))
forrest@0
   154
    end if
forrest@0
   155
    if(return_type.eq."float") then
forrest@0
   156
      return_array(i) = stringtofloat(chartostring(cstrings(i,ibeg:iend)))
forrest@0
   157
    end if
forrest@0
   158
    if(return_type.eq."double") then
forrest@0
   159
      return_array(i) = stringtodouble(chartostring(cstrings(i,ibeg:iend)))
forrest@0
   160
    end if
forrest@0
   161
    if(return_type.eq."string") then
forrest@0
   162
      return_array(i) = chartostring(cstrings(i,ibeg:iend))
forrest@0
   163
    end if
forrest@0
   164
    if(return_type.eq."character") then
forrest@0
   165
      if( (iend-ibeg+1) .gt. max_len) then
forrest@0
   166
        max_len = iend-ibeg+1
forrest@0
   167
      end if
forrest@0
   168
      tmp_return_array(i,0:iend-ibeg) = cstrings(i,ibeg:iend)
forrest@0
   169
    end if
forrest@0
   170
  end do
forrest@0
   171
forrest@0
   172
  if(return_type.eq."character") then
forrest@0
   173
    return_array = new((/nrows,max_len/),"character")
forrest@0
   174
    return_array = tmp_return_array(:,0:max_len-1)
forrest@0
   175
  end if
forrest@0
   176
forrest@0
   177
  return(return_array)
forrest@0
   178
end
forrest@0
   179
forrest@0
   180
forrest@0
   181
;----------------------------------------------------------------------
forrest@0
   182
; This function reads in string fields only to get the maximum string
forrest@0
   183
; length.
forrest@0
   184
;----------------------------------------------------------------------
forrest@0
   185
function get_maxlen(strings,ifield,indices)
forrest@0
   186
local nstring, cstrings, nf, tmp_str
forrest@0
   187
begin
forrest@0
   188
  nrows = dimsizes(strings)
forrest@0
   189
;
forrest@0
   190
; Handle special case if we only have one string. Make sure it
forrest@0
   191
; is put into a 2D array.
forrest@0
   192
;
forrest@0
   193
  if(nrows.eq.1) then
forrest@0
   194
    cstrings = new((/1,strlen(strings)+1/),character)
forrest@0
   195
  end if
forrest@0
   196
forrest@0
   197
  cstrings = stringtochar(strings)
forrest@0
   198
  nf       = dimsizes(indices(0,:))+1     ; indices is nrows x (nfields-1)
forrest@0
   199
forrest@0
   200
;
forrest@0
   201
; Error checking. Make sure user has entered a valid field.
forrest@0
   202
;
forrest@0
   203
  if(ifield.le.0.or.ifield.gt.nf) then
forrest@0
   204
    print("read_field: fatal: you've selected a field that is")
forrest@0
   205
    print("out-of-range of the number of fields that you have (" + nf + ").")
forrest@0
   206
    exit
forrest@0
   207
  end if
forrest@0
   208
;
forrest@0
   209
; We don't know what the biggest character array is at this point, so
forrest@0
   210
; make it bigger than necessary, and then resize later as necessary.
forrest@0
   211
;
forrest@0
   212
  tmp_return_array = new((/nrows,dimsizes(cstrings(0,:))/),"character")
forrest@0
   213
forrest@0
   214
  max_len = 0     ; Use to keep track of max lengths of strings.
forrest@0
   215
forrest@0
   216
  do i = 0,nrows-1
forrest@0
   217
;
forrest@0
   218
; Special case of first field in row.
forrest@0
   219
;
forrest@0
   220
    if(ifield.eq.1) then
forrest@0
   221
      ibeg = 0
forrest@0
   222
      iend = indices(i,ifield-1)-1
forrest@0
   223
    else
forrest@0
   224
;
forrest@0
   225
; Special case of first field in row.
forrest@0
   226
;
forrest@0
   227
      if(ifield.eq.nf) then
forrest@0
   228
        ibeg = indices(i,ifield-2)+1
forrest@0
   229
        iend = dimsizes(cstrings(i,:))-1
forrest@0
   230
;
forrest@0
   231
; Any field between first and last field.
forrest@0
   232
;
forrest@0
   233
      else
forrest@0
   234
        ibeg = indices(i,ifield-2)+1
forrest@0
   235
        iend = indices(i,ifield-1)-1
forrest@0
   236
      end if  
forrest@0
   237
    end if
forrest@0
   238
    if( (iend-ibeg+1) .gt. max_len) then
forrest@0
   239
      max_len = iend-ibeg+1
forrest@0
   240
    end if
forrest@0
   241
  end do
forrest@0
   242
forrest@0
   243
  return(max_len)
forrest@0
   244
end
forrest@0
   245
forrest@0
   246
;----------------------------------------------------------------------
forrest@0
   247
; Main code.
forrest@0
   248
;----------------------------------------------------------------------
forrest@0
   249
begin
forrest@0
   250
;
forrest@0
   251
; Set up defaults here.  We are hard-coding the field types here.
forrest@0
   252
; You can set up this script to try to determine the field types 
forrest@0
   253
; automatically, but this is a bit tedious. Maybe later.
forrest@0
   254
;
forrest@0
   255
  filename  = "data.933"                ; ASCII" file to read.
forrest@0
   256
  cdf_file  = filename + ".nc"          ; netCDF file to write.
forrest@0
   257
  nfields   = 45                        ; # of fields
forrest@0
   258
  delimiter = ","                       ; field delimiter
forrest@0
   259
;
forrest@0
   260
; In this case, fields #6-#8 are strings, fields #2, #3, and #11
forrest@0
   261
; are float, and the rest of the fields are integers.
forrest@0
   262
;
forrest@0
   263
  var_types      = new(nfields,string)
forrest@0
   264
  var_strlens    = new(nfields,integer)   ; var to hold strlens, just in case.
forrest@0
   265
forrest@0
   266
  var_types      = "integer"       ; Most are ints.
forrest@0
   267
  var_types(5:7) = "character"     ; Corresponds to fields 6-8.
forrest@0
   268
  var_types(1:2) = "float"
forrest@0
   269
  var_types(10)  = "float"
forrest@0
   270
forrest@0
   271
  if(isfilepresent(cdf_file))
forrest@0
   272
    print("Warning: '" + cdf_file + "' exists. Will remove it.")
forrest@0
   273
    system("/bin/rm " + cdf_file)
forrest@0
   274
  end if
forrest@0
   275
forrest@0
   276
;
forrest@0
   277
; Read in data as strings. This will create a string array that has the
forrest@0
   278
; same number of strings as there are rows in the file. We will then need
forrest@0
   279
; to parse each string later.
forrest@0
   280
;
forrest@0
   281
  read_data = asciiread(filename,-1,"string")
forrest@0
   282
  header    = read_data(0)        ; Header. Use for variable names.
forrest@0
   283
  data      = read_data(1:)       ; Get rid of first line which is a header.
forrest@0
   284
  nrows     = dimsizes(data)      ; Number of rows.
forrest@0
   285
forrest@0
   286
;
forrest@0
   287
; Read in locations of delimiters in each string row.
forrest@0
   288
;
forrest@0
   289
  hindices = delim_indices(header,nfields,delimiter)   ; header row
forrest@0
   290
  dindices = delim_indices(data,nfields,delimiter)     ; rest of file
forrest@0
   291
forrest@0
   292
;
forrest@0
   293
; Read in the field names which will become variable names on
forrest@0
   294
; the netCDF file.
forrest@0
   295
;
forrest@0
   296
  var_names = new(nfields,string)
forrest@0
   297
forrest@0
   298
  do i=0,nfields-1
forrest@0
   299
    var_names(i) = read_field(header,i+1,hindices,"string")
forrest@0
   300
  end do
forrest@0
   301
forrest@0
   302
;
forrest@0
   303
; Write out this netCDF file efficiently so it will be faster.
forrest@0
   304
; Try to predefine everything before you write to it.
forrest@0
   305
; 
forrest@0
   306
  f = addfile(cdf_file,"c")
forrest@0
   307
  setfileoption(f,"DefineMode",True)       ; Enter predefine phase.
forrest@0
   308
forrest@0
   309
;
forrest@0
   310
; Write global attributes to file. It's okay to do this before 
forrest@0
   311
; predefining the file's variables. We are still in "define" mode.
forrest@0
   312
;
forrest@0
   313
  fAtt               = True
forrest@0
   314
  fAtt@description   = "Data read in from " + filename + " ASCII file."
forrest@0
   315
  fAtt@creation_date = systemfunc ("date")        
forrest@0
   316
  fileattdef( f, fAtt )        
forrest@0
   317
forrest@0
   318
;
forrest@0
   319
; Write dimension names to file. If there are no character variables,
forrest@0
   320
; then there's only one dimension name ("nvalues").
forrest@0
   321
;
forrest@0
   322
; Otherwise, we need to write a dimension name for every character
forrest@0
   323
; variable, which will indicate the maximum string length for that
forrest@0
   324
; variable.
forrest@0
   325
;
forrest@0
   326
  indc = ind(var_types.eq."character")
forrest@0
   327
  if(.not.any(ismissing(indc))) then
forrest@0
   328
;
forrest@0
   329
; We have to treat the character arrays special here. We need to
forrest@0
   330
; know their sizes so we can write the maximum size of each char
forrest@0
   331
; array to the netCDF file as a dimension name. This means we
forrest@0
   332
; need to read in the character variables once to get the string
forrest@0
   333
; lengths, then we'll read them again later to get the actual values.
forrest@0
   334
;
forrest@0
   335
    do i=0,dimsizes(indc)-1
forrest@0
   336
      var_strlens(indc(i)) = get_maxlen(data,indc(i)+1,dindices)
forrest@0
   337
    end do
forrest@0
   338
forrest@0
   339
    ndims    = dimsizes(indc) + 1
forrest@0
   340
    dimNames = new(ndims,string)
forrest@0
   341
    dimSizes = new(ndims,integer)
forrest@0
   342
    dimUnlim = new(ndims,logical)
forrest@0
   343
forrest@0
   344
    dimUnlim            = False
forrest@0
   345
    dimUnlim(0)         = True
forrest@0
   346
    dimNames(0)         = "nvalues"
forrest@0
   347
    dimNames(1:ndims-1) = var_names(indc) + "_StrLen"
forrest@0
   348
    dimSizes(0)         = -1
forrest@0
   349
    dimSizes(1:ndims-1) = var_strlens(indc)
forrest@0
   350
    filedimdef(f,dimNames,dimSizes,dimUnlim)
forrest@0
   351
  else
forrest@0
   352
;
forrest@0
   353
; No character variables, so just write the one dimension name.
forrest@0
   354
;
forrest@0
   355
    filedimdef(f,"nvalues",-1,True)
forrest@0
   356
  end if
forrest@0
   357
forrest@0
   358
;
forrest@0
   359
; Define each variable on the file.
forrest@0
   360
;
forrest@0
   361
; Don't deal with variables that are of type string.
forrest@0
   362
;
forrest@0
   363
  do i=0,nfields-1
forrest@0
   364
    if(var_types(i).ne."string") then
forrest@0
   365
      if(var_types(i).ne."character") then
forrest@0
   366
        filevardef(f, var_names(i), var_types(i), "nvalues")
forrest@0
   367
      else
forrest@0
   368
        filevardef(f, var_names(i), var_types(i),  \
forrest@0
   369
                   (/"nvalues",var_names(i)+"_StrLen"/))
forrest@0
   370
      end if
forrest@0
   371
    end if
forrest@0
   372
  end do
forrest@0
   373
forrest@0
   374
;
forrest@0
   375
; Loop through each field, read the values for that field, print 
forrest@0
   376
; information about the variable, and then write it to the netCDF
forrest@0
   377
; file.
forrest@0
   378
;
forrest@0
   379
  do i=0,nfields-1
forrest@0
   380
    ifield = i+1           ; Fields start at #1, not #0.
forrest@0
   381
;
forrest@0
   382
; Note: you can't write strings to a netCDF file, so these have
forrest@0
   383
; to be written out as character arrays.
forrest@0
   384
;
forrest@0
   385
    tmp_data = read_field(data,ifield,dindices,var_types(i))
forrest@0
   386
;
forrest@0
   387
; Print some info about the variable.
forrest@0
   388
;
forrest@0
   389
    print("")
forrest@0
   390
    print("Writing variable '" + var_names(i) + "' (field #" + ifield + ").")
forrest@0
   391
    print("Type is " + var_types(i) + ".")
forrest@0
   392
    if(var_types(i).ne."string".and.var_types(i).ne."character") then
forrest@0
   393
      print("min/max = " + min(tmp_data) + "/" + max(tmp_data))
forrest@0
   394
    end if
forrest@0
   395
forrest@0
   396
    if(any(ismissing(tmp_data))) then
forrest@0
   397
      print("This variable does contain missing values.")
forrest@0
   398
    else
forrest@0
   399
      print("This variable doesn't contain missing values.")
forrest@0
   400
    end if
forrest@0
   401
forrest@0
   402
    f->$var_names(i)$ = tmp_data       ; Write to netCDF file.
forrest@0
   403
forrest@0
   404
    delete(tmp_data)                   ; Delete for next round.
forrest@0
   405
  end do
forrest@0
   406
end