1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/npp/01.read_81.ncl.x Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,271 @@
1.4 +;----------------------------------------------------------------------
1.5 +; This example reads an ASCII file that is formatted a specific way, and
1.6 +; writes out the results to a netCDF file.
1.7 +;
1.8 +; The first line in the ASCII file must be a header, with each field
1.9 +; separated by a single character delimiter (like a ","). The rest of
1.10 +; the file must be such that each row contains all fields, each
1.11 +; separated by the designated delimiter.
1.12 +;
1.13 +; The fields can be integer, float, double, or string.
1.14 +;----------------------------------------------------------------------
1.15 +
1.16 +;----------------------------------------------------------------------
1.17 +; This function returns the index locations of the given delimiter
1.18 +; in a row or several rows of strings.
1.19 +;----------------------------------------------------------------------
1.20 +function delim_indices(strings,nfields,delimiter)
1.21 +local cstrings, cdelim
1.22 +begin
1.23 + nrows = dimsizes(strings)
1.24 +;
1.25 +; Handle special case if we only have one string. Make sure it
1.26 +; is put into a 2D array.
1.27 +;
1.28 + if(nrows.eq.1) then
1.29 + cstrings = new((/1,strlen(strings)+1/),character)
1.30 + end if
1.31 +
1.32 + cstrings = stringtochar(strings) ; Convert to characters.
1.33 + cdelim = stringtochar(delimiter) ; Convert delimiter to character.
1.34 +;
1.35 +; Som error checking here. Make sure delimiter is one character.
1.36 +;
1.37 + nc = dimsizes(cdelim)
1.38 + rank = dimsizes(nc)
1.39 + if(rank.ne.1.or.(rank.eq.1.and.nc.ne.2)) then
1.40 + print("delim_indices: fatal: the delimiter you've selected")
1.41 + print("must be a single character. Can't continue.")
1.42 + exit
1.43 + end if
1.44 +
1.45 +;
1.46 +; Create array to hold indices of delimiter locations, and then loop
1.47 +; through each row and find all the delimiters. Make sure each row has
1.48 +; the correct number of delimiters.
1.49 +;
1.50 + ndelims = nfields-1
1.51 + cindices = new((/nrows,ndelims/),integer)
1.52 + do i = 0, nrows-1
1.53 + ii = ind(cstrings(i,:).eq.cdelim(0))
1.54 +;
1.55 +; Make sure there were delimiters on this row. If not, we just quit.
1.56 +; This could probably be modified to do this more gracefully.
1.57 +;
1.58 + if(any(ismissing(ii))) then
1.59 + print("delim_indices: fatal: I didn't find any delimiters")
1.60 + print("('" + delimiter + "') on row " + i + ". Can't continue.")
1.61 + exit
1.62 + end if
1.63 + if(dimsizes(ii).ne.ndelims) then
1.64 + print("delim_indices: fatal: I expected to find " + ndelims)
1.65 + print("delimiters on row " + i + ". Instead, I found " + dimsizes(ii) + ".")
1.66 + print("Can't continue.")
1.67 + exit
1.68 + end if
1.69 +
1.70 + cindices(i,:) = ii
1.71 +
1.72 + delete(ii) ; For next time through loop
1.73 + end do
1.74 +
1.75 + return(cindices)
1.76 +end
1.77 +
1.78 +;----------------------------------------------------------------------
1.79 +; This function reads in a particular field from a string array,
1.80 +; given the field number to read (fields start at 1 and go to nfield),
1.81 +; and the indices of the delimiters.
1.82 +;
1.83 +; It returns either an integer, float, double, or a string,
1.84 +; depending on the input flag "return_type".
1.85 +;----------------------------------------------------------------------
1.86 +function read_field(strings,ifield,indices,return_type)
1.87 +local nstring, cstrings, nf, tmp_str
1.88 +begin
1.89 + nrows = dimsizes(strings)
1.90 +;
1.91 +; Handle special case if we only have one string. Make sure it
1.92 +; is put into a 2D array.
1.93 +;
1.94 + if(nrows.eq.1) then
1.95 + cstrings = new((/1,strlen(strings)+1/),character)
1.96 + end if
1.97 +
1.98 + cstrings = stringtochar(strings)
1.99 + nf = dimsizes(indices(0,:))+1 ; indices is nrows x (nfields-1)
1.100 +
1.101 +;
1.102 +; Error checking. Make sure user has entered a valid field.
1.103 +;
1.104 + if(ifield.le.0.or.ifield.gt.nf) then
1.105 + print("read_field: fatal: you've selected a field that is")
1.106 + print("out-of-range of the number of fields that you have (" + nf + ").")
1.107 + exit
1.108 + end if
1.109 +
1.110 +;
1.111 +; Set up array to return.
1.112 +;
1.113 + return_array = new(nrows,return_type)
1.114 +
1.115 + do i = 0,nrows-1
1.116 +;
1.117 +; Special case of first field in row.
1.118 +;
1.119 + if(ifield.eq.1) then
1.120 + ibeg = 0
1.121 + iend = indices(i,ifield-1)-1
1.122 + else
1.123 +;
1.124 +; Special case of first field in row.
1.125 +;
1.126 + if(ifield.eq.nf) then
1.127 + ibeg = indices(i,ifield-2)+1
1.128 + iend = dimsizes(cstrings(i,:))-1
1.129 +;
1.130 +; Any field between first and last field.
1.131 +;
1.132 + else
1.133 + ibeg = indices(i,ifield-2)+1
1.134 + iend = indices(i,ifield-1)-1
1.135 + end if
1.136 + end if
1.137 +;
1.138 +; Here's the code that pulls off the correct string, and converts it
1.139 +; to float if desired.
1.140 +;
1.141 + if(return_type.eq."integer") then
1.142 + return_array(i) = stringtointeger(chartostring(cstrings(i,ibeg:iend)))
1.143 + end if
1.144 + if(return_type.eq."float") then
1.145 + return_array(i) = stringtofloat(chartostring(cstrings(i,ibeg:iend)))
1.146 + end if
1.147 + if(return_type.eq."double") then
1.148 + return_array(i) = stringtodouble(chartostring(cstrings(i,ibeg:iend)))
1.149 + end if
1.150 + if(return_type.eq."string") then
1.151 + return_array(i) = chartostring(cstrings(i,ibeg:iend))
1.152 + end if
1.153 + end do
1.154 +
1.155 + return(return_array)
1.156 +end
1.157 +
1.158 +;----------------------------------------------------------------------
1.159 +; Main code.
1.160 +;----------------------------------------------------------------------
1.161 +begin
1.162 +;
1.163 +; Set up defaults here. We have to hard-code the field types,
1.164 +; because we can't really tell otherwise.
1.165 +
1.166 + filename = "data.81" ; ASCII" file to read.
1.167 + cdf_file = filename + ".nc" ; netCDF file to write.
1.168 + nfields = 22 ; # of fields
1.169 + delimiter = "," ; field delimiter
1.170 +;
1.171 +; In this case, fields #6-#8 are strings, fields #2, #3, and #11
1.172 +; are float, and the rest of the fields are integers.
1.173 +;
1.174 + var_types = new(nfields,string)
1.175 + var_types = "integer" ; Most are ints.
1.176 + var_types(5:7) = "string" ; corresponds to fields 6-8
1.177 + var_types(1:2) = "float"
1.178 + var_types(10) = "float"
1.179 +
1.180 + if(isfilepresent(cdf_file))
1.181 + print("Warning: '" + cdf_file + "' exists. Will remove it.")
1.182 + system("/bin/rm " + cdf_file)
1.183 + end if
1.184 +
1.185 +;
1.186 +; Read in data as strings. This will create a string array that has the
1.187 +; same number of strings as there are rows in the file. We will then need
1.188 +; to parse each string later.
1.189 +;
1.190 + read_data = asciiread(filename,-1,"string")
1.191 + header = read_data(0) ; Header. Use for variable names.
1.192 + data = read_data(1:) ; Get rid of first line which is a header.
1.193 + nrows = dimsizes(data) ; Number of rows.
1.194 +
1.195 +;
1.196 +; Read in locations of delimiters in each string row.
1.197 +;
1.198 + hindices = delim_indices(header,nfields,delimiter) ; header row
1.199 + dindices = delim_indices(data,nfields,delimiter) ; rest of file
1.200 +
1.201 +;
1.202 +; Read in the field names which will become variable names on
1.203 +; the netCDF file.
1.204 +;
1.205 + var_names = new(nfields,string)
1.206 +
1.207 + do i=0,nfields-1
1.208 + var_names(i) = read_field(header,i+1,hindices,"string")
1.209 + end do
1.210 +;
1.211 +; Write out this netCDF file efficiently so it will be faster.
1.212 +; Try to predefine everything before you write to it.
1.213 +;
1.214 + f = addfile(cdf_file,"c")
1.215 + setfileoption(f,"DefineMode",True) ; Enter predefine phase.
1.216 +
1.217 +;
1.218 +; Write global attributes to file. It's okay to do this before
1.219 +; predefining the file's variables. We are still in "define" mode.
1.220 +;
1.221 + fAtt = True
1.222 + fAtt@description = "Data read in from '" + filename + "' ASCII file."
1.223 + fAtt@creation_date = systemfunc ("date")
1.224 + fileattdef( f, fAtt )
1.225 +
1.226 +;
1.227 +; Write dimension name. There's only one here.
1.228 +;
1.229 + filedimdef(f,"nvalues",-1,True)
1.230 +
1.231 +;
1.232 +; Define each variable on the file, giving it the dimension name
1.233 +; we just created above.
1.234 +;
1.235 +; Don't deal with variables that are of type string.
1.236 +;
1.237 + do i=0,nfields-1
1.238 + if(var_types(i).ne."string") then
1.239 + filevardef(f, var_names(i), var_types(i), "nvalues")
1.240 + end if
1.241 + end do
1.242 +
1.243 +;
1.244 +; Loop through each field, read the values for that field, print
1.245 +; information about the variable, and then write it to the netCDF
1.246 +; file.
1.247 +;
1.248 + do i=0,nfields-1
1.249 + ifield = i+1 ; Fields start at #1, not #0.
1.250 +;
1.251 +; You can't write strings to a netCDF file, so skip the string fields.
1.252 +;
1.253 + if(var_types(i).ne."string") then
1.254 + tmp_data = read_field(data,ifield,dindices,var_types(i))
1.255 +;
1.256 +; Print some info about the variable.
1.257 +;
1.258 + print("")
1.259 + print("Writing variable " + var_names(i) + " (field #" + ifield + ").")
1.260 + print("Type is " + var_types(i) + ", min/max = " + min(tmp_data) + \
1.261 + "/" + max(tmp_data))
1.262 +
1.263 + if(any(ismissing(tmp_data))) then
1.264 + print("This variable does contain missing values.")
1.265 + else
1.266 + print("This variable doesn't contain missing values.")
1.267 + end if
1.268 +
1.269 + f->$var_names(i)$ = tmp_data ; Write to netCDF file.
1.270 +
1.271 + delete(tmp_data) ; Delete for next round.
1.272 + end if
1.273 + end do
1.274 +end