fast_io + test

README.md

@@ -16,6 +16,9 @@ General Time Series Data Format, a binary hdf5 data format for storing time seri
 - [gtsdf](wetb/gtsdf/gtsdf.py): read/write/append gtsdf files
 - [unix_time](wetb/gtsdf/unix_time.py): convert between datetime and unix time (seconds since 1/1/1970)
 
+### [FAST](wetb/fast)
+Tools for working with NREL's FAST code (An aeroelastic computer-aided engineering (CAE) tool for horizontal axis wind turbines)
+- [fast_io](wetb/fast/fast_io.py): Read binary and ascii result files
 
 ### [functions](wetb/functions)
 Other functions

wetb/fast/ReadFASTbinary.m

+function [Channels, ChanName, ChanUnit, FileID] = ReadFASTbinary(FileName)
+%[Channels, ChannelNames, ChannelUnits] = ReadFASTbinary(FileName)
+% Author: Bonnie Jonkman, National Renewable Energy Laboratory
+% (c) 2012, National Renewable Energy Laboratory
+%
+%  Edited for FAST v7.02.00b-bjj  22-Oct-2012
+%
+% Input:
+%  FileName      - string: contains file name to open
+%
+% Output:
+%  Channels      - 2-D array: dimension 1 is time, dimension 2 is channel 
+%  ChanName      - cell array containing names of output channels
+%  ChanUnit      - cell array containing unit names of output channels
+%  FileID        - constant that determines if the time is stored in the
+%                  output, indicating possible non-constant time step
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%
+LenName = 10;  % number of characters per channel name
+LenUnit = 10;  % number of characters per unit name
+
+FileFmtID = struct( 'WithTime',   1, ...               % File identifiers used in FAST
+                    'WithoutTime',2 );
+
+fid  = fopen( FileName );
+if fid > 0
+    %----------------------------        
+    % get the header information
+    %----------------------------
+    
+    FileID       = fread( fid, 1, 'int16');             % FAST output file format, INT(2)
+
+    NumOutChans  = fread( fid, 1, 'int32');             % The number of output channels, INT(4)
+    NT           = fread( fid, 1, 'int32');             % The number of time steps, INT(4)
+
+    if FileID == FileFmtID.WithTime
+        TimeScl  = fread( fid, 1, 'float64');           % The time slopes for scaling, REAL(8)
+        TimeOff  = fread( fid, 1, 'float64');           % The time offsets for scaling, REAL(8)
+    else
+        TimeOut1 = fread( fid, 1, 'float64');           % The first time in the time series, REAL(8)
+        TimeIncr = fread( fid, 1, 'float64');           % The time increment, REAL(8)
+    end
+    
+    ColScl       = fread( fid, NumOutChans, 'float32'); % The channel slopes for scaling, REAL(4)
+    ColOff       = fread( fid, NumOutChans, 'float32'); % The channel offsets for scaling, REAL(4)
+
+    LenDesc      = fread( fid, 1,           'int32' );  % The number of characters in the description string, INT(4)
+    DescStrASCII = fread( fid, LenDesc,     'uint8' );  % DescStr converted to ASCII
+    DescStr      = char( DescStrASCII' );                     
+    
+    
+    ChanName = cell(NumOutChans+1,1);                   % initialize the ChanName cell array
+    for iChan = 1:NumOutChans+1 
+        ChanNameASCII = fread( fid, LenName, 'uint8' ); % ChanName converted to numeric ASCII
+        ChanName{iChan}= strtrim( char(ChanNameASCII') );
+    end
+    
+    ChanUnit = cell(NumOutChans+1,1);                   % initialize the ChanUnit cell array
+    for iChan = 1:NumOutChans+1
+        ChanUnitASCII = fread( fid, LenUnit, 'uint8' ); % ChanUnit converted to numeric ASCII
+        ChanUnit{iChan}= strtrim( char(ChanUnitASCII') );
+    end            
+
+    disp( ['Reading from the file ' FileName ' with heading: ' ] );
+    disp( ['   "' DescStr '".' ] ) ;
+    
+    %-------------------------        
+    % get the channel time series
+    %-------------------------
+
+    nPts        = NT*NumOutChans;           % number of data points in the file   
+    Channels    = zeros(NT,NumOutChans+1);  % output channels (including time in column 1)
+    
+    if FileID == FileFmtID.WithTime
+        [PackedTime cnt] = fread( fid, NT, 'int32' ); % read the time data
+        if ( cnt < NT ) 
+            fclose(fid);
+            error(['Could not read entire ' FileName ' file: read ' num2str( cnt ) ' of ' num2str( NT ) ' time values.']);
+        end
+    end
+        
+        
+    [PackedData cnt] = fread( fid, nPts, 'int16' ); % read the channel data
+    if ( cnt < nPts ) 
+        fclose(fid);
+        error(['Could not read entire ' FileName ' file: read ' num2str( cnt ) ' of ' num2str( nPts ) ' values.']);
+    end
+    
+    fclose(fid);
+    
+    %-------------------------
+    % Scale the packed binary to real data
+    %-------------------------
+    
+    ip = 1;
+    for it = 1:NT
+        for ic = 1:NumOutChans
+            Channels(it,ic+1) = ( PackedData(ip) - ColOff(ic) ) / ColScl(ic) ;
+            ip = ip + 1;
+        end % ic       
+    end %it
+
+    if FileID == FileFmtID.WithTime
+        Channels(:,1) = ( PackedTime - TimeOff ) ./ TimeScl;
+    else
+        Channels(:,1) = TimeOut1 + TimeIncr*(0:(NT-1))';
+    end
+    
+else
+    error(['Could not open the FAST binary file: ' FileName]) ;
+end
+
+return;
+

wetb/fast/fast_io.py

+'''
+Created on 03/09/2015
+
+@author: MMPE
+'''
+import os
+import numpy as np
+import struct
+
+def load_output(filename):
+    """Load a FAST binary or ascii output file
+
+    Parameters
+    ----------
+    filename : str
+        filename
+
+    Returns
+    -------
+    data : ndarray
+        data values
+    info : dict
+        info containing:
+            - name: filename
+            - description: description of dataset
+            - attribute_names: list of attribute names
+            - attribute_units: list of attribute units
+    """
+
+    assert os.path.isfile(filename), "File, %s, does not exists" % filename
+    with open(filename, 'r') as f:
+        try:
+            f.readline()
+        except UnicodeDecodeError:
+            return load_binary_output(filename)
+    return load_ascii_output(filename)
+
+def load_ascii_output(filename):
+    with open(filename) as f:
+        info = {}
+        info['name'] = os.path.splitext(os.path.basename(filename))[0]
+        try:
+            header = [f.readline() for _ in range(8)]
+            info['description'] = header[4].strip()
+            info['attribute_names'] = header[6].split()
+            info['attribute_units'] = [unit[1:-1] for unit in header[7].split()]  #removing "()"
+            data = np.array([line.split() for line in f.readlines()]).astype(np.float)
+
+            return data, info
+        except (ValueError, AssertionError):
+
+            raise
+
+
+def load_binary_output(filename):
+    """Ported from ReadFASTbinary.m by Mads M Pedersen, DTU Wind
+    Info about ReadFASTbinary.m:
+    % Author: Bonnie Jonkman, National Renewable Energy Laboratory
+    % (c) 2012, National Renewable Energy Laboratory
+    %
+    %  Edited for FAST v7.02.00b-bjj  22-Oct-2012
+    """
+    def fread(fid, n, type):
+        fmt, nbytes = {'uint8': ('B', 1), 'int16':('h', 2), 'int32':('i', 4), 'float32':('f', 4), 'float64':('d', 8)}[type]
+
+        return struct.unpack(fmt * n, fid.read(nbytes * n))
+
+    FileFmtID_WithTime = 1  #% File identifiers used in FAST
+    FileFmtID_WithoutTime = 2
+    LenName = 10  #;  % number of characters per channel name
+    LenUnit = 10  #;  % number of characters per unit name
+
+    with open(filename, 'rb') as fid:
+        FileID = fread(fid, 1, 'int16')  #;             % FAST output file format, INT(2)
+
+        NumOutChans = fread(fid, 1, 'int32')[0]  #;             % The number of output channels, INT(4)
+        NT = fread(fid, 1, 'int32')[0]  #;             % The number of time steps, INT(4)
+
+        if FileID == FileFmtID_WithTime:
+            TimeScl = fread(fid, 1, 'float64')  #;           % The time slopes for scaling, REAL(8)
+            TimeOff = fread(fid, 1, 'float64')  #;           % The time offsets for scaling, REAL(8)
+        else:
+            TimeOut1 = fread(fid, 1, 'float64')  #;           % The first time in the time series, REAL(8)
+            TimeIncr = fread(fid, 1, 'float64')  #;           % The time increment, REAL(8)
+
+
+
+
+        ColScl = fread(fid, NumOutChans, 'float32')  #; % The channel slopes for scaling, REAL(4)
+        ColOff = fread(fid, NumOutChans, 'float32')  #; % The channel offsets for scaling, REAL(4)
+
+        LenDesc = fread(fid, 1, 'int32')[0]  #;  % The number of characters in the description string, INT(4)
+        DescStrASCII = fread(fid, LenDesc, 'uint8')  #;  % DescStr converted to ASCII
+        DescStr = "".join(map(chr, DescStrASCII)).strip()
+
+
+
+        ChanName = []  # initialize the ChanName cell array
+        for iChan in range(NumOutChans + 1):
+            ChanNameASCII = fread(fid, LenName, 'uint8')  #; % ChanName converted to numeric ASCII
+            ChanName.append("".join(map(chr, ChanNameASCII)).strip())
+
+
+        ChanUnit = []  # initialize the ChanUnit cell array
+        for iChan in range(NumOutChans + 1):
+            ChanUnitASCII = fread(fid, LenUnit, 'uint8')  #; % ChanUnit converted to numeric ASCII
+            ChanUnit.append("".join(map(chr, ChanUnitASCII)).strip()[1:-1])
+
+
+        #    %-------------------------
+        #    % get the channel time series
+        #    %-------------------------
+
+        nPts = NT * NumOutChans  #;           % number of data points in the file
+
+
+        if FileID == FileFmtID_WithTime:
+            PackedTime = fread(fid, NT, 'int32')  #; % read the time data
+            cnt = len(PackedTime)
+            if cnt < NT:
+                raise Exception('Could not read entire %s file: read %d of %d time values' % (filename, cnt, NT))
+        PackedData = fread(fid, nPts, 'int16')  #; % read the channel data
+        cnt = len(PackedData)
+        if cnt < nPts:
+            raise Exception('Could not read entire %s file: read %d of %d values' % (filename, cnt, nPts))
+
+    #    %-------------------------
+    #    % Scale the packed binary to real data
+    #    %-------------------------
+    #
+
+
+    data = np.array(PackedData).reshape(NT, NumOutChans)
+    data = (data - ColOff) / ColScl
+
+    if FileID == FileFmtID_WithTime:
+        time = (np.array(PackedTime) - TimeOff) / TimeScl;
+    else:
+        time = TimeOut1 + TimeIncr * np.arange(NT)
+
+    data = np.concatenate([time.reshape(NT, 1), data], 1)
+
+    info = {'name': os.path.splitext(os.path.basename(filename))[0],
+            'description': DescStr,
+            'attribute_names': ChanName,
+            'attribute_units': ChanUnit}
+    return data, info
+

wetb/fast/tests/test_fast_io.py

+'''
+Created on 03/09/2015
+
+@author: MMPE
+'''
+import unittest
+
+from wetb.fast.fast_io import load_output
+
+
+class Test(unittest.TestCase):
+
+
+    def test_load_output(self):
+        data, info = load_output('testfiles/DTU10MW.out')
+        self.assertAlmostEqual(data[4, 3], 4.295E-04)
+        self.assertEqual(info['name'], "DTU10MW")
+        self.assertEqual(info['attribute_names'][1], "RotPwr")
+        self.assertEqual(info['attribute_units'][1], "kW")
+
+
+
+    def test_load_binary(self):
+        data, info = load_output('testfiles/test_binary.outb')
+        self.assertEqual(info['name'], 'test_binary')
+        self.assertEqual(info['description'], 'Modified by mwDeriveSensors on 27-Jul-2015 16:32:06')
+        self.assertEqual(info['attribute_names'][4], 'RotPwr')
+        self.assertEqual(info['attribute_units'][7], 'deg/s^2')
+        self.assertAlmostEqual(data[10, 4], 138.822277739535)
+
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testload_output']
+    unittest.main()