windIO.py

        # calculate the statistics values:
        stats['max'] = sig[i0:i1,:].max(axis=0)
        stats['min'] = sig[i0:i1,:].min(axis=0)
        stats['mean'] = sig[i0:i1,:].mean(axis=0)
        stats['std'] = sig[i0:i1,:].std(axis=0)
        stats['range'] = stats['max'] - stats['min']
        stats['absmax'] = np.absolute(sig[i0:i1,:]).max(axis=0)
        stats['rms'] = np.sqrt(np.mean(sig[i0:i1,:]*sig[i0:i1,:], axis=0))
        stats['int'] = integrate.trapz(sig[i0:i1,:], x=sig[i0:i1,0], axis=0)
        return stats

    # TODO: general signal method, this is not HAWC2 specific, move out
    def calc_fatigue(self, signal, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1):
        """
        Parameters
        ----------

        signal: 1D array
            One dimentional array containing the signal.
        no_bins: int
            Number of bins for the binning of the amplitudes.
        m: list
            Values of the slope of the SN curve.
        neq: int
            Number of equivalent cycles

        Returns
        -------
        eq: list
            Damage equivalent loads for each m value.
        """

        try:
            sig_rf = rainflow_astm(signal)
        except (TypeError) as e:
            print(e)
            return []

        if len(sig_rf) < 1 and not sig_rf:
            return []

        hist_data, x, bin_avg =  rfc_hist(sig_rf, no_bins)
        m = np.atleast_1d(m)

        eq = []
        for i in range(len(m)):
            eq.append(np.power(np.sum(0.5 * hist_data *\
                                    np.power(bin_avg, m[i])) / neq, 1. / m[i]))
        return eq

    # TODO: general signal method, this is not HAWC2 specific, move out
    def cycle_matrix(self, signal, no_bins=46, m=[3, 4, 6, 8, 10, 12]):

#        import fatigue_tools.fatigue as ft
#        cycles, ampl_bin_mean, ampl_bin_edges, mean_bin_mean, mean_edges \
#            = ft.cycle_matrix(signal, ampl_bins=no_bins, mean_bins=1,
#                              rainflow_func=ft.rainflow_windap)
#        # in this case eq = sum( n_i*S_i^m )
#        return [np.sum(cycles * ampl_bin_mean ** _m) for _m in m]

        try:
            sig_rf = rainflow_astm(signal)
        except:
            return []

        if len(sig_rf) < 1 and not sig_rf:
            return []

        hist_data, x, bin_avg =  rfc_hist(sig_rf, no_bins)
        m = np.atleast_1d(m)
        return [np.sum(0.5 * hist_data * bin_avg ** _m) for _m in m]

    def blade_deflection(self):
        """
        """

        # select all the y deflection channels
        db = misc.DictDB(self.ch_dict)

        db.search({'sensortype' : 'state pos', 'component' : 'z'})
        # sort the keys and save the mean values to an array/list
        chiz, zvals = [], []
        for key in sorted(db.dict_sel.keys()):
            zvals.append(-self.sig[:,db.dict_sel[key]['chi']].mean())
            chiz.append(db.dict_sel[key]['chi'])

        db.search({'sensortype' : 'state pos', 'component' : 'y'})
        # sort the keys and save the mean values to an array/list
        chiy, yvals = [], []
        for key in sorted(db.dict_sel.keys()):
            yvals.append(self.sig[:,db.dict_sel[key]['chi']].mean())
            chiy.append(db.dict_sel[key]['chi'])

        return np.array(zvals), np.array(yvals)

    def save_csv(self, fname, fmt='%.18e', delimiter=','):
        """
        Save to csv and use the unified channel names as columns
        """
        map_sorting = {}
        # first, sort on channel index
        for ch_key, ch in self.ch_dict.items():
            map_sorting[ch['chi']] = ch_key

        header = []
        # not all channels might be present...iterate again over map_sorting
        for chi in map_sorting:
            try:
                sensortag = self.ch_dict[map_sorting[chi]]['sensortag']
                header.append(map_sorting[chi] + ' // ' + sensortag)
            except:
                header.append(map_sorting[chi])

        # and save
        print('saving...', end='')
        np.savetxt(fname, self.sig[:,list(map_sorting.keys())], fmt=fmt,
                   delimiter=delimiter, header=delimiter.join(header))
        print(fname)

    def save_df(self, fname):
        """
        Save the HAWC2 data and sel file in a DataFrame that contains all the
        data, and all the channel information (the one from the sel file
        and the parsed from this function)
        """

        self.sig
        self.ch_details
        self.ch_dict


def ReadOutputAtTime(fname):
    """Distributed blade loading as generated by the HAWC2 output_at_time
    command.
    """
    # because the formatting is really weird, we need to sanatize it a bit
    with opent(fname, 'r') as f:
        # read the header from line 3
        f.readline()
        f.readline()
        header = f.readline().replace('\r', '').replace('\n', '')
        cols = [k.strip().replace(' ', '_') for k in header.split('#')[1:]]

#    data = pd.read_fwf(fname, skiprows=3, header=None)
#    pd.read_table(fname, sep='  ', skiprows=3)
#    data.index.names = cols

    data = np.loadtxt(fname, skiprows=3)
    return pd.DataFrame(data, columns=cols)


def ReadEigenBody(fname, debug=False):
    """
    Read HAWC2 body eigenalysis result file
    =======================================

    Parameters
    ----------

    file_path : str

    file_name : str


    Returns
    -------

    results : DataFrame
        Columns: body, Fd_hz, Fn_hz, log_decr_pct

    """

    #Body data for body number : 3 with the name :nacelle
    #Results:         fd [Hz]       fn [Hz]       log.decr [%]
    #Mode nr:  1:   1.45388E-21    1.74896E-03    6.28319E+02
    FILE = opent(fname)
    lines = FILE.readlines()
    FILE.close()

    df_dict = {'Fd_hz':[], 'Fn_hz':[], 'log_decr_pct':[], 'body':[]}
    for i, line in enumerate(lines):
        if debug: print('line nr: %5i' % i)
        # identify for which body we will read the data
        if line[:25] == 'Body data for body number':
            body = line.split(':')[2].rstrip().lstrip()
            # remove any annoying characters
            body = body.replace('\n','').replace('\r','')
            if debug: print('modes for body: %s' % body)
        # identify mode number and read the eigenfrequencies
        elif line[:8] == 'Mode nr:':
            linelist = line.replace('\n','').replace('\r','').split(':')
            #modenr = linelist[1].rstrip().lstrip()
            # text after Mode nr can be empty
            try:
                eigenmodes = linelist[2].rstrip().lstrip().split('   ')
            except IndexError:
                eigenmodes = ['0', '0', '0']

            if debug: print(eigenmodes)
            # in case we have more than 3, remove all the empty ones
            # this can happen when there are NaN values
            if not len(eigenmodes) == 3:
                eigenmodes = linelist[2].rstrip().lstrip().split(' ')
                eigmod = []
                for k in eigenmodes:
                    if len(k) > 1:
                        eigmod.append(k)
                #eigenmodes = eigmod
            else:
                eigmod = eigenmodes
            # remove any trailing spaces for each element
            for k in range(len(eigmod)):
                eigmod[k] = float(eigmod[k])#.lstrip().rstrip()

            df_dict['body'].append(body)
            df_dict['Fd_hz'].append(eigmod[0])
            df_dict['Fn_hz'].append(eigmod[1])
            df_dict['log_decr_pct'].append(eigmod[2])

    return pd.DataFrame(df_dict)


def ReadEigenStructure(file_path, file_name, debug=False, max_modes=500):
    """
    Read HAWC2 structure eigenalysis result file
    ============================================

    The file looks as follows:
    #0 Version ID : HAWC2MB 11.3
    #1 ___________________________________________________________________
    #2 Structure eigenanalysis output
    #3 ___________________________________________________________________
    #4 Time : 13:46:59
    #5 Date : 28:11.2012
    #6 ___________________________________________________________________
    #7 Results:         fd [Hz]       fn [Hz]       log.decr [%]
    #8 Mode nr:  1:   3.58673E+00    3.58688E+00    5.81231E+00
    #...
    #302  Mode nr:294:   0.00000E+00    6.72419E+09    6.28319E+02

    Parameters
    ----------

    file_path : str

    file_name : str

    debug : boolean, default=False

    max_modes : int
        Stop evaluating the result after max_modes number of modes have been
        identified

    Returns
    -------

    modes_arr : ndarray(3,n)
        An ndarray(3,n) holding Fd, Fn [Hz] and the logarithmic damping
        decrement [%] for n different structural eigenmodes

    """

    #0 Version ID : HAWC2MB 11.3
    #1 ___________________________________________________________________
    #2 Structure eigenanalysis output
    #3 ___________________________________________________________________
    #4 Time : 13:46:59
    #5 Date : 28:11.2012
    #6 ___________________________________________________________________
    #7 Results:         fd [Hz]       fn [Hz]       log.decr [%]
    #8 Mode nr:  1:   3.58673E+00    3.58688E+00    5.81231E+00
    #  Mode nr:294:   0.00000E+00    6.72419E+09    6.28319E+02

    FILE = opent(os.path.join(file_path, file_name))
    lines = FILE.readlines()
    FILE.close()

    header_lines = 8

    # we now the number of modes by having the number of lines
    nrofmodes = len(lines) - header_lines

    modes_arr = np.ndarray((3,nrofmodes))

    for i, line in enumerate(lines):
        if i > max_modes:
            # cut off the unused rest
            modes_arr = modes_arr[:,:i]
            break

        # ignore the header
        if i < header_lines:
            continue

        # split up mode nr from the rest
        parts = line.split(':')
        #modenr = int(parts[1])
        # get fd, fn and damping, but remove all empty items on the list
        modes_arr[:,i-header_lines]=misc.remove_items(parts[2].split(' '),'')

    return modes_arr


class UserWind(object):
    """
    """

    def __init__(self):
        pass

    def __call__(self, z_h, r_blade_tip, a_phi=None, shear_exp=None, nr_hor=3,
                 nr_vert=20, h_ME=500.0, fname=None, wdir=None):
        """

        Parameters
        ----------

        z_h : float
            Hub height

        r_blade_tip : float
            Blade tip radius

        a_phi : float, default=None
            :math:`a_{\\varphi} \\approx 0.5` parameter for the modified
            Ekman veer distribution. Values vary between -1.2 and 0.5.

        shear_exp : float, default=None

        nr_vert : int, default=3

        nr_hor : int, default=20

        h_ME : float, default=500
            Modified Ekman parameter. Take roughly 500 for off shore sites,
            1000 for on shore sites.

        fname : str, default=None
            When specified, the HAWC2 user defined veer input file will be
            written.

        wdir : float, default=None
            A constant veer angle, or yaw angle. Equivalent to setting the
            wind direction. Angle in degrees.

        Returns
        -------

        None

        """

        x, z = self.create_coords(z_h, r_blade_tip, nr_vert=nr_vert,
                                  nr_hor=nr_hor)
        if a_phi is not None:
            phi_rad = self.veer_ekman_mod(z, z_h, h_ME=h_ME, a_phi=a_phi)
            assert len(phi_rad) == nr_vert
        else:
            nr_vert = len(z)
            phi_rad = np.zeros((nr_vert,))
        # add any yaw error on top of
        if wdir is not None:
            # because wdir cw positive, and phi veer ccw positive
            phi_rad -= (wdir*np.pi/180.0)
        u, v, w, xx, zz = self.decompose_veer(phi_rad, x, z)
        # scale the shear on top of that
        if shear_exp is not None:
            shear = self.shear_powerlaw(zz, z_h, shear_exp)
            uu = u*shear[:,np.newaxis]
            vv = v*shear[:,np.newaxis]
            ww = w*shear[:,np.newaxis]
        # and write to a file
        if fname is not None:
            self.write_user_defined_shear(fname, uu, vv, ww, xx, zz)

    def create_coords(self, z_h, r_blade_tip, nr_vert=3, nr_hor=20):
        """
        Utility to create the coordinates of the wind field based on hub heigth
        and blade length.
        """
        # take 15% extra space after the blade tip
        z = np.linspace(0, z_h + r_blade_tip*1.15, nr_vert)
        # along the horizontal, coordinates with 0 at the rotor center
        x = np.linspace(-r_blade_tip*1.15, r_blade_tip*1.15, nr_hor)

        return x, z

    def shear_powerlaw(self, z, z_ref, a):
        profile = np.power(z/z_ref, a)
        # when a negative, make sure we return zero and not inf
        profile[np.isinf(profile)] = 0.0
        return profile

    def veer_ekman_mod(self, z, z_h, h_ME=500.0, a_phi=0.5):
        """
        Modified Ekman veer profile, as defined by Mark C. Kelly in email on
        10 October 2014 15:10 (RE: veer profile)

        .. math::
            \\varphi(z) - \\varphi(z_H) \\approx a_{\\varphi}
            e^{-\sqrt{z_H/h_{ME}}}
            \\frac{z-z_H}{\sqrt{z_H*h_{ME}}}
            \\left( 1 - \\frac{z-z_H}{2 \sqrt{z_H h_{ME}}}
            - \\frac{z-z_H}{4z_H} \\right)

        where:
        :math:`h_{ME} \\equiv \\frac{\\kappa u_*}{f}`
        and :math:`f = 2 \Omega \sin \\varphi` is the coriolis parameter,
        and :math:`\\kappa = 0.41` as the von Karman constant,
        and :math:`u_\\star = \\sqrt{\\frac{\\tau_w}{\\rho}}` friction velocity.

        For on shore, :math:`h_{ME} \\approx 1000`, for off-shore,
        :math:`h_{ME} \\approx 500`

        :math:`a_{\\varphi} \\approx 0.5`

        Parameters
        ----------

        :math:`a_{\\varphi} \\approx 0.5` parameter for the modified
            Ekman veer distribution. Values vary between -1.2 and 0.5.

        returns
        -------

        phi_rad : ndarray
            veer angle in radians

        """

        t1 = np.exp(-math.sqrt(z_h / h_ME))
        t2 = (z - z_h) / math.sqrt(z_h * h_ME)
        t3 = ( 1.0 - (z-z_h)/(2.0*math.sqrt(z_h*h_ME)) - (z-z_h)/(4.0*z_h) )

        return a_phi * t1 * t2 * t3

    def decompose_veer(self, phi_rad, x, z):
        """
        Convert a veer angle into u, v, and w components, ready for the
        HAWC2 user defined veer input file.

        Paramters
        ---------

        phi_rad : ndarray
            veer angle in radians

        method : str, default=linear
            'linear' for a linear veer, 'ekman_mod' for modified ekman method

        Returns
        -------

        u, v, w, v_coord, w_coord

        """

        nr_hor = len(x)
        nr_vert = len(z)
        assert len(phi_rad) == nr_vert

        tan_phi = np.tan(phi_rad)

        # convert veer angles to veer components in v, u. Make sure the
        # normalized wind speed remains 1!
#        u = sympy.Symbol('u')
#        v = sympy.Symbol('v')
#        tan_phi = sympy.Symbol('tan_phi')
#        eq1 = u**2.0 + v**2.0 - 1.0
#        eq2 = (tan_phi*u/v) - 1.0
#        sol = sympy.solvers.solve([eq1, eq2], [u,v], dict=True)
#        # proposed solution is:
#        u2 = np.sqrt(tan_phi**2/(tan_phi**2 + 1.0))/tan_phi
#        v2 = np.sqrt(tan_phi**2/(tan_phi**2 + 1.0))
#        # but that gives the sign switch wrong, simplify/rewrite to:
        u = np.sqrt(1.0/(tan_phi**2 + 1.0))
        v = np.sqrt(1.0/(tan_phi**2 + 1.0))*tan_phi
        # verify they are actually the same but the sign:
#        assert np.allclose(np.abs(u), np.abs(u2))
#        assert np.allclose(np.abs(v), np.abs(v2))

        u_full = u[:,np.newaxis] + np.zeros((3,))[np.newaxis,:]
        v_full = v[:,np.newaxis] + np.zeros((3,))[np.newaxis,:]
        w_full = np.zeros((nr_vert,nr_hor))

        return u_full, v_full, w_full, x, z

    def load_user_defined_veer(self, fname):
        """
        Load a user defined veer and shear file as used for HAWC2

        Returns
        -------

        u_comp, v_comp, w_comp, v_coord, w_coord, phi_deg
        """
        blok = 0
        bloks = {}
        with opent(fname) as f:
            for i, line in enumerate(f.readlines()):
                if line.strip()[0] == '#' and blok > 0:
                    bloks[blok] = i
                    blok += 1
                elif line.strip()[0] == '#':
                    continue
                elif blok == 0:
                    items = line.split(' ')
                    items = misc.remove_items(items, '')
                    nr_hor, nr_vert = int(items[0]), int(items[1])
                    blok += 1
#        nr_lines = i

        k = nr_hor + 4*nr_vert + 7
        v_comp = np.genfromtxt(fname, skiprows=3, skip_footer=i-3-3-nr_vert)
        u_comp = np.genfromtxt(fname, skiprows=3+1+nr_vert,
                               skip_footer=i-3-3-nr_vert*2)
        w_comp = np.genfromtxt(fname, skiprows=3+2+nr_vert*2,
                               skip_footer=i-3-3-nr_vert*3)
        v_coord = np.genfromtxt(fname, skiprows=3+3+nr_vert*3,
                               skip_footer=i-3-3-nr_vert*3-3)
        w_coord = np.genfromtxt(fname, skiprows=3+3+nr_vert*3+4,
                               skip_footer=i-k)
        phi_deg = np.arctan(v_comp[:,0]/u_comp[:,0])*180.0/np.pi

        return u_comp, v_comp, w_comp, v_coord, w_coord, phi_deg

    def write_user_defined_shear(self, fname, u, v, w, v_coord, w_coord,
                                 fmt_uvw='% 08.05f', fmt_coord='% 8.02f'):
        """
        """
        nr_hor = len(v_coord)
        nr_vert = len(w_coord)

        try:
            assert u.shape == v.shape
            assert u.shape == w.shape
            assert u.shape[0] == nr_vert
            assert u.shape[1] == nr_hor
        except AssertionError:
            raise ValueError('u, v, w shapes should be consistent with '
                             'nr_hor and nr_vert: u.shape: %s, nr_hor: %i, '
                             'nr_vert: %i' % (str(u.shape), nr_hor, nr_vert))

        # and create the input file
        with open(fname, 'wb') as fid:
            fid.write(b'# User defined shear file\n')
            fid.write(b'%i %i # nr_hor (v), nr_vert (w)\n' % (nr_hor, nr_vert))
            h1 = b'normalized with U_mean, nr_hor (v) rows, nr_vert (w) columns'
            fid.write(b'# v component, %s\n' % h1)
            np.savetxt(fid, v, fmt=fmt_uvw, delimiter='  ')
            fid.write(b'# u component, %s\n' % h1)
            np.savetxt(fid, u, fmt=fmt_uvw, delimiter='  ')
            fid.write(b'# w component, %s\n' % h1)
            np.savetxt(fid, w, fmt=fmt_uvw, delimiter='  ')
            h2 = b'# v coordinates (along the horizontal, nr_hor, 0 rotor center)'
            fid.write(b'%s\n' % h2)
            np.savetxt(fid, v_coord.reshape((v_coord.size,1)), fmt=fmt_coord)
            h3 = b'# w coordinates (zero is at ground level, height, nr_hor)'
            fid.write(b'%s\n' % h3)
            np.savetxt(fid, w_coord.reshape((w_coord.size,1)), fmt=fmt_coord)


class WindProfiles(object):

    def __init__(self):
        pass

    def powerlaw(self, z, z_ref, a):
        profile = np.power(z/z_ref, a)
        # when a negative, make sure we return zero and not inf
        profile[np.isinf(profile)] = 0.0
        return profile

    def veer_ekman_mod(self, z, z_h, h_ME=500.0, a_phi=0.5):
        """
        Modified Ekman veer profile, as defined by Mark C. Kelly in email on
        10 October 2014 15:10 (RE: veer profile)

        .. math::
            \\varphi(z) - \\varphi(z_H) \\approx a_{\\varphi}
            e^{-\sqrt{z_H/h_{ME}}}
            \\frac{z-z_H}{\sqrt{z_H*h_{ME}}}
            \\left( 1 - \\frac{z-z_H}{2 \sqrt{z_H h_{ME}}}
            - \\frac{z-z_H}{4z_H} \\right)

        where:
        :math:`h_{ME} \\equiv \\frac{\\kappa u_*}{f}`
        and :math:`f = 2 \Omega \sin \\varphi` is the coriolis parameter,
        and :math:`\\kappa = 0.41` as the von Karman constant,
        and :math:`u_\\star = \\sqrt{\\frac{\\tau_w}{\\rho}}` friction velocity.

        For on shore, :math:`h_{ME} \\approx 1000`, for off-shore,
        :math:`h_{ME} \\approx 500`

        :math:`a_{\\varphi} \\approx 0.5`

        Parameters
        ----------

        :math:`a_{\\varphi} \\approx 0.5` parameter for the modified
            Ekman veer distribution. Values vary between -1.2 and 0.5.

        returns
        -------

        phi_rad : ndarray
            veer angle in radians as function of height

        """

        t1 = np.exp(-math.sqrt(z_h / h_ME))
        t2 = (z - z_h) / math.sqrt(z_h * h_ME)
        t3 = ( 1.0 - (z-z_h)/(2.0*math.sqrt(z_h*h_ME)) - (z-z_h)/(4.0*z_h) )

        return a_phi * t1 * t2 * t3


class Turbulence(object):

    def __init__(self):

        pass

    def read_hawc2(self, fpath, shape):
        """
        Read the HAWC2 turbulence format
        """

        fid = open(fpath, 'rb')
        tmp = np.fromfile(fid, 'float32', shape[0]*shape[1]*shape[2])
        turb = np.reshape(tmp, shape)

        return turb

    def read_bladed(self, fpath, basename):

        fid = open(fpath + basename + '.wnd', 'rb')
        R1 = struct.unpack('h', fid.read(2))[0]
        R2 = struct.unpack('h', fid.read(2))[0]
        turb = struct.unpack('i', fid.read(4))[0]
        lat = struct.unpack('f', fid.read(4))[0]
        rough = struct.unpack('f', fid.read(4))[0]
        refh = struct.unpack('f', fid.read(4))[0]
        longti = struct.unpack('f', fid.read(4))[0]
        latti = struct.unpack('f', fid.read(4))[0]
        vertti = struct.unpack('f', fid.read(4))[0]
        dv = struct.unpack('f', fid.read(4))[0]
        dw = struct.unpack('f', fid.read(4))[0]
        du = struct.unpack('f', fid.read(4))[0]
        halfalong = struct.unpack('i', fid.read(4))[0]
        mean_ws = struct.unpack('f', fid.read(4))[0]
        VertLongComp = struct.unpack('f', fid.read(4))[0]
        LatLongComp = struct.unpack('f', fid.read(4))[0]
        LongLongComp = struct.unpack('f', fid.read(4))[0]
        Int = struct.unpack('i', fid.read(4))[0]
        seed = struct.unpack('i', fid.read(4))[0]
        VertGpNum = struct.unpack('i', fid.read(4))[0]
        LatGpNum = struct.unpack('i', fid.read(4))[0]
        VertLatComp = struct.unpack('f', fid.read(4))[0]
        LatLatComp = struct.unpack('f', fid.read(4))[0]
        LongLatComp = struct.unpack('f', fid.read(4))[0]
        VertVertComp = struct.unpack('f', fid.read(4))[0]
        LatVertComp = struct.unpack('f', fid.read(4))[0]
        LongVertComp = struct.unpack('f', fid.read(4))[0]

        points = np.fromfile(fid, 'int16', 2*halfalong*VertGpNum*LatGpNum*3)
        fid.close()
        return points

    def convert2bladed(self, fpath, basename, shape=(4096,32,32)):
        """
        Convert turbulence box to BLADED format
        """

        u = self.read_hawc2(fpath + basename + 'u.bin', shape)
        v = self.read_hawc2(fpath + basename + 'v.bin', shape)
        w = self.read_hawc2(fpath + basename + 'w.bin', shape)

        # mean velocity components at the center of the box
        v1, v2 = (shape[1]/2)-1, shape[1]/2
        w1, w2 = (shape[2]/2)-1, shape[2]/2
        ucent = (u[:,v1,w1] + u[:,v1,w2] + u[:,v2,w1] + u[:,v2,w2]) / 4.0
        vcent = (v[:,v1,w1] + v[:,v1,w2] + v[:,v2,w1] + v[:,v2,w2]) / 4.0
        wcent = (w[:,v1,w1] + w[:,v1,w2] + w[:,v2,w1] + w[:,v2,w2]) / 4.0

        # FIXME: where is this range 351:7374 coming from?? The original script
        # considered a box of lenght 8192
        umean = np.mean(ucent[351:7374])
        vmean = np.mean(vcent[351:7374])
        wmean = np.mean(wcent[351:7374])

        ustd = np.std(ucent[351:7374])
        vstd = np.std(vcent[351:7374])
        wstd = np.std(wcent[351:7374])

        # gives a slight different outcome, but that is that significant?
#        umean = np.mean(u[351:7374,15:17,15:17])
#        vmean = np.mean(v[351:7374,15:17,15:17])
#        wmean = np.mean(w[351:7374,15:17,15:17])

        # this is wrong since we want the std on the center point
#        ustd = np.std(u[351:7374,15:17,15:17])
#        vstd = np.std(v[351:7374,15:17,15:17])
#        wstd = np.std(w[351:7374,15:17,15:17])

        iu = np.zeros(shape)
        iv = np.zeros(shape)
        iw = np.zeros(shape)

        iu[:,:,:] = (u - umean)/ustd*1000.0
        iv[:,:,:] = (v - vmean)/vstd*1000.0
        iw[:,:,:] = (w - wmean)/wstd*1000.0

        # because MATLAB and Octave do a round when casting from float to int,
        # and Python does a floor, we have to round first
        np.around(iu, decimals=0, out=iu)
        np.around(iv, decimals=0, out=iv)
        np.around(iw, decimals=0, out=iw)

        return iu.astype(np.int16), iv.astype(np.int16), iw.astype(np.int16)

    def write_bladed(self, fpath, basename, shape):
        """
        Write turbulence BLADED file
        """
        # TODO: get these parameters from a HAWC2 input file
        seed = 6
        mean_ws = 11.4
        turb = 3
        R1 = -99
        R2 = 4

        du = 0.974121094
        dv = 4.6875
        dw = 4.6875

        longti = 14
        latti = 9.8
        vertti = 7

        iu, iv, iw = self.convert2bladed(fpath, basename, shape=shape)

        fid = open(fpath + basename + '.wnd', 'wb')
        fid.write(struct.pack('h', R1)) # R1
        fid.write(struct.pack('h', R2)) # R2
        fid.write(struct.pack('i', turb)) # Turb
        fid.write(struct.pack('f', 999)) # Lat
        fid.write(struct.pack('f', 999)) # rough
        fid.write(struct.pack('f', 999)) # refh
        fid.write(struct.pack('f', longti)) # LongTi
        fid.write(struct.pack('f', latti)) # LatTi
        fid.write(struct.pack('f', vertti)) # VertTi
        fid.write(struct.pack('f', dv)) # VertGpSpace
        fid.write(struct.pack('f', dw)) # LatGpSpace
        fid.write(struct.pack('f', du)) # LongGpSpace
        fid.write(struct.pack('i', shape[0]/2)) # HalfAlong
        fid.write(struct.pack('f', mean_ws)) # meanWS
        fid.write(struct.pack('f', 999.)) # VertLongComp
        fid.write(struct.pack('f', 999.)) # LatLongComp
        fid.write(struct.pack('f', 999.)) # LongLongComp
        fid.write(struct.pack('i', 999)) # Int
        fid.write(struct.pack('i', seed)) # Seed
        fid.write(struct.pack('i', shape[1])) # VertGpNum
        fid.write(struct.pack('i', shape[2])) # LatGpNum
        fid.write(struct.pack('f', 999)) # VertLatComp
        fid.write(struct.pack('f', 999)) # LatLatComp
        fid.write(struct.pack('f', 999)) # LongLatComp
        fid.write(struct.pack('f', 999)) # VertVertComp
        fid.write(struct.pack('f', 999)) # LatVertComp
        fid.write(struct.pack('f', 999)) # LongVertComp
#        fid.flush()

#        bladed2 = np.ndarray((shape[0], shape[2], shape[1], 3), dtype=np.int16)
#        for i in xrange(shape[0]):
#            for k in xrange(shape[1]):
#                for j in xrange(shape[2]):
#                    fid.write(struct.pack('i', iu[i, shape[1]-j-1, k]))
#                    fid.write(struct.pack('i', iv[i, shape[1]-j-1, k]))
#                    fid.write(struct.pack('i', iw[i, shape[1]-j-1, k]))
#                    bladed2[i,k,j,0] = iu[i, shape[1]-j-1, k]
#                    bladed2[i,k,j,1] = iv[i, shape[1]-j-1, k]
#                    bladed2[i,k,j,2] = iw[i, shape[1]-j-1, k]

        # re-arrange array for bladed format
        bladed = np.ndarray((shape[0], shape[2], shape[1], 3), dtype=np.int16)
        bladed[:,:,:,0] = iu[:,::-1,:]
        bladed[:,:,:,1] = iv[:,::-1,:]
        bladed[:,:,:,2] = iw[:,::-1,:]
        bladed_swap_view = bladed.swapaxes(1,2)
        bladed_swap_view.tofile(fid, format='%int16')

        fid.flush()
        fid.close()


class Bladed(object):

    def __init__(self):
        """
        Some BLADED results I have seen are just weird text files. Convert
        them to a more convienent format.

        path/to/file
        channel 1 description
        col a name/unit col b name/unit
        a0 b0
        a1 b1
        ...
        path/to/file
        channel 2 description
        col a name/unit col b name/unit
        ...
        """
        pass

    def infer_format(self, lines):
        """
        Figure out how many channels and time steps are included
        """
        count = 1
        for line in lines[1:]:
            if line == lines[0]:
                break
            count += 1
        iters = count - 3
        chans = len(lines) / (iters + 3)
        return int(chans), int(iters)

    def read(self, fname, chans=None, iters=None, enc='cp1252'):
        """
        Parameters
        ----------

        fname : str

        chans : int, default=None

        iters : int, default=None

        enc : str, default='cp1252'
            character encoding of the source file. Usually BLADED is used on
            windows so Western-European windows encoding is a safe bet.
        """

        with codecs.opent(fname, 'r', enc) as f:
            lines = f.readlines()
        nrl = len(lines)
        if chans is None and iters is None:
            chans, iters = self.infer_format(lines)
        if iters is not None:
            chans = int(nrl / (iters + 3))
        if chans is not None:
            iters = int((nrl / chans) - 3)
#        file_head = [ [k[:-2],0] for k in lines[0:nrl:iters+3] ]
#        chan_head = [ [k[:-2],0] for k in lines[1:nrl:iters+3] ]
#        cols_head = [ k.split('\t')[:2] for k in lines[2:nrl:iters+3] ]

        data = {}
        for k in range(chans):
            # take the column header from the 3 comment line, but
            head = lines[2 + (3 + iters)*k][:-2].split('\t')[1].encode('utf-8')
            i0 = 3 + (3 + iters)*k
            i1 = i0 + iters
            data[head] = np.array([k[:-2].split('\t')[1] for k in lines[i0:i1:1]])
            data[head] = data[head].astype(np.float64)
        time = np.array([k[:-2].split('\t')[0] for k in lines[i0:i1:1]])
        df = pd.DataFrame(data, index=time.astype(np.float64))
        df.index.name = lines[0][:-2]
        return df


if __name__ == '__main__':

    pass