wetb/utils/euler.py

@@ -6,41 +6,46 @@ Created on 15/01/2014
 
 import numpy as np
 from wetb.utils.geometry import deg
+import warnings
 
 
 def Ax(angle):
     cos = np.cos(angle)
     sin = np.sin(angle)
     return np.array([[1, 0, 0],
-                   [0, cos, -sin],
-                   [0, sin, cos]])
+                     [0, cos, -sin],
+                     [0, sin, cos]])
+
 
 def Ay(angle):
     cos = np.cos(angle)
     sin = np.sin(angle)
     return np.array([[cos, 0, sin],
-                   [0, 1, 0       ],
-                   [-sin, 0, cos  ]])
+                     [0, 1, 0],
+                     [-sin, 0, cos]])
+
 
 def Az(angle):
     cos = np.cos(angle)
     sin = np.sin(angle)
     return np.array([[cos, -sin, 0],
-                   [sin, cos, 0],
-                   [0, 0, 1]])
-
-
+                     [sin, cos, 0],
+                     [0, 0, 1]])
 
 
 def euler2A(euler_param):
+    warnings.warn("deprecated, use wetb.rotation.quaternion2matrix instead", DeprecationWarning)
     assert len(euler_param) == 4
-    e = euler_param
-    return np.array([[e[0] ** 2 + e[1] ** 2 - e[2] ** 2 - e[3] ** 2, 2 * (e[1] * e[2] + e[0] * e[3]) , 2 * (e[1] * e[3] - e[0] * e[2]) ],
-                         [2 * (e[1] * e[2] - e[0] * e[3]), e[0] ** 2 - e[1] ** 2 + e[2] ** 2 - e[3] ** 2, 2 * (e[2] * e[3] + e[0] * e[1]) ],
-                         [2 * (e[1] * e[3] + e[0] * e[2]), 2 * (e[2] * e[3] - e[0] * e[1]), e[0] ** 2 - e[1] ** 2 - e[2] ** 2 + e[3] ** 2]]).T
+    e = euler_param / np.sqrt(np.sum(euler_param**2))
+    e2 = e**2
+    return np.array([[e2[0] + e2[1] - e2[2] - e2[3], 2 * (e[1] * e[2] + e[0] * e[3]), 2 * (e[1] * e[3] - e[0] * e[2])],
+                     [2 * (e[1] * e[2] - e[0] * e[3]), e2[0] - e2[1] + e2[2] - e2[3], 2 * (e[2] * e[3] + e[0] * e[1])],
+                     [2 * (e[1] * e[3] + e[0] * e[2]), 2 * (e[2] * e[3] - e[0] * e[1]), e2[0] - e2[1] - e2[2] + e2[3]]]).T
+
 
 def A2euler(A):
     # method from http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm
+    warnings.warn("deprecated, use wetb.rotation.matrix2quaternion instead", DeprecationWarning)
     sqrt = np.sqrt
     (m00, m01, m02), (m10, m11, m12), (m20, m21, m22) = A
     tr = m00 + m11 + m22
@@ -69,11 +74,11 @@ def A2euler(A):
         qx = (m02 + m20) / S
         qy = (m12 + m21) / S
         qz = 0.25 * S
-    return np.array([qw, qx, qy, qz])
-
+    e = np.array([qw, qx, qy, qz])
+    return e / np.sqrt(np.sum(e**2))
 
 
-#def A2xyz(A):
+# def A2xyz(A):
 #    if abs(A[2, 0]) != 1:
 #        y = -np.arcsin(A[2, 0])
 #        x = np.arctan2(A[2, 1] / np.cos(y), A[2, 2] / np.cos(y))
@@ -88,19 +93,19 @@ def A2euler(A):
 #            x = -z + np.arctan(-A[0, 1], -A[0, 2])
 #    return np.array([x, y, z])
 #
-#def zxz2euler(z1, x, z2):
+# def zxz2euler(z1, x, z2):
 #    return np.array([np.cos(.5 * (z1 + z2)) * np.cos(.5 * x),
 #                     np.cos(.5 * (z1 - z2)) * np.sin(.5 * x),
 #                     np.sin(.5 * (z1 - z2)) * np.sin(.5 * x),
 #                     np.sin(.5 * (z1 + z2)) * np.cos(.5 * x)])
 #
-#def xyz2A(x, y, z):
+# def xyz2A(x, y, z):
 #    return np.dot(Ax(x), np.dot(Ay(y), Az(z)))
 
-#def euler2xyz(euler):
+# def euler2xyz(euler):
 #    return A2xyz(euler2A(euler))
 
-#def A2euler(A):
+# def A2euler(A):
 #    return xyz2euler(*A2xyz(A))
 
 def euler2angle(euler):
@@ -111,5 +116,6 @@ def euler2angle(euler):
 
     return np.arccos(euler[0]) * 2
 
+
 def euler2gl(euler):
     return np.r_[deg(euler2angle(euler)), euler[1:]]

wetb/utils/geometry.py

-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import absolute_import
-from builtins import map
-from future import standard_library
-standard_library.install_aliases()
 import numpy as np
 
 def rad(deg):

wetb/utils/installer.py

+import glob
+import io
+import json
+import os
+import platform
+import shutil
+import sys
+import zipfile
+from pathlib import Path
+from platform import architecture
+from urllib.request import Request
+from urllib import request
+import ssl
+import certifi
+
+
+def urlopen(*args, **kwargs):
+    return request.urlopen(
+        *args, **kwargs, context=ssl.create_default_context(cafile=certifi.where())
+    )
+
+
+def chmod_x(exe_path: str):
+    """Utility function to change the file mode of a file to allow execution
+
+    Parameters
+    ----------
+    path : str
+        path to the file to apply the +x policy to.
+    """
+    st = os.stat(exe_path)
+    os.chmod(exe_path, st.st_mode | 0o111)
+
+
+def install_wind_tool(
+    tool: str = None, version: str = None, platform: str = None, destination: str = None
+):
+    """Utility function to download and install any tool available from the HAWC2 tools website
+
+    Parameters
+    ----------
+    tool : str, required
+        The tool to install. If None, the function will return immediately without doing anything. By default None
+    version : str, optional
+        Version of the software. If None, the latest version will be downloaded. By default None
+    platform : str, optional
+        The build platform of the tool. If None, the platform of the executing machine will be used. By default None
+    destination : str, optional
+        Destination path for the download / installation. If None, the destination is set to cwd. By default None
+    """
+
+    # Escape backslash if windowspath is given
+    destination = Path(destination.encode("unicode_escape").decode()).as_posix()
+
+    if tool is None:
+        print("No tool has been given for install. Nothing has been installed.")
+        return
+
+    # If platform is not set, get the platform of the executing machine
+    if platform is None:
+        if os.name == "nt":
+            if architecture()[0][:2] == "32":
+                platform = "win32"
+            else:
+                platform = "win64"
+        else:
+            platform = "linux"
+    else:
+        platform = platform.lower()
+
+    # Check if tool is available on the tools website
+    req = Request("https://tools.windenergy.dtu.dk/version_inventory.json")
+    versions = json.loads(urlopen(req).read())
+    if tool not in versions:
+        print(
+            f"The tool '{tool}' is not available in the inventory on tools.windenergy.dtu.dk."
+        )
+        return
+
+    # Check if requested version is available, and default it is not.
+    if version is not None and version not in versions[tool]["available"]:
+        print(
+            f"Version '{version}' of '{tool}' is not available - defaulting to the latest version: '{versions[tool]['latest']}'"
+        )
+        version = versions[tool]["latest"]
+    elif version is None:
+        version = versions[tool]["latest"]
+        print(f"Using latest version of '{tool}': {version}")
+
+    # Get a list of the different types of executables available for the given platform
+    req = Request("https://tools.windenergy.dtu.dk/product_inventory.json")
+    download = json.loads(urlopen(req).read())
+    exes = download[tool][version][platform]
+    types = [exes[app]["type"] for app in exes]
+
+    # Select the standalone executable if that is available, otherwise select the zip archive for download
+    for ind, file in enumerate(exes):
+        if "executable" in types[ind].lower() and "standalone" in types[ind].lower():
+            req = Request(f"https://tools.windenergy.dtu.dk{exes[file]['link']}")
+            break
+        elif file.split(".")[-1].lower() == "zip":
+            req = Request(f"https://tools.windenergy.dtu.dk{exes[file]['link']}")
+            break
+
+    # Download the file to a buffer object
+    buffer = io.BytesIO(urlopen(req).read())
+
+    # If destination is not provided, set the installation destination to cwd.
+    if destination is None:
+        destination = os.getcwd()
+
+    print(f"Installing {tool} version {version} for {platform}")
+    os.makedirs(destination, exist_ok=True)
+    # If the download is a zip archive, extract the archive to the destination directory
+    if req.full_url.endswith(".zip"):
+        zipfile.ZipFile(buffer).extractall(path=destination)
+    # If the download is not a zip archive, but an executable, save the bytes object to a file
+    else:
+        with open(f"{destination}/{req.full_url.split('/')[-1]}", "wb") as f:
+            f.write(buffer.getvalue())
+
+    # Add execution policy to the executable files (files with .exe or no extension)
+    for file in glob.glob(f"{destination}/*"):
+        if file.endswith(".exe") or os.path.splitext(file)[-1] == "":
+            chmod_x(file)
+
+    print(f"{tool} version {version} succesfully installed in {destination}")
+
+    return
+
+
+def install_hawc2_dtu_license():
+    """Function to install the DTU HAWC2 license. In order to install the license, you must be logged in to the DTU network."""
+    install_dtu_license("hawc2")
+
+
+def install_hawcstab2_dtu_license():
+    """Function to install the DTU HAWCStab2 license. In order to install the license, you must be logged in to the DTU network."""
+    install_dtu_license("hawcstab2")
+
+
+def install_ellipsys_dtu_license():
+    """Function to install the DTU HAWCStab2 license. In order to install the license, you must be logged in to the DTU network."""
+    install_dtu_license("ellipsys")
+
+
+def install_dtu_license(software : str):
+    """Function to install the DTU online license for HAWC2, HAWCStab2 and Ellipsys. In order to install the license, you must be logged in to the DTU network."""
+    software = software.lower()
+    assert software in ["hawc2", "hawcstab2", "ellipsys"], "Argument 'software' must be one of ['hawc2', 'hawcstab2,' 'ellipsys']"
+    if sys.platform.lower() == "win32":
+        f = Path(os.getenv("APPDATA")) / f"DTU Wind Energy/{software}/license.cfg"
+    else:
+        f = Path.home() / f".config/{software}/license.cfg"
+    if not f.exists():
+        f.parent.mkdir(parents=True, exist_ok=True)
+        r = urlopen("http://license-internal.windenergy.dtu.dk:34523").read()
+        if b"LICENSE SERVER RUNNING" in r:
+            f.write_text(
+                "[licensing]\nhost = http://license-internal.windenergy.dtu.dk\nport = 34523"
+            )
+        else:
+            raise ConnectionError(
+                f"Could not connect to the DTU license server. You must be connected to the DTU network to use this function."
+            )
+
+
+
+
+def install_keygen_license(software: str, cfg_file: str, force: bool = False):
+    """Install license file for HAWC2, HAWCStab2 or Ellipsys on your machine
+
+    Parameters
+    ----------
+    software : str
+        Must be one of HAWC2, HAWCStab2 or Ellipsys. The argument is case insensitive.
+    cfg_file : str
+        Path to the license file to install
+    force : bool, optional
+        Switch to force the installation, overwriting any existing , by default False
+
+    Returns
+    -------
+    NoneType
+        None
+
+    Raises
+    ------
+    ValueError
+        A ValueError is raised if the name of the software argument is not supported.
+    """
+
+    SUPPORTED_SOFTWARES = ["hawc2", "hawcstab2", "ellipsys"]
+    if software.lower() not in SUPPORTED_SOFTWARES:
+        raise ValueError(f"'software' must be one of {SUPPORTED_SOFTWARES}")
+
+    USER_PLATFORM = platform.uname().system
+
+    if USER_PLATFORM == "Windows":
+        APPDATA = f"{os.environ['APPDATA']}"
+    else:
+        APPDATA = "None"
+
+    def local_license_dir(platform, software):
+        return {
+            "Windows": os.path.join(
+                APPDATA,
+                "DTU Wind Energy",
+                f"{software}",
+            ),
+            "Linux": os.path.join(f"{Path.home()}", ".config", f"{software}"),
+            "Darwin": os.path.join(f"{Path.home()}", "Library", "Application Support"),
+        }[platform]
+
+    def local_license_file(software):
+        return {
+            "hawc2": "license.cfg",
+            "hawcstab2": "license.cfg",
+            "pywasp": "",
+            "ellipsys": "license.cfg",
+        }[software.lower()]
+
+    license_path = local_license_dir(USER_PLATFORM, software)
+    lic_name = local_license_file(software)
+
+    os.makedirs(license_path, exist_ok=True)
+
+    if (
+        os.path.exists(os.path.join(license_path, lic_name))
+        and os.path.isfile(os.path.join(license_path, lic_name))
+        and (not force)
+    ):
+        print(
+            f"License already installed for {software}, use 'force=True' to overwrite installation"
+        )
+    else:
+        shutil.copy(f"{cfg_file}", f"{os.path.join(license_path,lic_name)}")

wetb/utils/plot_rotation.py

+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+from wetb.utils.rotation import axis2axis_angle, axis_angle2quaternion, quaternion2matrix, s2matrix, matrix2quaternion,\
+    quaternion2axis_angle, axis_angle2axis, matrix2axis
+
+
+def init_plot(ref_coo='gl', rot=(180, 180)):
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+    plot_coo(ref_coo, np.eye(3) * 1.5, lw=1, alpha=1)
+    for xyz in 'xyz':
+        getattr(ax, 'set_%slim' % xyz)([-2, 2])
+        getattr(ax, 'set_%slabel' % xyz)(xyz)
+
+    ax.view_init(rot[0] + 10, rot[1] - 10)
+
+
+def plot_coo(coo, A, origo=(0, 0, 0), lw=1, alpha=1):
+    ax = plt.gca()
+    for v, c, xyz in zip(np.asarray(A).T, 'rgb', 'xyz'):
+        ax.quiver(*origo, *v, color=c, lw=lw, alpha=alpha, arrow_length_ratio=.2, zorder=-32)
+        ax.text(*(np.array(origo) + v * 1.1), "$%s_{%s}$" % (xyz, coo),
+                alpha=alpha, fontsize=12, fontweight='bold', zorder=32)
+
+
+def plot_axis_rotation(rot_coo, axis, origo=(0, 0, 0), ref_coo=np.eye(3), deg=False):
+
+    axis_angle = axis2axis_angle(axis)
+    mat = quaternion2matrix(axis_angle2quaternion(axis_angle, deg))
+    mat = np.dot(ref_coo, mat)
+    axis = np.dot(ref_coo, axis_angle[:3])
+    xyz = np.array([origo, np.array(origo) + axis * 1.5])
+    plt.gca().plot(*xyz.T, '-.k', lw=1)
+    plot_coo(rot_coo, mat, origo)
+    return mat
+
+
+def plot_matrix_rotation(rot_name, matrix, origo, ref_coo=np.eye(3)):
+    mat = np.dot(ref_coo, matrix)
+    plot_coo(rot_name, mat, origo)
+
+
+def plot_body(nodes):
+    ax = plt.gca()
+    ax.plot(*np.array(nodes).T, alpha=0.2, lw=10)
+    ax.plot(*np.array(nodes).T, '.k', alpha=1)
+
+
+def set_aspect_equal():
+    ax = plt.gca()
+    xyz_lim = [getattr(ax, 'get_%slim' % xyz)() for xyz in 'xyz']
+    max_range = np.max([np.abs(np.subtract(*lim)) for lim in xyz_lim]) / 2
+    mid = np.mean(xyz_lim, 1)
+    for xyz, m in zip('xyz', mid):
+        getattr(ax, "set_%slim" % xyz)([m - max_range, m + max_range])
+    plt.tight_layout()
+
+
+def vs_array(s, shape):
+    arr = np.array(s.strip().replace("D", "E").split("\n"), dtype=float)
+    if len(np.atleast_1d(shape)) == 2:
+        arr = arr.reshape(shape[::-1]).T
+    return arr
+
+
+if __name__ == '__main__':
+    TBD = vs_array("""-0.991775138730704
+-1.644972253859245D-002
+0.126931008126843
+1.644972253859245D-002
+0.967100554922815
+0.253862016253685
+-0.126931008126843
+0.253862016253685
+-0.958875693653519""", (3, 3))
+    TUB = vs_array("""-1.00000000000000
+0.000000000000000D+000
+0.000000000000000D+000
+0.000000000000000D+000
+1.00000000000000
+0.000000000000000D+000
+0.000000000000000D+000
+0.000000000000000D+000
+-1.00000000000000""", (3, 3))
+    TDU = vs_array("""0.991775138730704
+-1.644972253859245D-002
+0.126931008126843
+-1.644972253859245D-002
+0.967100554922815
+0.253862016253685
+-0.126931008126843
+-0.253862016253685
+0.958875693653519""", (3, 3))
+    print(np.dot(TBD.T, TUB.T) - TDU)
+    print(TDU)
+    rot_bd = matrix2axis(TBD)
+    rot_bu = matrix2axis(TUB)
+    rot_du = matrix2axis(TDU)
+
+    print(np.dot(TUB, matrix2axis(np.dot(TUB.T, TBD.T))))
+    print(matrix2axis(TDU))
+
+    print(np.dot(TUB.T, TBD.T).T)

wetb/utils/postprocs.py

+# -*- coding: utf-8 -*-
+"""
+Created on Wed Sep 18 11:21:18 2024
+
+@author: nicgo
+"""
+
+import numpy as np
+import xarray as xr
+from wetb.fatigue_tools.fatigue import eq_load, cycle_matrix
+from wetb.utils.envelope import projected_extremes
+from wetb.utils.rotation import projection_2d
+                                        
+def statistics(time, data, info,
+               statistics=['min', 'mean', 'max', 'std']) -> xr.DataArray:
+    """
+    Calculate statistics from different sensors.
+
+    Parameters
+    ----------
+    time : array (Ntimesteps)
+        Array containing the simulation time from start to end of writting output
+    data : array (Ntimesteps x Nsensors)
+        Array containing the time series of all sensors
+    info : dict
+        Dictionary that must contain the following entries:\n
+            attribute_names: list of sensor names\n
+            attribute_units: list of sensor units\n
+            attribute_descriptions: list of sensor descriptions
+    statistics : list (Nstatistics), optional
+        List containing the different types of statistics to calculate for each sensor.
+        The default is ['min', 'mean', 'max', 'std'].
+
+    Returns
+    -------
+    xarray.DataArray (Nsensors x Nstatistics)
+        DataArray containing statistics for all sensors.\n
+        Dims: sensor_name, statistic\n
+        Coords: sensor_name, statistic, sensor_unit, sensor_description
+
+    """
+    def get_stat(stat):
+        if hasattr(np, stat):
+            return getattr(np, stat)(data, 0)
+    data = np.array([get_stat(stat) for stat in statistics]).T
+    dims = ['sensor_name', 'statistic']
+    coords = {'statistic': statistics,
+              'sensor_name': info['attribute_names'],
+              'sensor_unit': ('sensor_name', info['attribute_units']),
+              'sensor_description': ('sensor_name', info['attribute_descriptions'])}
+    return xr.DataArray(data=data, dims=dims, coords=coords)
+
+
+def extreme_loads(data, sensors_info) -> xr.DataArray:
+    """
+    Calculate the extreme load matrix (Fx, Fy, Fz, Mx, My, Mz) for different sensors
+    using as criteria driving force or moment along x, y and z directions,
+    including also maximum x-y in-plane loads Fres and Mres.
+
+    Parameters
+    ----------
+    data : array (Ntimesteps x Nsensors_all)
+        Array containing the time series of all sensors
+    sensors_info : list of tuples (Nsensors)
+        List of 7-element tuples containing the name of a load sensor and the indices of the 6 of their components.\n
+        Example: [('Tower base', 0, 1, 2, 3, 4, 5), ('Blade 1 root', 6, 7, 8, 9, 10, 11)]
+
+    Returns
+    -------
+    xarray.DataArray (Nsensors x 14)
+        DataArray containing extreme load states (Fx, Fy, Fz, Mx, My, Mz) for each sensor in sensors_info
+        for each case where a load is maximum or minimum (6 x 2 = 12) plus the 2 cases where Fres and Mres are maximum.\n
+        Dims: sensor_name, driver, load\n
+        Coords: sensor_name, driver, load, sensor_unit
+
+    """
+    extreme_loads_data = []
+    for component, i_fx, i_fy, i_fz, i_mx, i_my, i_mz in sensors_info:
+        extreme_loads_data.append([])
+        for i in i_fx, i_fy, i_fz, i_mx, i_my, i_mz:
+            time_step_max = np.argmax(data[:, i])
+            time_step_min = np.argmin(data[:, i])
+            extreme_loads_data[-1].append([data[time_step_max, i_fx], data[time_step_max, i_fy], data[time_step_max, i_fz],
+                                     data[time_step_max, i_mx], data[time_step_max, i_my], data[time_step_max, i_mz]])
+            extreme_loads_data[-1].append([data[time_step_min, i_fx], data[time_step_min, i_fy], data[time_step_min, i_fz],
+                                     data[time_step_min, i_mx], data[time_step_min, i_my], data[time_step_min, i_mz]])
+        fres = np.sqrt(data[:, i_fx]**2 + data[:, i_fy]**2)
+        mres = np.sqrt(data[:, i_mx]**2 + data[:, i_my]**2)
+        time_step_max_fres = np.argmax(fres)
+        time_step_max_mres = np.argmax(mres)
+        extreme_loads_data[-1].append([data[time_step_max_fres, i_fx], data[time_step_max_fres, i_fy], data[time_step_max_fres, i_fz],
+                                 data[time_step_max_fres, i_mx], data[time_step_max_fres, i_my], data[time_step_max_fres, i_mz]])
+        extreme_loads_data[-1].append([data[time_step_max_mres, i_fx], data[time_step_max_mres, i_fy], data[time_step_max_mres, i_fz],
+                                 data[time_step_max_mres, i_mx], data[time_step_max_mres, i_my], data[time_step_max_mres, i_mz]])       
+    data = np.array(extreme_loads_data)
+    dims = ['sensor_name', 'driver', 'load']
+    coords = {'sensor_name': [s[0] for s in sensors_info],
+              'driver': ['Fx_max', 'Fx_min', 'Fy_max', 'Fy_min', 'Fz_max', 'Fz_min',
+                         'Mx_max', 'Mx_min', 'My_max', 'My_min', 'Mz_max', 'Mz_min',
+                         'Fres_max', 'Mres_max'],
+              'load': ['Fx', 'Fy', 'Fz', 'Mx', 'My', 'Mz'],
+              'sensor_unit': ('load', ['kN', 'kN', 'kN', 'kNm', 'kNm', 'kNm']),
+                }
+    return xr.DataArray(data=data, dims=dims, coords=coords)
+
+
+def directional_extreme_loads(data, info, sensors_info, angles=np.linspace(-150,180,12), sweep_angle=30, degrees=True) -> xr.DataArray:
+    """
+    Calculate extreme loads from different sensors along different directions.
+
+    Parameters
+    ----------
+    data : array (Ntimesteps x Nsensors_all)
+        Array containing the time series of all sensors
+    info : dict
+        Dictionary that must contain the following entries:\n
+            attribute_names: list of sensor names\n
+            attribute_units: list of sensor units
+    sensors_info : list of tuples (Nsensors)
+        List of 3-element tuples containing the name of a load sensor and the indices of 2 of their components.\n
+        Example: [('Tower base shear force', 0, 1), ('Blade 1 root bending moment', 9, 10)]
+    angles : array (Ndirections), optional
+        Array containing the directions along which extreme loads should be computed.
+        The default is np.linspace(-150,180,12).
+    sweep_angle : float, optional
+        Angle to be swept to both sides of a direction in the search for its extreme load.
+        It should match the angular step in angles. The default is 30.
+    degrees : bool, optional
+        Whether angles and sweep_angle are in degrees (True) or radians (False). The default is True.
+
+    Returns
+    -------
+    xarray.DataArray (Nsensors x Ndirections)
+        DataArray containing extreme loads for each sensor in sensors_info in each direction in angles.\n
+        Dims: sensor_name, angle\n
+        Coords: sensor_name, angle, sensor_unit
+
+    """
+    directional_extreme_loads_data = []
+    for sensor, ix, iy in sensors_info:
+        directional_extreme_loads_data.append(projected_extremes(np.vstack([data[:,ix],data[:,iy]]).T, angles, sweep_angle, degrees)[:,1])
+    data = np.array(directional_extreme_loads_data)
+    dims = ['sensor_name', 'angle']
+    coords = {'angle': angles,
+              'sensor_name': [s[0] for s in sensors_info],
+              'sensor_unit': ('sensor_name', [info['attribute_units'][s[1]] for s in sensors_info]),
+               }
+    return xr.DataArray(data=data, dims=dims, coords=coords)
+
+
+def extremes_and_contemporaneous(data, info, sensor_indices) -> xr.DataArray:
+    """
+    Calculate the maximum and minimum values for different sensors as well as
+    the contemporaneous values of the other sensors at such instants.
+    
+    Parameters
+    ----------
+    data : array (Ntimesteps x Nsensors_all)
+        Array containing the time series of all sensors
+    sensor_indices : list
+        List containing the indices of the desired output sensors
+
+    Returns
+    -------
+    xarray.DataArray ((Nsensors*2) x Nsensors)
+        DataArray containing the maximum and minimum values for each sensor
+        and the contemporaneous values of the other sensors.
+    """
+    table = []
+    for s in sensor_indices:
+        time_step_max = np.argmax(data[:, s])
+        time_step_min = np.argmin(data[:, s])
+        table.append([data[time_step_max, s2] for s2 in sensor_indices])
+        table.append([data[time_step_min, s2] for s2 in sensor_indices])
+    data = np.array(table)
+    dims = ['driver', 'sensor_description']
+    sensor_descriptions = [info['attribute_descriptions'][s] for s in sensor_indices]
+    coords = {'driver': [s + '_' + e for s in sensor_descriptions for e in ['max', 'min']],
+              'sensor_description': sensor_descriptions,
+              'sensor_name': ('sensor_description', [info['attribute_names'][s] for s in sensor_indices]),
+              'sensor_unit': ('sensor_description', [info['attribute_units'][s] for s in sensor_indices])}
+    return xr.DataArray(data=data, dims=dims, coords=coords)
+
+
+def equivalent_loads(data, time, info, m_list=[3, 4, 6, 8, 10, 12], neq=None, no_bins=46) -> xr.DataArray:
+    """
+    Calculate fatigue equivalent loads for different sensors
+    for different Woehler slopes.
+
+    Parameters
+    ----------
+    data : array (Ntimesteps x Nsensors)
+        Array containing the time series of all sensors
+    info : dict
+        Dictionary that must contain the following entries:\n
+            attribute_names: list of sensor names\n
+            attribute_units: list of sensor units
+    m_list : list (Nwoehlerslopes), optional
+        List containing the different woehler slopes. The default is [3, 4, 6, 8, 10, 12].
+    neq : int or float, optional
+        Number of equivalent load cycles. The default is the time duration in seconds.
+    no_bins : int, optional
+        Number of bins in rainflow count histogram. The default is 46.
+
+    Returns
+    -------
+    xarray.DataArray (Nsensors x Nwoehlerslopes)
+        DataArray containing fatigue equivalent loads for
+        each Woehler slope in m_list.\n
+        Dims: sensor_name, m\n
+        Coords: sensor_name, sensor_unit, sensor_description, m
+
+    """
+    if neq is None:
+        neq = time[-1] - time[0] + time[1] - time[0]
+    equivalent_loads_data = [eq_load(sensor, no_bins=no_bins, m=m_list, neq=neq)[0] for sensor in data.T]
+    data = np.array(equivalent_loads_data)
+    dims = ['sensor_name', 'm']
+    coords = {'sensor_name': info['attribute_names'],
+              'sensor_unit': ('sensor_name', info['attribute_units']),
+              'sensor_description': ('sensor_name', info['attribute_descriptions']),
+              'm': m_list,
+               }
+    return xr.DataArray(data=data, dims=dims, coords=coords)
+
+
+def directional_equivalent_loads(data, time, info, sensors_info, m_list=[3, 4, 6, 8, 10, 12], neq=None, no_bins=46, angles=np.linspace(-150,180,12), degrees=True) -> xr.DataArray:
+    """
+    Calculate directional fatigue equivalent loads for different load sensors
+    for different Woehler slopes for different directions.
+
+    Parameters
+    ----------
+    data : array (Ntimesteps x Nsensors_all)
+        Array containing the time series of all sensors
+    info : dict
+        Dictionary that must contain the following entries:\n
+            attribute_names: list of sensor names\n
+            attribute_units: list of sensor units
+    sensors_info : list of tuples (Nsensors)
+        List of 3-element tuples containing the name of a load sensor and the indices of 2 of their components.\n
+        Example: [('Tower base shear force', 0, 1), ('Blade 1 root bending moment', 9, 10)]
+    m_list : list (Nwoehlerslopes), optional
+        List containing the different woehler slopes. The default is [3, 4, 6, 8, 10, 12].
+    neq : int or float, optional
+        Number of equivalent load cycles. The default is the time duration in seconds.
+    no_bins : int, optional
+        Number of bins in rainflow count histogram. The default is 46.
+    angles : array (Ndirections), optional
+        Array containing the directions along which fatigue loads should be computed.
+        The default is np.linspace(-150,180,12).
+    degrees : bool, optional
+        Whether angles are in degrees (True) or radians (False). The default is True.
+
+    Returns
+    -------
+    xarray.DataArray (Nsensors x Ndirections x Nwoehlerslopes)
+        DataArray containing directional fatigue loads for each sensor in sensors_info, for each Woehler slope in m_list
+        and for each direction in angles.\n
+        Dims: sensor_name, angle, m\n
+        Coords: sensor_name, angle, m, sensor_unit
+
+    """
+    if neq is None:
+        neq = time[-1] - time[0] + time[1] - time[0]
+    directional_equivalent_loads_data = []
+    for sensor, ix, iy in sensors_info:
+        directional_equivalent_loads_data.append([])
+        for angle in angles:
+            directional_equivalent_loads_data[-1].append([])
+            for m in m_list:
+                directional_equivalent_loads_data[-1][-1].append(eq_load(data[:, [ix, iy]] @ projection_2d(angle, degrees=degrees),
+                                             no_bins=no_bins,
+                                             m=m,
+                                             neq=neq)[0][0])
+    data = np.array(directional_equivalent_loads_data)
+    dims = ['sensor_name', 'angle', 'm']
+    coords = {'sensor_name': [s[0] for s in sensors_info],
+              'angle': angles,
+              'm': m_list,
+              'sensor_unit': ('sensor_name', [info['attribute_units'][s[1]] for s in sensors_info]),
+               }
+    return xr.DataArray(data=data, dims=dims, coords=coords)
+
+
+def markov_matrices(data, info, no_bins=46) -> xr.DataArray:
+    """
+    Calculate the Markov matrices for different sensors.
+
+    Parameters
+    ----------
+    data : array (Ntimesteps x Nsensors)
+        Array containing the time series of all sensors
+    info : dict
+        Dictionary that must contain the following entries:\n
+            attribute_names: list of sensor names\n
+            attribute_units: list of sensor units
+    no_bins : int, optional
+        Number of bins in rainflow count histogram. The default is 46.
+
+    Returns
+    -------
+    xarray.DataArray (Nsensors x Nbins x 2)
+        DataArray containing number of cycles and load amplitude 
+        for each sensor.\n
+        Dims: sensor_name, bin, (cycles, amplitude)\n
+        Coords: sensor_name, sensor_unit, sensor_description
+
+    """
+    cycles_and_amplitudes = []
+    for sensor in data.T:
+        try:
+            cycles, ampl_bin_mean, _, _, _ = cycle_matrix(sensor,
+                                                          ampl_bins=no_bins,
+                                                          mean_bins=1)
+            cycles_and_amplitudes.append([[cycles.flatten()[i],
+                                         ampl_bin_mean.flatten()[i]]
+                                         for i in range(no_bins)])
+        except TypeError:
+            cycles_and_amplitudes.append([[np.nan, np.nan]]*no_bins)
+    data = np.array(cycles_and_amplitudes)
+    dims = ['sensor_name', 'bin', '(cycles, amplitude)']
+    coords = {'sensor_name': info['attribute_names'],
+              'sensor_unit': ('sensor_name', info['attribute_units']),
+              'sensor_description': ('sensor_name', info['attribute_descriptions']),
+               }
+    return xr.DataArray(data=data, dims=dims, coords=coords)
+
+
+def directional_markov_matrices(data, info, sensors_info, no_bins=46, angles=np.linspace(-150,180,12), degrees=True) -> xr.DataArray:
+    """
+    Calculate the Markov matrices for different load sensors for different directions.
+
+    Parameters
+    ----------
+    data : array (Ntimesteps x Nsensors_all)
+        Array containing the time series of all sensors
+    info : dict
+        Dictionary that must contain the following entries:\n
+            attribute_names: list of sensor names\n
+            attribute_units: list of sensor units
+    sensors_info : list of tuples (Nsensors)
+        List of 3-element tuples containing the name of a load sensor and the indices of 2 of their components.\n
+        Example: [('Tower base shear force', 0, 1), ('Blade 1 root bending moment', 9, 10)]
+    no_bins : int, optional
+        Number of bins in rainflow count histogram. The default is 46.
+    angles : array (Ndirections), optional
+        Array containing the directions along which fatigue loads should be computed.
+        The default is np.linspace(-150,180,12).
+    degrees : bool, optional
+        Whether angles are in degrees (True) or radians (False). The default is True.
+
+    Returns
+    -------
+    xarray.DataArray (Nsensors x Ndirections x Nbins x 2)
+        DataArray containing directional number of cycles and load amplitude 
+        for each sensor in sensors_info, for each direction in angles and
+        for each bin in no_bins.\n
+        Dims: sensor_name, angle, bin, (cycles, amplitude)\n
+        Coords: sensor_name, angle, sensor_unit
+
+    """
+    cycles_and_amplitudes = []
+    for sensor, ix, iy in sensors_info:
+        cycles_and_amplitudes.append([])
+        for angle in angles:
+            cycles, ampl_bin_mean, _, _, _ = cycle_matrix(data[:, [ix, iy]] @ projection_2d(angle, degrees=degrees),
+                                         ampl_bins=no_bins,
+                                         mean_bins=1)
+            cycles_and_amplitudes[-1].append([[cycles.flatten()[i],
+                                             ampl_bin_mean.flatten()[i]]
+                                             for i in range(no_bins)])
+    data = np.array(cycles_and_amplitudes)
+    dims = ['sensor_name', 'angle', 'bin', '(cycles, amplitude)']
+    coords = {'sensor_name': [s[0] for s in sensors_info],
+              'angle': angles,
+              'sensor_unit': ('sensor_name', [info['attribute_units'][s[1]] for s in sensors_info]),
+               }
+    return xr.DataArray(data=data, dims=dims, coords=coords)
+
+    
+def remove_postproc(file_h5py, remove_all=False, postproc_list=[]) -> None:
+    """
+    Remove postproc data from hdf5 file.
+
+    Parameters
+    ----------
+    file_h5py : h5py.File object
+        h5py.File object in append mode from hdf5 file
+    remove_all : bool, optional
+        Whether all postprocs should be removed from hdf5 file. The default is False.
+    postproc_list : list of functions, optional
+        List containing which postproc functions to remove their data from hdf5 file. The default is [].
+
+    Returns
+    -------
+    None.
+
+    """
+    if remove_all:
+        try:
+            del file_h5py['postproc']
+            print("All postprocs removed'")
+        except KeyError:
+            print("No postprocs already")
+    else:
+        for postproc in postproc_list:
+            try:
+                del file_h5py['postproc'][postproc.__name__]
+                print(f"{postproc.__name__} removed'")
+            except KeyError:
+                print(f"No {postproc.__name__} already'") 

wetb/utils/process_exec.py

@@ -3,26 +3,26 @@ Created on 10/03/2014
 
 @author: MMPE
 '''
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import division
-from __future__ import absolute_import
-from builtins import range
-from builtins import str
-from future import standard_library
 import glob
 import re
-standard_library.install_aliases()
 
 import os
 
 
 DEBUG = False
-def pexec(args, cwd=None):
+
+
+def pexec(args, cwd=None, shell=False):
     """
     usage: errorcode, stdout, stderr, cmd = pexec("MyProgram.exe arg1, arg2", r"c:\tmp\")
 
     """
+    proc, cmd = process(args, cwd)
+    errorcode, stdout, stderr = exec_process(proc)
+    return errorcode, stdout, stderr, cmd
+
+
+def process(args, cwd=None, shell=False):
     import subprocess
     if not isinstance(args, (list, tuple)):
         args = [args]
@@ -31,33 +31,24 @@ def pexec(args, cwd=None):
         if os.path.exists(args[i]):
             args[i] = str(args[i]).replace('/', os.path.sep).replace('\\', os.path.sep).replace('"', '')
 
-    cmd = "%s" % '{} /c "{}"'.format (os.environ.get("COMSPEC", "cmd.exe"), subprocess.list2cmdline(args))
+    cmd_args = subprocess.list2cmdline(args)
     if cwd and os.path.isfile(cwd):
         cwd = os.path.dirname(cwd)
-    proc = subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True, cwd=cwd)
-    stdout, stderr = proc.communicate()
-    errorcode = proc.returncode
+    if cwd:
+        cmd_cwd = "cd %s && " % cwd
+    else:
+        cmd_cwd = ""
 
-    return errorcode, stdout.decode('cp1252'), stderr.decode('cp1252'), cmd
+    if os.name == 'nt':
+        cmd = '%s /c "%s%s"' % (os.environ.get("COMSPEC", "cmd.exe"), cmd_cwd, cmd_args)
+    else:
+        cmd = '%s%s' % (cmd_cwd, cmd_args)
 
+    return subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=shell, cwd=cwd), cmd
 
-def process(args, cwd=None):
-    import subprocess
-    if not isinstance(args, (list, tuple)):
-        args = [args]
-    args = [str(arg) for arg in args]
-    for i in range(len(args)):
-        if os.path.exists(args[i]):
-            args[i] = str(args[i]).replace('/', os.path.sep).replace('\\', os.path.sep).replace('"', '')
-
-    cmd = "%s" % '{} /c "{}"'.format (os.environ.get("COMSPEC", "cmd.exe"), subprocess.list2cmdline(args))
-    if cwd is not None and os.path.isfile(cwd):
-        cwd = os.path.dirname(cwd)
-    return subprocess.Popen(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=False, cwd=cwd)
 
 def exec_process(process):
     stdout, stderr = process.communicate()
     errorcode = process.returncode
 
     return errorcode, stdout.decode(), stderr.decode()
-

wetb/utils/rotation.py

-'''
+"""
 Created on 21/07/2014
 
 @author: mmpe
-'''
+"""
 import numpy as np
+import numpy.typing as npt
+
+
+def s2matrix(s):
+    """Creates a transformation matrix from string
+
+    Parameters
+    ----------
+    s : string
+        x is
+    Returns
+    -------
+    matrix : array_like
+        3x3 transformation matrix
+
+    Examples
+    --------
+    >> s2matrix('x,y,z') # identity matrix
+    >> s2matrix('x,-z,y) # 90 deg rotation around x
+
+    """
+    d = {xyz: v for xyz, v in zip("xyz", np.eye(3))}
+    return np.array([eval(xyz, d) for xyz in s.split(",")]).T
 
 
 def transformation_matrix(angles, xyz):
@@ -26,21 +49,26 @@ def transformation_matrix(angles, xyz):
     transformation_matrix : array_like, shape = (no_angles, 3,3)
         Rotation matrix(es)
     """
-    indexes = [(0, 4, 8, 5, 7), (4, 0, 8, 6, 2), (8, 0, 4, 1, 3)]  # 1, cos,cos,sin,-sin for x,y and z rotation matrix
+    indexes = [
+        (0, 4, 8, 5, 7),
+        (4, 0, 8, 6, 2),
+        (8, 0, 4, 1, 3),
+    ]  # 1, cos,cos,sin,-sin for x,y and z rotation matrix
     if isinstance(angles, (int, float)):
         n = 1
     else:
         n = len(angles)
-    m = np.zeros(n * 9, np.float)
+    m = np.zeros(n * 9, float)
     cosx = np.cos(angles)
     sinx = np.sin(angles)
-    m[indexes[xyz][0]::9] = 1
-    m[indexes[xyz][1]::9] = cosx
-    m[indexes[xyz][2]::9] = cosx
-    m[indexes[xyz][3]::9] = sinx
-    m[indexes[xyz][4]::9] = -sinx
+    m[indexes[xyz][0] :: 9] = 1
+    m[indexes[xyz][1] :: 9] = cosx
+    m[indexes[xyz][2] :: 9] = cosx
+    m[indexes[xyz][3] :: 9] = sinx
+    m[indexes[xyz][4] :: 9] = -sinx
     return m.reshape(n, 3, 3)
 
+
 def rotmat(angles, xyz):
     """Create rotation matrix(es)
     !!!Note that the columns of the returned matrix(es) is rotated xyz-axes in original(unrotated) coordinates\n
@@ -61,21 +89,26 @@ def rotmat(angles, xyz):
     transformation_matrix : array_like, shape = (no_angles, 3,3)
         Rotation matrix(es)
     """
-    indexes = [(0, 4, 8, 5, 7), (4, 0, 8, 6, 2), (8, 0, 4, 1, 3)]  # 1, cos,cos,sin,-sin for x,y and z rotation matrix
+    indexes = [
+        (0, 4, 8, 5, 7),
+        (4, 0, 8, 6, 2),
+        (8, 0, 4, 1, 3),
+    ]  # 1, cos,cos,sin,-sin for x,y and z rotation matrix
     if isinstance(angles, (int, float)):
         n = 1
     else:
         n = len(angles)
-    m = np.zeros(n * 9, np.float)
+    m = np.zeros(n * 9, float)
     cosx = np.cos(angles)
     sinx = np.sin(angles)
-    m[indexes[xyz][0]::9] = 1
-    m[indexes[xyz][1]::9] = cosx
-    m[indexes[xyz][2]::9] = cosx
-    m[indexes[xyz][3]::9] = -sinx
-    m[indexes[xyz][4]::9] = sinx
+    m[indexes[xyz][0] :: 9] = 1
+    m[indexes[xyz][1] :: 9] = cosx
+    m[indexes[xyz][2] :: 9] = cosx
+    m[indexes[xyz][3] :: 9] = -sinx
+    m[indexes[xyz][4] :: 9] = sinx
     return m.reshape(n, 3, 3)
 
+
 def mdot(m1, m2):
     """Multiplication of matrix pairs
 
@@ -108,7 +141,9 @@ def mdot(m1, m2):
         for i in range(m1.shape[0]):
             mprod[i, :] = np.dot(m1[i, :], m2[i, :])
     else:
-        raise Exception("m1 and m2 must have same first dimension or one of the matrices must have first dimension 1")
+        raise Exception(
+            "m1 and m2 must have same first dimension or one of the matrices must have first dimension 1"
+        )
     return mprod
 
 
@@ -144,7 +179,7 @@ def dots(rotmats, v):
     else:
         if m.shape[0] != v.shape[1]:
             raise Exception("m must have same dimension as v has number of columns")
-        res = np.zeros_like(v, dtype=np.float)
+        res = np.zeros_like(v, dtype=float)
         for i in range(v.shape[1]):
             res[:, i] = np.dot(m[i], v[:, i])
         return res
@@ -160,10 +195,16 @@ def rotate(rotmats, v):
             return rotate(rotmats_red, v)
         else:
             rotmats = rotmats[0]
-    assert rotmats.shape[-1] == rotmats.shape[-2] == v.shape[-1] == 3, "rotmats is %d x %x and v is i x %d, but must be 3" % (rotmats.shape[-1], rotmats.shape[-2], v.shape[-1])
+    assert (
+        rotmats.shape[-1] == rotmats.shape[-2] == v.shape[-1] == 3
+    ), "rotmats is %d x %x and v is i x %d, but must be 3" % (
+        rotmats.shape[-1],
+        rotmats.shape[-2],
+        v.shape[-1],
+    )
 
-#    if isinstance(v, tuple):
-#        v = np.array([v])
+    #    if isinstance(v, tuple):
+    #        v = np.array([v])
     n = (1, v.shape[0])[len(v.shape) == 2]
 
     if rotmats.shape[0] == 1:
@@ -172,11 +213,12 @@ def rotate(rotmats, v):
         if v.shape[0] != rotmats.shape[0]:
             raise Exception("V and rotmats must have same first dimension")
         v = v.T
-        v_rot = np.zeros_like(v, dtype=np.float)
+        v_rot = np.zeros_like(v, dtype=float)
         for i in range(n):
             v_rot[:, i] = np.dot(rotmats[i], v[:, i])
         return v_rot.T
 
+
 def rotate_x(v, angle):
     """Rotate vector(s) around x axis
 
@@ -192,7 +234,10 @@ def rotate_x(v, angle):
     y : array_like
         Rotated vector(s)
     """
-    return _rotate(v, angle, lambda x, y, z, cos, sin : [x, cos * y - sin * z, sin * y + cos * z])
+    return _rotate(
+        v, angle, lambda x, y, z, cos, sin: [x, cos * y - sin * z, sin * y + cos * z]
+    )
+
 
 def rotate_y(v, angle):
     """Rotate vector(s) around y axis
@@ -209,7 +254,10 @@ def rotate_y(v, angle):
     y : array_like
         Rotated vector(s)
     """
-    return _rotate(v, angle, lambda x, y, z, cos, sin : [cos * x + sin * z, y, -sin * x + cos * z])
+    return _rotate(
+        v, angle, lambda x, y, z, cos, sin: [cos * x + sin * z, y, -sin * x + cos * z]
+    )
+
 
 def rotate_z(v, angle):
     """Rotate vector(s) around z axis
@@ -226,7 +274,10 @@ def rotate_z(v, angle):
     y : array_like
         Rotated vector(s)
     """
-    return _rotate(v, angle, lambda x, y, z, cos, sin : [cos * x - sin * y, sin * x + cos * y, z])
+    return _rotate(
+        v, angle, lambda x, y, z, cos, sin: [cos * x - sin * y, sin * x + cos * y, z]
+    )
+
 
 def _rotate(v, angle, rotfunc):
     angle = np.atleast_1d(angle)
@@ -235,8 +286,242 @@ def _rotate(v, angle, rotfunc):
     if len(v.shape) == 1:
         assert angle.shape[0] == 1
         assert v.shape[0] == 3
-        return np.array(rotfunc(v[0], v[1], v[2], cos, sin), dtype=np.float).T
+        return np.array(rotfunc(v[0], v[1], v[2], cos[0], sin[0]), dtype=float).T
     else:
         assert angle.shape[0] == 1 or angle.shape[0] == v.shape[0]
         assert v.shape[1] == 3
-        return np.array(rotfunc(v[:, 0], v[:, 1], v[:, 2], cos, sin), dtype=np.float).T
+        return np.array(rotfunc(v[:, 0], v[:, 1], v[:, 2], cos, sin), dtype=float).T
+
+
+# =======================================================================================================================
+# Conversions
+# =======================================================================================================================
+# http://www.euclideanspace.com/maths/geometry/rotations/conversions/
+# - axis: [x,y,z]*angle_deg
+# - axis_angle: [x,y,z,angle_rad]
+# - Quaternion: [qw,qx,qy,qz]
+# - Matrix 3x3 transformation matrix
+
+
+def norm(vector):
+    return np.sqrt(np.sum(np.asarray(vector) ** 2))
+
+
+# =======================================================================================================================
+# axis to ...
+# =======================================================================================================================
+def axis2axis_angle(axis):
+    """
+    deg : boolean
+        if True, axis length is assumed to be angle in deg
+    """
+    axis = np.asarray(axis)
+    angle = np.sqrt(((axis**2).sum()))
+    return np.r_[axis / np.sqrt((axis**2).sum()), angle]
+
+
+def axis2matrix(axis, deg=False):
+    # http://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToMatrix/index.htm
+    axis = np.asarray(axis)
+    angle = np.sqrt(((axis**2).sum()))
+    if deg:
+        angle = np.deg2rad(angle)
+    if angle == 0:
+        return np.eye(3)
+    x, y, z = xyz = axis / np.sqrt((axis**2).sum())
+
+    c, s = np.cos(angle), np.sin(angle)
+    t = 1 - c
+
+    return np.array(
+        [
+            [t * x * x + c, t * x * y - z * s, t * x * z + y * s],
+            [t * x * y + z * s, t * y * y + c, t * y * z - x * s],
+            [t * x * z - y * s, t * y * z + x * s, t * z * z + c],
+        ]
+    )
+    # alternative implementation
+    # asym = np.array([[0, -z, y], [z, 0, -x], [-y, x, 0]])
+    # return c * np.eye(3, 3) + s * asym + (1 - c) * xyz[np.newaxis] * xyz[:, np.newaxis]
+
+
+def axis2quaternion(axis):
+    axis = np.asarray(axis)
+    q = np.zeros([axis.shape[0], 4])
+
+    if len(axis.shape) > 1:
+        yaw = axis[:, 2]
+        pitch = axis[:, 1]
+        roll = axis[:, 0]
+    else:
+        yaw = axis[2]
+        pitch = axis[1]
+        roll = axis[0]
+
+    cy = np.cos(yaw * 0.5)
+    sy = np.sin(yaw * 0.5)
+    cp = np.cos(pitch * 0.5)
+    sp = np.sin(pitch * 0.5)
+    cr = np.cos(roll * 0.5)
+    sr = np.sin(roll * 0.5)
+
+    qw = cr * cp * cy + sr * sp * sy
+    qx = sr * cp * cy - cr * sp * sy
+    qy = cr * sp * cy + sr * cp * sy
+    qz = cr * cp * sy - sr * sp * cy
+
+    q[:, 0] = qw
+    q[:, 1] = qx
+    q[:, 2] = qy
+    q[:, 3] = qz
+
+    return q
+
+
+# =======================================================================================================================
+# # axis_angle to ...
+# =======================================================================================================================
+
+
+def axis_angle2axis(axis_angle):
+    x, y, z, angle = axis_angle
+    return np.array([x, y, z]) * angle
+
+
+def axis_angle2quaternion(axis_angle, deg=False):
+    x, y, z, angle = axis_angle
+    if deg:
+        angle = np.deg2rad(angle)
+    s = np.sin(angle / 2)
+    x = x * s
+    y = y * s
+    z = z * s
+    w = np.cos(angle / 2)
+    return w, x, y, z
+
+
+# =======================================================================================================================
+# # quaternion to ...
+# =======================================================================================================================
+
+
+def quaternion2axis_angle(quaternion, deg=False):
+    qw, qx, qy, qz = quaternion / norm(quaternion)
+    angle = 2 * np.arccos(qw)
+    if deg:
+        angle = np.rad2deg(angle)
+    t = np.sqrt(1 - qw**2)
+    x = qx / t
+    y = qy / t
+    z = qz / t
+    return x, y, z, angle
+
+
+def quaternion2matrix(quaternion):
+    q = quaternion / np.linalg.norm(quaternion, axis=0)
+    qw, qx, qy, qz = q
+    sqw, sqx, sqy, sqz = q**2
+
+    qxw = qx * qw
+    qxy = qx * qy
+    qxz = qx * qz
+    qyw = qy * qw
+    qyz = qy * qz
+    qzw = qz * qw
+
+    return np.array(
+        [
+            [sqx - sqy - sqz + sqw, 2.0 * (qxy - qzw), 2.0 * (qxz + qyw)],
+            [2.0 * (qxy + qzw), -sqx + sqy - sqz + sqw, 2.0 * (qyz - qxw)],
+            [2.0 * (qxz - qyw), 2.0 * (qyz + qxw), -sqx - sqy + sqz + sqw],
+        ]
+    )
+
+
+# =======================================================================================================================
+# Matrix to ...
+# =======================================================================================================================
+
+
+def matrix2quaternion(matrix):
+    # method from http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm
+    sqrt = np.sqrt
+    (m00, m01, m02), (m10, m11, m12), (m20, m21, m22) = matrix
+    tr = m00 + m11 + m22
+
+    if tr > 0:
+        S = sqrt(tr + 1.0) * 2  # // S=4*qw
+        qw = 0.25 * S
+        qx = (m21 - m12) / S
+        qy = (m02 - m20) / S
+        qz = (m10 - m01) / S
+    elif (m00 > m11) and (m00 > m22):
+        S = sqrt(1.0 + m00 - m11 - m22) * 2  # // S=4*qx
+        qw = (m21 - m12) / S
+        qx = 0.25 * S
+        qy = (m01 + m10) / S
+        qz = (m02 + m20) / S
+    elif m11 > m22:
+        S = sqrt(1.0 + m11 - m00 - m22) * 2  # // S=4*qy
+        qw = (m02 - m20) / S
+        qx = (m01 + m10) / S
+        qy = 0.25 * S
+        qz = (m12 + m21) / S
+    else:
+        S = sqrt(1.0 + m22 - m00 - m11) * 2  # // S=4*qz
+        qw = (m10 - m01) / S
+        qx = (m02 + m20) / S
+        qy = (m12 + m21) / S
+        qz = 0.25 * S
+    e = np.array([qw, qx, qy, qz])
+    return e / np.sqrt(np.sum(e**2))
+
+
+def matrix2axis_angle(matrix, deg=False):
+    return quaternion2axis_angle(matrix2quaternion(matrix), deg)
+
+
+def matrix2axis(matrix, deg=False):
+    return axis_angle2axis(matrix2axis_angle(matrix, deg))
+
+
+def rotation_matrix_2d(theta: float, degrees=True) -> npt.NDArray:
+    """Generate 2D rotation matrix, rotating with angle theta. Theta is by default degrees, but is calculated in radians if keyword degrees=False.
+
+    Parameters
+    ----------
+    theta : float
+        Angle of the resulting rotation matrix
+    degrees : bool, optional
+        Switch to set the unit of theta. If False, unit is radians. By default True
+
+    Returns
+    -------
+    npt.NDArray
+        The 2D rotation matrix for the angle theta
+    """
+    if degrees:
+        theta = np.deg2rad(theta)
+    R = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
+    return R
+
+
+def projection_2d(theta: float, degrees=True) -> npt.NDArray:
+    """Returns the rotation vector which will project a 2D vector onto the angle theta. Theta is by default degrees, but is calculated in radians if keyword degrees=False.
+
+    Parameters
+    ----------
+    theta : float
+        Angle of the resulting projection vector
+    degrees : bool, optional
+        Switch to set the units of theta. If False, unit is radians. By default True
+
+    Returns
+    -------
+    npt.NDArray
+        The projection vector of angle theta.
+    """
+    if degrees:
+        theta = np.deg2rad(theta)
+    R = np.array([np.cos(theta), np.sin(theta)])
+    return R

wetb/utils/test_envelope.py

+#%%
+
+from wetb.utils.envelope import projected_extremes, compute_ensemble_2d_envelope
+import numpy as np
+import numpy.testing as npt
+
+
+def test_projected_extremes_basic():
+    x = np.cos(np.deg2rad(np.linspace(0,360,12, endpoint=False)))
+    y = np.sin(np.deg2rad(np.linspace(0,360,12, endpoint=False)))
+    
+    res = projected_extremes(
+        np.vstack([x,y]).T,
+        angles= np.linspace(-150, 180, 12),
+        degrees=True
+    )
+    
+    npt.assert_allclose(res[:,1],1.0)
+    
+    
+    
+def test_projected_extremes_radians():
+    x = np.cos(np.linspace(np.pi/4,9*np.pi/4,4, endpoint=False))
+    y = np.sin(np.linspace(np.pi/4,9*np.pi/4,4, endpoint=False))
+
+    res = projected_extremes(
+        signal=np.vstack([x,y]).T,
+        angles=np.linspace(0,2*np.pi,4,endpoint=False),
+        degrees=False
+    )
+    print(res)
+
+    npt.assert_allclose(res[:,1],np.sqrt(2)/2)
+
+
+def test_projected_extremes_sectors():
+    x = np.cos(np.linspace(0,360,12, endpoint=False))
+    y = np.sin(np.linspace(0,360,12, endpoint=False))
+
+    x[x>0] = 0
+
+    res = projected_extremes(
+        signal=np.vstack([x,y]).T,
+        degrees=True,
+        sweep_angle=30
+    )
+
+    npt.assert_allclose(res[(res[:,0]> - 75) * (res[:,0] < 75) ][:,1],0.0)
+
+
+
+# %%

wetb/utils/test_files.py

@@ -6,7 +6,11 @@ Created on 20. jul. 2017
 import os
 import wetb
 import inspect
-wetb_rep_path = os.path.join(os.path.dirname(wetb.__file__), "../")                                   
+from urllib.request import urlretrieve
+
+wetb_rep_path = os.path.abspath(os.path.dirname(wetb.__file__) + "/../").replace("\\", "/") + "/"
+local_TestFiles_path = wetb_rep_path + "TestFiles/"
+remote_TestFiles_url = "https://gitlab.windenergy.dtu.dk/toolbox/TestFiles/raw/master/"
 
 
 def _absolute_filename(filename):
@@ -15,16 +19,17 @@ def _absolute_filename(filename):
         caller_module_path = os.path.dirname(inspect.stack()[index][1])
         tfp = caller_module_path + "/test_files/"
         filename = tfp + filename
-    return filename
+    return os.path.abspath(filename).replace("\\", "/")
 
-def get_test_file(filename):
-    filename = _absolute_filename(filename) 
-    if os.path.exists(filename):
-        return filename
-    else:
-        return os.path.join(wetb_rep_path, 'TestFiles', os.path.relpath(filename, wetb_rep_path))
-        
 
+def get_test_file(filename):
+    filename = _absolute_filename(filename)
+    if not os.path.exists(filename):
+        rel_path = os.path.relpath(filename, wetb_rep_path).replace("\\", "/")
+        filename = local_TestFiles_path + rel_path
+        if not os.path.exists(filename):
+            urlretrieve(remote_TestFiles_url + rel_path, filename)
+    return filename
 
 
 def move2test_files(filename):
@@ -35,8 +40,3 @@ def move2test_files(filename):
     if not os.path.exists(folder):
         os.makedirs(folder)
     os.rename(filename, dst_filename)
-    
-    
-    
-    
-    
\ No newline at end of file

wetb/utils/tests/test_caching.py

@@ -3,14 +3,8 @@ Created on 08/11/2013
 
 @author: mmpe
 '''
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import division
-from __future__ import absolute_import
-from future import standard_library
 from collections import OrderedDict
 import os
-standard_library.install_aliases()
 import multiprocessing
 import time
 import unittest

wetb/utils/tests/test_geometry.py

@@ -3,12 +3,6 @@ Created on 15/01/2014
 
 @author: MMPE
 '''
-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import absolute_import
-from future import standard_library
-standard_library.install_aliases()
 import unittest
 
 

wetb/utils/tests/test_installer.py

+import platform
+import pathlib
+import os
+
+
+import pytest
+
+from wetb.utils.installer import install_wind_tool, install_hawc2_dtu_license, install_keygen_license, install_hawcstab2_dtu_license, install_ellipsys_dtu_license
+import shutil
+
+DESTINATION="/tmp/hawc2"
+
+USER_PLATFORM = platform.uname().system
+if USER_PLATFORM == "Windows":
+    APPDATA = f"{os.environ['APPDATA']}"
+else:
+    APPDATA = "None"
+
+
+TEST_LICENSE_FILE = "/tmp/license.cfg"
+with open(TEST_LICENSE_FILE, "w") as file:
+    file.writelines(["\n[licensing]", "\nhost: www.fakehost.com", "\nport=360"])
+    
+
+def local_license_dir(platform, software):
+    return {
+        "Windows": os.path.join(
+            APPDATA,
+            "DTU Wind Energy",
+            f"{software}",
+        ),
+        "Linux": os.path.join(f"{pathlib.Path.home()}", ".config", f"{software}"),
+        "Darwin": os.path.join(
+            f"{pathlib.Path.home()}", "Library", "Application Support"
+        ),
+    }[platform]
+
+def local_license_file(software):
+    return {
+        "hawc2": "license.cfg",
+        "hawcstab2": "license.cfg",
+        "pywasp": "",
+        "ellipsys": "license.cfg",
+    }[software.lower()]
+
+
+def test_installer_zip():
+    # Install a program distributed as a zip file
+    try:
+        install_wind_tool(
+            tool="HAWC2",
+            destination=DESTINATION
+        )
+    except:
+        raise
+    finally:
+        shutil.rmtree(DESTINATION)
+
+def test_installer_executable():
+    # Install a program distributed as a standalone executable
+    try:
+        install_wind_tool(
+            tool="license_manager",
+            destination=DESTINATION
+        )
+    except:
+        raise
+    finally:
+        shutil.rmtree(DESTINATION)
+
+def test_version_not_available():
+    # Test fallback to latest version
+    try:
+        install_wind_tool(
+            tool="HAWC2",
+            version="1.2.3",
+            destination=DESTINATION
+        )
+    except:
+        raise
+    finally:
+        shutil.rmtree(DESTINATION)
+
+
+@pytest.mark.parametrize("software", ["HAwC2","HAWCStab2", "ellipsys"])
+def test_install_dtu_license(software):
+    software = software.lower()
+    license_path = local_license_dir(USER_PLATFORM, "HAWC2")
+
+    if software == "hawc2":
+        install_hawc2_dtu_license()
+    elif software == "hawcstab2":
+        install_hawcstab2_dtu_license()
+    elif software == "ellipsys":
+        install_ellipsys_dtu_license()
+        
+    assert os.path.exists(f"{local_license_dir(USER_PLATFORM, software)}/{local_license_file(software)}")
+
+    shutil.rmtree(license_path, ignore_errors=True)
+
+
+
+
+TEST_LICENSE_FILE = "/tmp/license.cfg"
+with open(TEST_LICENSE_FILE, "w") as file:
+    file.writelines(["\n[licensing]", "\nhost: www.fakehost.com", "\nport=360"])
+@pytest.mark.parametrize("software,license", [("HAwC2", TEST_LICENSE_FILE),("HAWCStab2", TEST_LICENSE_FILE),("ellipSYS", TEST_LICENSE_FILE)])
+def test_install_keygen_license(software, license):
+    license_path = local_license_dir(USER_PLATFORM, software)
+    try:
+        install_keygen_license(software=software, cfg_file=license)
+        assert os.path.exists(f"{local_license_dir(USER_PLATFORM, software)}/{local_license_file(software)}")
+    except Exception as exc:
+        raise exc
+    finally:
+        shutil.rmtree(license_path, ignore_errors=True)
\ No newline at end of file

wetb/utils/tests/test_pbs_file.py

+from wetb.utils.cluster_tools.pbsfile import PBSFile, Template, PBSMultiRunner
+from wetb.utils.cluster_tools import ssh_client
+from wetb.utils.cluster_tools.ssh_client import SSHClient
+import io
+from wetb.utils.cluster_tools.pbsjob import SSHPBSJob, DONE
+import time
+import pytest
+try:
+    import x
+except ImportError:
+    x = None
+
+
+def test_template():
+    t = Template('[a]B[c]')
+    assert t(a="A", c="C") == "ABC"
+    assert t(a="[c]", c="C") == "CBC", "%s!=%s" % (t(a="[c]", c="C"), 'CBC')
+
+
+def test_pbs_file_str():
+    pbsfile = PBSFile('/home/user/tmp', "test", '''python -c "print('hello world')"''', 'workq')
+    ref = """### Jobid
+#PBS -N test
+### Standard Output
+#PBS -o /home/user/tmp/stdout/test.out
+### merge stderr into stdout
+#PBS -j oe
+#PBS -W umask=0003
+### Maximum wallclock time format HOURS:MINUTES:SECONDS
+#PBS -l walltime=00:10:00
+#PBS -l nodes=1:ppn=1
+### Queue name
+#PBS -q workq
+cd "/home/user/tmp"
+mkdir -p "stdout"
+if [ -z "$PBS_JOBID" ]; then echo "Run using qsub"; exit ; fi
+pwd
+python -c "print('hello world')"
+exit
+"""
+    assert str(pbsfile) == ref
+
+
+def test_pbs_file():
+    if x is None:
+        pytest.xfail("Password missing")
+    pbsfile = PBSFile("/home/mmpe/tmp", "test", '''python -c "print('hello world')"''', 'workq')
+    ssh = SSHClient("jess.dtu.dk", 'mmpe', x.mmpe)
+    pbs_job = SSHPBSJob(ssh)
+    pbs_job.submit(pbsfile, "./tmp")
+    with pbs_job.ssh:
+        start = time.time()
+        while time.time() < start + 10:
+            time.sleep(.1)
+            if pbs_job.status == DONE:
+                break
+        else:
+            raise Exception("job not finished within 10 s")
+        _, out, _ = ssh.execute('cat ./tmp/stdout/test.out')
+    assert "hello world" in out
+
+
+@pytest.mark.parametrize('i,s', [("01:02:03", "01:02:03"),
+                                 (5, "00:00:05"),
+                                 (4000, '01:06:40')])
+def test_pbs_walltime(i, s):
+    pbsfile = PBSFile("./tmp", "test", '', 'workq', walltime=i)
+    assert pbsfile.walltime == s
+
+
+def test_pbs_multirunner():
+    pbs = PBSMultiRunner("/home/user/tmp", )
+    ref = r"""### Jobid
+#PBS -N pbs_multirunner
+### Standard Output
+#PBS -o /home/user/tmp/stdout/pbs_multirunner.out
+### merge stderr into stdout
+#PBS -j oe
+#PBS -W umask=0003
+### Maximum wallclock time format HOURS:MINUTES:SECONDS
+#PBS -l walltime=01:00:00
+#PBS -l nodes=1:ppn=1
+### Queue name
+#PBS -q workq
+cd "/home/user/tmp"
+mkdir -p "stdout"
+if [ -z "$PBS_JOBID" ]; then echo "Run using qsub"; exit ; fi
+pwd
+echo "import os
+import glob
+import numpy as np
+import re
+
+# find available nodes
+with open(os.environ['PBS_NODEFILE']) as fid:
+    nodes = set([f.strip() for f in fid.readlines() if f.strip() != ''])
+pbs_files = [os.path.join(root, f) for root, folders, f_lst in os.walk('.') for f in f_lst if f.endswith('.in')]
+
+# Make a list of [(pbs_in_filename, stdout_filename, walltime),...]
+pat = re.compile(r'[\s\S]*#\s*PBS\s+-o\s+(.*)[\s\S]*(\d\d:\d\d:\d\d)[\s\S]*')
+
+def get_info(f):
+    try:
+        with open(f) as fid:
+            return (f,) + pat.match(fid.read()).groups()
+    except Exception:
+        return (f, f.replace('.in', '.out'), '00:30:00')
+pbs_info_lst = map(get_info, pbs_files)
+
+# sort wrt walltime
+pbs_info_lst = sorted(pbs_info_lst, key=lambda fow: tuple(map(int, fow[2].split(':'))))[::-1]
+# make dict {node1: pbs_info_lst1, ...} and save
+d = dict([(f, pbs_info_lst[i::len(nodes)]) for i, f in enumerate(nodes)])
+with open('pbs.dict', 'w') as fid:
+    fid.write(str(d))
+
+" | python
+
+for node in `cat $PBS_NODEFILE | sort | uniq`
+do
+
+     ssh -T $node << EOF &
+cd "/home/user/tmp"
+python -c "import os
+import multiprocessing
+import platform
+import time
+with open('pbs.dict') as fid:
+    pbs_info_lst = eval(fid.read())[platform.node()]
+arg_lst = ['echo starting %s && mkdir -p "%s" && env PBS_JOBID=$PBS_JOBID "%s" &> "%s" && echo finished %s' %
+           (f, os.path.dirname(o), f, o, f) for f, o, _ in pbs_info_lst]
+print(arg_lst[0])
+print('Starting %d jobs on %s' % (len(arg_lst), platform.node()))
+pool = multiprocessing.Pool(int('$PBS_NUM_PPN'))
+res = pool.map_async(os.system, arg_lst)
+t = time.time()
+for (f, _, _), r in zip(pbs_info_lst, res.get()):
+    print('%-50s\t%s' % (f, ('Errorcode %d' % r, 'Done')[r == 0]))
+print('Done %d jobs on %s in %ds' % (len(arg_lst), platform.node(), time.time() - t))
+
+"
+EOF
+done
+wait
+
+exit
+"""
+    assert str(pbs) == ref

wetb/utils/tests/test_rotation.py

@@ -8,119 +8,188 @@ import unittest
 import numpy as np
 from wetb.utils.geometry import rad
 from wetb.utils.rotation import transformation_matrix, mdot, dots, rotate, rotate_x, \
-    rotmat, rotate_y, rotate_z
+    rotmat, rotate_y, rotate_z, norm, axis2axis_angle, axis_angle2axis, axis2matrix, axis_angle2quaternion,\
+    quaternion2matrix, quaternion2axis_angle, matrix2quaternion, s2matrix
+from tests import npt
+import pytest
 
 x, y, z = 0, 1, 2
-class Test(unittest.TestCase):
 
 
-    def test_transformation_matrix(self):
-        np.testing.assert_array_almost_equal(transformation_matrix(rad(30), x), [[[ 1., 0. , 0.],
-                                                                                  [ 0., 0.8660254, 0.5],
-                                                                                  [ 0., -0.5 , 0.8660254]]])
-        np.testing.assert_array_almost_equal(transformation_matrix(rad(30), y), [[[ 0.8660254, 0. , -0.5],
-                                                                     [ 0., 1., 0 ],
-                                                                     [ 0.5, 0        , 0.8660254]]])
-        np.testing.assert_array_almost_equal(transformation_matrix(rad(30), z), [[[ 0.8660254, 0.5 , 0.],
-                                                                                  [ -0.5, 0.8660254, 0],
-                                                                                  [ 0., 0 , 1]]])
+def test_s2matrix():
+    npt.assert_array_equal(s2matrix('x,-y,-z'), np.array([[1, 0, 0], [0, -1, 0], [0, 0, -1]]).T)
+    npt.assert_array_equal(s2matrix('x,-z,y'), np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]).T)
 
-    def test_rotation_matrix(self):
-        np.testing.assert_array_almost_equal(rotmat(rad(30), x), [[[ 1., 0. , 0.],
-                                                                   [ 0., 0.8660254, -0.5],
-                                                                   [ 0., 0.5 , 0.8660254]]])
-        np.testing.assert_array_almost_equal(rotmat(rad(30), y), [[[ 0.8660254, 0. , 0.5],
-                                                                   [ 0., 1., 0 ],
-                                                                   [ -0.5, 0        , 0.8660254]]])
-        np.testing.assert_array_almost_equal(rotmat(rad(30), z), [[[ 0.8660254, -0.5 , 0.],
-                                                                   [ 0.5, 0.8660254, 0],
-                                                                   [ 0., 0 , 1]]])
-    def test_rotation_matrixes(self):
-        np.testing.assert_array_almost_equal(rotmat([rad(30), rad(60)], 0), [[[ 1., 0. , 0.],
-                                                                     [ 0., 0.8660254, -0.5],
-                                                                     [ 0., 0.5 , 0.8660254]],
-                                                                             [[ 1., 0. , 0.],
-                                                                     [ 0., 0.5, -0.8660254, ],
-                                                                     [ 0., +0.8660254, 0.5 ]]])
 
-    def test_mdot(self):
-        x, y, z = 0, 1, 2
-        m1 = transformation_matrix(np.pi, x)
-        m2 = transformation_matrix([np.pi, np.pi / 2], y)
+def test_transformation_matrix():
+    npt.assert_array_almost_equal(transformation_matrix(rad(30), x), [[[1., 0., 0.],
+                                                                       [0., 0.8660254, 0.5],
+                                                                       [0., -0.5, 0.8660254]]])
+    npt.assert_array_almost_equal(transformation_matrix(rad(30), y), [[[0.8660254, 0., -0.5],
+                                                                       [0., 1., 0],
+                                                                       [0.5, 0, 0.8660254]]])
+    npt.assert_array_almost_equal(transformation_matrix(rad(30), z), [[[0.8660254, 0.5, 0.],
+                                                                       [-0.5, 0.8660254, 0],
+                                                                       [0., 0, 1]]])
 
-        np.testing.assert_array_almost_equal(m1, [[[1, 0, 0], [0, -1, 0], [0, 0, -1]]])
-        np.testing.assert_array_almost_equal(mdot(m1, m1), [[[1, 0, 0], [0, 1, 0], [0, 0, 1]]])
-        np.testing.assert_array_almost_equal(mdot(m1, m2), [[[-1, 0, 0], [0, -1, 0], [0, 0, 1]], [[0, 0, -1], [0, -1, 0], [-1, 0, 0]]])
-        np.testing.assert_array_almost_equal(mdot(m2, m2), [[[1, 0, 0], [0, 1, 0], [0, 0, 1]], [[-1, 0, 0], [0, 1, 0], [0, 0, -1]]])
-        np.testing.assert_array_almost_equal(mdot(m2, m1), [[[-1, 0, 0], [0, -1, 0], [0, 0, 1]], [[0, 0, 1], [0, -1, 0], [1, 0, 0]]])
 
+def test_rotation_matrix():
+    npt.assert_array_almost_equal(rotmat(rad(30), x), [[[1., 0., 0.],
+                                                        [0., 0.8660254, -0.5],
+                                                        [0., 0.5, 0.8660254]]])
+    npt.assert_array_almost_equal(rotmat(rad(30), y), [[[0.8660254, 0., 0.5],
+                                                        [0., 1., 0],
+                                                        [-0.5, 0, 0.8660254]]])
+    npt.assert_array_almost_equal(rotmat(rad(30), z), [[[0.8660254, -0.5, 0.],
+                                                        [0.5, 0.8660254, 0],
+                                                        [0., 0, 1]]])
 
-    def test_dots(self):
-        x, y, z = 0, 1, 2
-        v1 = np.array([1, 2, 3])
-        v2 = np.array([1, 2, 4])
 
-        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi / 2, x), v1).T, [[1, -3, 2]])
-        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi , x), v1).T, [[1, -2, -3]])
+def test_rotation_matrixes():
+    npt.assert_array_almost_equal(rotmat([rad(30), rad(60)], 0), [[[1., 0., 0.],
+                                                                   [0., 0.8660254, -0.5],
+                                                                   [0., 0.5, 0.8660254]],
+                                                                  [[1., 0., 0.],
+                                                                   [0., 0.5, -0.8660254, ],
+                                                                   [0., +0.8660254, 0.5]]])
 
-        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi / 2, y), v1).T, [[3, 2, -1]])
-        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi , y), v1).T, [[-1, 2, -3]])
 
-        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi / 2, z), v1).T, [[-2, 1, 3]])
-        np.testing.assert_array_almost_equal(dots(transformation_matrix(-np.pi , z), v1).T, [[-1, -2, 3]])
+def test_mdot():
+    x, y, z = 0, 1, 2
+    m1 = transformation_matrix(np.pi, x)
+    m2 = transformation_matrix([np.pi, np.pi / 2], y)
 
-        v = np.array([v1, v2]).T
+    npt.assert_array_almost_equal(m1, [[[1, 0, 0], [0, -1, 0], [0, 0, -1]]])
+    npt.assert_array_almost_equal(mdot(m1, m1), [[[1, 0, 0], [0, 1, 0], [0, 0, 1]]])
+    npt.assert_array_almost_equal(
+        mdot(m1, m2), [[[-1, 0, 0], [0, -1, 0], [0, 0, 1]], [[0, 0, -1], [0, -1, 0], [-1, 0, 0]]])
+    npt.assert_array_almost_equal(mdot(m2, m2), [[[1, 0, 0], [0, 1, 0], [0, 0, 1]], [
+                                  [-1, 0, 0], [0, 1, 0], [0, 0, -1]]])
+    npt.assert_array_almost_equal(
+        mdot(m2, m1), [[[-1, 0, 0], [0, -1, 0], [0, 0, 1]], [[0, 0, 1], [0, -1, 0], [1, 0, 0]]])
 
-        np.testing.assert_array_almost_equal(dots(transformation_matrix([-np.pi / 2, np.pi], x), v).T, [[1, -3, 2], [1, -2, -4]])
-        np.testing.assert_array_almost_equal(dots(transformation_matrix([-np.pi / 2, np.pi], y), v).T, [[3, 2, -1], [-1, 2, -4]])
-        np.testing.assert_array_almost_equal(dots(transformation_matrix([-np.pi / 2, np.pi], z), v).T, [[-2, 1, 3], [-1, -2, 4]])
-        np.testing.assert_array_almost_equal(dots([transformation_matrix(np.pi / 2, z), transformation_matrix(np.pi / 2, y), transformation_matrix(np.pi / 2, x)], v1).T, [[3, -2, 1]])
 
+def test_dots():
+    x, y, z = 0, 1, 2
+    v1 = np.array([1, 2, 3])
+    v2 = np.array([1, 2, 4])
 
-    def test_rotate(self):
-        x, y, z = 0, 1, 2
-        v = np.array([1, 2, 3])
+    npt.assert_array_almost_equal(dots(transformation_matrix(-np.pi / 2, x), v1).T, [[1, -3, 2]])
+    npt.assert_array_almost_equal(dots(transformation_matrix(-np.pi, x), v1).T, [[1, -2, -3]])
 
-        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi / 2, x), v), [1, -3, 2])
-        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi, x), v), [1, -2, -3])
+    npt.assert_array_almost_equal(dots(transformation_matrix(-np.pi / 2, y), v1).T, [[3, 2, -1]])
+    npt.assert_array_almost_equal(dots(transformation_matrix(-np.pi, y), v1).T, [[-1, 2, -3]])
 
-        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi / 2, y), v), [3, 2, -1])
-        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi, y), v), [-1, 2, -3])
+    npt.assert_array_almost_equal(dots(transformation_matrix(-np.pi / 2, z), v1).T, [[-2, 1, 3]])
+    npt.assert_array_almost_equal(dots(transformation_matrix(-np.pi, z), v1).T, [[-1, -2, 3]])
 
-        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi / 2, z), v), [-2, 1, 3])
-        np.testing.assert_array_almost_equal(rotate(rotmat(np.pi, z), v), [-1, -2, 3])
+    v = np.array([v1, v2]).T
 
-        v = np.array([[1, 2, 3], [1, 2, 4]])
+    npt.assert_array_almost_equal(dots(transformation_matrix([-np.pi / 2, np.pi], x), v).T, [[1, -3, 2], [1, -2, -4]])
+    npt.assert_array_almost_equal(dots(transformation_matrix([-np.pi / 2, np.pi], y), v).T, [[3, 2, -1], [-1, 2, -4]])
+    npt.assert_array_almost_equal(dots(transformation_matrix([-np.pi / 2, np.pi], z), v).T, [[-2, 1, 3], [-1, -2, 4]])
+    npt.assert_array_almost_equal(dots([transformation_matrix(
+        np.pi / 2, z), transformation_matrix(np.pi / 2, y), transformation_matrix(np.pi / 2, x)], v1).T, [[3, -2, 1]])
 
-        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, -np.pi], x), v)[0], [1, -3, 2])
-        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, -np.pi], x), v)[1], [1, -2, -4])
 
-        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], y), v)[0], [3, 2, -1])
-        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], y), v)[1], [-1, 2, -4])
+def test_rotate():
+    x, y, z = 0, 1, 2
+    v = np.array([1, 2, 3])
 
-        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], z), v)[0], [-2, 1, 3])
-        np.testing.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], z), v)[1], [-1, -2, 4])
+    npt.assert_array_almost_equal(rotate(rotmat(np.pi / 2, x), v), [1, -3, 2])
+    npt.assert_array_almost_equal(rotate(rotmat(np.pi, x), v), [1, -2, -3])
 
+    npt.assert_array_almost_equal(rotate(rotmat(np.pi / 2, y), v), [3, 2, -1])
+    npt.assert_array_almost_equal(rotate(rotmat(np.pi, y), v), [-1, 2, -3])
 
-    def test_rotate_xyz(self):
-        x, y, z = 0, 1, 2
-        v = np.array([1, 2, 3])
+    npt.assert_array_almost_equal(rotate(rotmat(np.pi / 2, z), v), [-2, 1, 3])
+    npt.assert_array_almost_equal(rotate(rotmat(np.pi, z), v), [-1, -2, 3])
 
-        np.testing.assert_array_almost_equal(rotate_x(v, np.pi / 2), [1, -3, 2])
-        np.testing.assert_array_almost_equal(rotate_x(v, -np.pi), [1, -2, -3])
+    v = np.array([[1, 2, 3], [1, 2, 4]])
 
-        np.testing.assert_array_almost_equal(rotate_y(v, np.pi / 2), [3, 2, -1])
-        np.testing.assert_array_almost_equal(rotate_y(v, np.pi), [-1, 2, -3])
+    npt.assert_array_almost_equal(rotate(rotmat([np.pi / 2, -np.pi], x), v)[0], [1, -3, 2])
+    npt.assert_array_almost_equal(rotate(rotmat([np.pi / 2, -np.pi], x), v)[1], [1, -2, -4])
 
-        np.testing.assert_array_almost_equal(rotate_z(v, np.pi / 2), [-2, 1, 3])
-        np.testing.assert_array_almost_equal(rotate_z(v, np.pi), [-1, -2, 3])
+    npt.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], y), v)[0], [3, 2, -1])
+    npt.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], y), v)[1], [-1, 2, -4])
 
-        v = np.array([[1, 2, 3], [4, 5, 6]])
-        np.testing.assert_array_almost_equal(rotate_z(v, np.pi / 2), [[-2, 1, 3], [-5, 4, 6]])
-        np.testing.assert_array_almost_equal(rotate_z(v, [np.pi / 2]), [[-2, 1, 3], [-5, 4, 6]])
-        np.testing.assert_array_almost_equal(rotate_z(v, [-np.pi / 2, np.pi]), [[2, -1, 3], [-4, -5, 6]])
+    npt.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], z), v)[0], [-2, 1, 3])
+    npt.assert_array_almost_equal(rotate(rotmat([np.pi / 2, np.pi], z), v)[1], [-1, -2, 4])
 
 
-if __name__ == "__main__":
-    #import sys;sys.argv = ['', 'Test.test_rad']
-    unittest.main()
+def test_rotate_xyz():
+    x, y, z = 0, 1, 2
+    v = np.array([1, 2, 3])
+
+    npt.assert_array_almost_equal(rotate_x(v, np.pi / 2), [1, -3, 2])
+    npt.assert_array_almost_equal(rotate_x(v, -np.pi), [1, -2, -3])
+
+    npt.assert_array_almost_equal(rotate_y(v, np.pi / 2), [3, 2, -1])
+    npt.assert_array_almost_equal(rotate_y(v, np.pi), [-1, 2, -3])
+
+    npt.assert_array_almost_equal(rotate_z(v, np.pi / 2), [-2, 1, 3])
+    npt.assert_array_almost_equal(rotate_z(v, np.pi), [-1, -2, 3])
+
+    v = np.array([[1, 2, 3], [4, 5, 6]])
+    npt.assert_array_almost_equal(rotate_z(v, np.pi / 2), [[-2, 1, 3], [-5, 4, 6]])
+    npt.assert_array_almost_equal(rotate_z(v, [np.pi / 2]), [[-2, 1, 3], [-5, 4, 6]])
+    npt.assert_array_almost_equal(rotate_z(v, [-np.pi / 2, np.pi]), [[2, -1, 3], [-4, -5, 6]])
+
+
+def test_norm():
+    npt.assert_equal(norm([3, 4]), 5)
+
+
+def test_axis2axis_angle():
+    axis_angle = np.array([2.504687681842697E-002, 0.654193239950699, 0.755912599950848, np.rad2deg(3.09150937138433)])
+    axis = np.array([4.43656429408, 115.877535983, 133.895130906])
+    npt.assert_array_almost_equal(axis2axis_angle(axis), axis_angle)
+    npt.assert_array_almost_equal(axis_angle2axis(axis2axis_angle(axis)), axis)
+
+
+axis_quaternion_matrix_lst = [([30, 0, 0], [0.9659258262890683, 0.25881904510252074, 0.0, 0.0], [[1., 0., 0.],
+                                                                                                 [0., 0.8660254, -0.5],
+                                                                                                 [0., 0.5, 0.8660254]]),
+                              ([0, 30, 0], [0.9659258262890683, 0.0, 0.25881904510252074, 0.0], [[0.8660254, 0., 0.5],
+                                                                                                 [0., 1, 0],
+                                                                                                 [-0.5, 0, 0.8660254]]),
+                              ([0, 0, 30], [0.9659258262890683, 0, 0, 0.25881904510252074], [[0.8660254, -0.5, 0.],
+                                                                                             [0.5, 0.8660254, 0],
+                                                                                             [0., 0, 1]])
+                              ]
+
+
+@pytest.mark.parametrize("axis,quaternion,matrix", axis_quaternion_matrix_lst)
+def test_axis2matrix(axis, quaternion, matrix):
+    npt.assert_array_almost_equal(axis2matrix(axis, deg=True), matrix)
+
+
+def test_axis_angle2axis():
+    axis_angle = np.array([2.504687681842697E-002, 0.654193239950699, 0.755912599950848, np.rad2deg(3.09150937138433)])
+    axis = np.array([4.43656429408, 115.877535983, 133.895130906])
+    npt.assert_array_almost_equal(axis_angle2axis(axis_angle), axis)
+    npt.assert_array_almost_equal(axis2axis_angle(axis_angle2axis(axis_angle)), axis_angle)
+
+
+def test_axis_angle2quaternion():
+    quaternion = np.array([2.518545283038265e-002, 2.518545283038264E-002, 0.657812672860355, 0.760094812138323])
+    axis_angle = np.array([2.504687681842697E-002, 0.654193239950699, 0.755912599950848, 3.09150937138433])
+    npt.assert_array_almost_equal(axis_angle2quaternion(axis_angle), quaternion / norm(quaternion))
+
+
+def test_quaternion2axis_angle():
+    axis_angle = np.array([2.504687681842697E-002, 0.654193239950699, 0.755912599950848, 3.09150937138433])
+    quaternion = np.array([2.518545283038265e-002, 2.518545283038264E-002, 0.657812672860355, 0.760094812138323])
+    npt.assert_array_almost_equal(quaternion2axis_angle(quaternion), axis_angle)
+    npt.assert_array_almost_equal(axis_angle2quaternion(
+        quaternion2axis_angle(quaternion)), quaternion / norm(quaternion))
+
+
+@pytest.mark.parametrize("axis,quaternion,matrix", axis_quaternion_matrix_lst)
+def test_quaternion2matrix(axis, quaternion, matrix):
+    npt.assert_array_almost_equal(quaternion2matrix(quaternion), matrix)
+
+
+@pytest.mark.parametrize("axis,quaternion,matrix", axis_quaternion_matrix_lst)
+def test_matrix2quaternion(axis, quaternion, matrix):
+    npt.assert_array_almost_equal(matrix2quaternion(matrix), quaternion)

wetb/utils/tests/test_test_files.py

@@ -4,7 +4,8 @@ Created on 20. jul. 2017
 @author: mmpe
 '''
 import unittest
-from wetb.utils.test_files import move2test_files, get_test_file
+from wetb.utils.test_files import move2test_files, get_test_file,\
+    local_TestFiles_path
 import os
 from wetb.utils import test_files
 import wetb
@@ -24,6 +25,9 @@ class Test_test_files(unittest.TestCase):
         self.assertTrue(os.path.isfile(dst))
 
     def test_test_files(self):
+        fn = local_TestFiles_path + "wetb/utils/tests/test_files/test_file.txt"
+        if os.path.isfile(fn):
+            os.remove(fn)
         fn1 = get_test_file(tfp+'test_file.txt')
         self.assertTrue(os.path.isfile(fn1))
         fn2 = get_test_file('test_file.txt')

wetb/utils/timing.py

-from __future__ import print_function
-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import absolute_import
-from builtins import dict
-from builtins import zip
-from future import standard_library
-standard_library.install_aliases()
-from six import exec_
 import time
 import inspect
 def get_time(f):

wetb/wind/shear.py

@@ -3,14 +3,8 @@ Created on 16/06/2014
 
 @author: MMPE
 '''
-from __future__ import division
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import absolute_import
-from future import standard_library
-standard_library.install_aliases()
-
-from scipy.optimize.optimize import fmin
+
+from scipy.optimize import fmin
 import numpy as np
 
 
@@ -19,6 +13,7 @@ def _z_u(z_u_lst):
     u = np.array([np.mean(np.array([u])[:]) for _, u in z_u_lst])
     return z, u
 
+
 def power_shear(alpha, z_ref, u_ref):
     """Power shear
 
@@ -41,7 +36,7 @@ def power_shear(alpha, z_ref, u_ref):
     >>> power_shear(.5, 70, 9)([20,50,70])
     [ 4.81070235  7.60638829  9.        ]
     """
-    return lambda z : u_ref * (np.array(z) / z_ref) ** alpha
+    return lambda z: u_ref * (np.array(z) / z_ref) ** alpha
 
 
 def fit_power_shear(z_u_lst):
@@ -70,6 +65,7 @@ def fit_power_shear(z_u_lst):
     alpha, _ = np.polyfit(np.log(z / z_hub), np.log(u / u_hub), 1)
     return alpha
 
+
 def fit_power_shear_ref(z_u_lst, z_ref, plt=None):
     """Estimate power shear parameter, alpha, from two or more specific reference heights using polynomial fit.
 
@@ -108,12 +104,11 @@ def fit_power_shear_ref(z_u_lst, z_ref, plt=None):
         z = np.linspace(min(z), max(z), 100)
         plt.plot(power_shear(alpha, z_ref, u_ref)(z), z)
         plt.margins(.1)
-    if alpha==.1 and u_ref==10: # Initial conditions
+    if alpha == .1 and u_ref == 10:  # Initial conditions
         return np.nan, np.nan
     return alpha, u_ref
 
 
-
 def log_shear(u_star, z0):
     """logarithmic shear
 
@@ -139,21 +134,25 @@ def log_shear(u_star, z0):
     [ 1.73286795  4.02359478  4.86477537]
     """
     K = 0.4  # von Karmans constant
-    def log_shear(z,Obukhov_length=None):
+
+    def log_shear(z, Obukhov_length=None):
         if Obukhov_length is None:
             return u_star / K * (np.log(np.array(z) / z0))
         else:
             return u_star / K * (np.log(z / z0) - stability_term(z, z0, Obukhov_length))
     return log_shear
 
+
 def stability_term(z, z0, L):
     """Calculate the stability term for the log shear
 
     Not validated!!!
     """
     zL = z / L
-    phi_us = lambda zL : (1 + 16 * np.abs(zL)) ** (-1 / 4)  #unstable
-    phi_s = lambda zL : 1 + 5 * zL  # stable
+
+    def phi_us(zL): return (1 + 16 * np.abs(zL)) ** (-1 / 4)  # unstable
+
+    def phi_s(zL): return 1 + 5 * zL  # stable
     phi = np.zeros_like(zL) + np.nan
 
     for m, f in [(((-2 <= zL) & (zL <= 0)), phi_us), (((0 < zL) & (zL <= 1)), phi_s)]:
@@ -199,6 +198,7 @@ def fit_log_shear(z_u_lst, include_R=False):
         return a * kappa, np.exp(-b / a), sum((U - (a * np.log(z) + b)) ** 2)
     return a * kappa, np.exp(-b / a)
 
+
 if __name__ == '__main__':
     from matplotlib.pyplot import plot, show
     z = np.arange(0, 211)

wetb/wind/stability.py

@@ -5,7 +5,6 @@ Created on 05/10/2015
 '''
 
 import numpy as np
-from scipy.integrate import trapz
 from scipy.signal.signaltools import detrend
 
 
@@ -38,7 +37,7 @@ def MoninObukhov_length(u,v,w, T):
 
 def L2category(L, full_category_name=False):
     """Stability category from Monin-Obukhov length
-    
+
     Categories:
     0>L>-50: Extreme unstable (eu)
     -50>L>-100: Very unstable (vu)
@@ -57,15 +56,15 @@ def L2category(L, full_category_name=False):
     full_category_name : bool, optional
         If False, default, category ids are returned, e.g. "n" for neutral
         If True, full name of category are returned, e.g. "Neutral"
-    
+
     Returns
     -------
     Stability category : str
-    
+
         Examples
     --------
     >>> L2category(1000)
-    n 
+    n
     """
     cat_limits = np.array([-1e-99,-50,-100,-200,-500,500,200,50,10,1e-99])
     index = np.searchsorted( 1/cat_limits, 1/np.array(L))-1
@@ -73,7 +72,7 @@ def L2category(L, full_category_name=False):
         return np.array(['Extreme unstable', 'Very unstable','Unstable','Near unstable','Neutral','Near stable','Stable','Very stable','Extreme stable','Undefined'])[index]
     else:
         return np.array(['eu', 'vu','u','nu','n','ns','s','vs','es','-'])[index]
-    
+
 def MoninObukhov_length2(u_star, w, T, specific_humidity=None):
     """Calculate the Monin Obukhov length
 

wetb/wind/tests/test_Shear.py

@@ -3,13 +3,6 @@ Created on 05/06/2012
 
 @author: Mads
 '''
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import division
-from __future__ import absolute_import
-from builtins import zip
-from future import standard_library
-standard_library.install_aliases()
 import os
 from wetb.wind.utils import xyz2uvw
 import wetb.gtsdf

wetb/wind/tests/test_weibull.py

@@ -8,8 +8,8 @@ import matplotlib.pyplot as plt
 import numpy as np
 from wetb.wind import weibull
 
-class TestWeibull(unittest.TestCase):
 
+class TestWeibull(unittest.TestCase):
 
     def test_weibull(self):
         wb = weibull.pdf(4, 2)
@@ -21,37 +21,39 @@ class TestWeibull(unittest.TestCase):
 
     def test_random_weibull(self):
         y = weibull.random(4, 2, 1000000)
-        pdf, x = np.histogram(y, 100, normed=True)
+        pdf, x = np.histogram(y, 100, density=True)
         x = (x[1:] + x[:-1]) / 2
         self.assertLess(sum((pdf - weibull.pdf(4, 2)(x)) ** 2), 0.0001)
         if 0:
             plt.plot(x, pdf)
             plt.plot(x, weibull.pdf(4, 2)(x))
             plt.show()
- 
+
     def test_fit_weibull(self):
+        np.random.seed(1)
         y = weibull.random(4, 2, 1000000)
         A, k = weibull.fit(y)
         self.assertAlmostEqual(A, 4, delta=0.01)
         self.assertAlmostEqual(k, 2, delta=0.01)
         if 0:
-            plt.hist(y, 100, normed=True)
+            plt.hist(y, 100, density=True)
             x = np.arange(0, 20, .1)
             plt.plot(x, weibull.pdf(4, 2)(x))
             plt.show()
-             
+
     def test_weibull_cdf(self):
         wb = weibull.pdf(4, 2)
         x = np.arange(0, 20, .01)
-        cdf = weibull.cdf(4,2)
-        self.assertEqual(cdf(999),1)
-        np.testing.assert_array_almost_equal(np.cumsum(wb(x))/100, cdf(x), 3)
-         
+        cdf = weibull.cdf(4, 2)
+        self.assertEqual(cdf(999), 1)
+        np.testing.assert_array_almost_equal(np.cumsum(wb(x)) / 100, cdf(x), 3)
+
         if 0:
-            plt.plot(x, np.cumsum(wb(x))*.01)
+            plt.plot(x, np.cumsum(wb(x)) * .01)
             plt.plot(x, cdf(x))
             plt.show()
-             
+
+
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.testweibull']
     unittest.main()