windIO.py

        nr_hor : int
            Number of horizontal grid points

        Returns
        -------

        u : ndarray(nr_hor, nr_vert)

        v : ndarray(nr_hor, nr_vert)

        w : ndarray(nr_hor, nr_vert)

        """

        nr_vert = len(phi_rad)
        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

    def read(self, fname):
        """
        Read a user defined shear input file as used for HAWC2.

        Returns
        -------

        u_comp, v_comp, w_comp, v_coord, w_coord, phi_deg
        """
        # read the header
        with opent(fname) as f:
            for i, line in enumerate(f.readlines()):
                if line.strip()[0] != '#':
                    nr_v, nr_w = misc.remove_items(line.split('#')[0].split(), '')
                    nr_hor, nr_vert = int(nr_v), int(nr_w)
                    i_header = i
                    break

        # u,v and w components on 2D grid
        tmp = np.genfromtxt(fname, skip_header=i_header+1, comments='#',
                            max_rows=nr_vert*3)
        if not tmp.shape == (nr_vert*3, nr_hor):
            raise AssertionError('user defined shear input file inconsistent')
        v_comp = tmp[:nr_vert,:]
        u_comp = tmp[nr_vert:nr_vert*2,:]
        w_comp = tmp[nr_vert*2:nr_vert*3,:]

        # coordinates of the 2D grid
        tmp = np.genfromtxt(fname, skip_header=3*(nr_vert+1)+2,
                            max_rows=nr_hor+nr_vert)
        if not tmp.shape == (nr_vert+nr_hor,):
            raise AssertionError('user defined shear input file inconsistent')
        v_coord = tmp[:nr_hor]
        w_coord = tmp[nr_hor:]

        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(self, fid, u, v, w, v_coord, w_coord, fmt_uvw='% 08.05f',
              fmt_coord='% 8.02f'):
        """Write a user defined shear input file for HAWC2.
        """
        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))

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

        return fid


class WindProfiles(object):

    def logarithmic(z, z_ref, r_0):
        return np.log10(z/r_0)/np.log10(z_ref/r_0)

    def powerlaw(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(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
        ----------

        z : ndarray(n)
            z-coordinates (height) of the grid on which the veer angle should
            be calculated.

        z_h : float
            Hub height in meters.

        :math:`a_{\\varphi}` : default=0.5
            Parameter for the modified Ekman veer distribution. Value varies
            between -1.2 and 0.5.

        Returns
        -------

        phi_rad : ndarray
            Veer angle in radians as function of z.


        """

        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