diff --git a/wind_tools/VTKbased/cutoffVisualizationTool.py b/wind_tools/VTKbased/cutoffVisualizationTool.py new file mode 100644 index 0000000..378efd8 --- /dev/null +++ b/wind_tools/VTKbased/cutoffVisualizationTool.py @@ -0,0 +1,152 @@ +#!/usr/bin/python3 +# -*- coding: utf-8 -*- + +import vtk +from PyQt5.QtWidgets import QWidget, QSlider, QLabel, QVBoxLayout +from PyQt5.QtCore import Qt +import numpy as np +from vtk.qt.QVTKRenderWindowInteractor import QVTKRenderWindowInteractor + + +class cutoffInterface(QWidget): + def __init__(self, fileLoc): + super().__init__() + self.OpacB = 3 + self.OpacV = .1 + + # Start the VTK viewer and the user interface + self.ren, self.oTF, self.scalar_range, self.scalar_bar = instantiateVTKviewer(fileLoc) + + self.initUI() + self.rW = self.vtkWidget.GetRenderWindow() + self.iren = self.rW.GetInteractor() + self.rW.AddRenderer(self.ren) + + # Change the VTK element interaction style and show the GUI + self.iren.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera()) + self.iren.Initialize() + self.show() + + def initUI(self): + layout = QVBoxLayout() + self.setLayout(layout) + self.l1 = QLabel('Cutoff value = 3 m/s') + self.l1.setAlignment(Qt.AlignLeft) + layout.addWidget(self.l1) + + sldCuttoff = QSlider(Qt.Horizontal) + sldCuttoff.setSingleStep(.1) + sldCuttoff.setGeometry(30, 40, 100, 30) + sldCuttoff.valueChanged[int].connect(self.changeOpacityBoundary) + layout.addWidget(sldCuttoff) + + self.l2 = QLabel('Transparancy value = 10%') + self.l2.setAlignment(Qt.AlignLeft) + layout.addWidget(self.l2) + + sldTransparancy = QSlider(Qt.Horizontal) + sldTransparancy.setSingleStep(.1) + sldTransparancy.setGeometry(30, 40, 100, 30) + sldTransparancy.valueChanged[int].connect(self.changeOpacity) + layout.addWidget(sldTransparancy) + + self.vtkWidget = QVTKRenderWindowInteractor() + layout.addWidget(self.vtkWidget) + + self.setGeometry(300, 100, 1000, 800) + self.setWindowTitle('Slicer Interface') + + def changeOpacityBoundary(self, value): + self.OpacB = (self.scalar_range[0] + + (self.scalar_range[1]-self.scalar_range[0])*value/100) + self.l1.setText('Cutoff value = %.3f m/s' % self.OpacB) + self.updateLut() + + def changeOpacity(self, value): + self.OpacV = value/1000 + self.l2.setText('Transparancy value = %.1f%%' % (value/10)) + self.updateLut() + + def updateLut(self): + self.oTF.RemoveAllPoints() + self.oTF.AddPoint(self.scalar_range[0], self.OpacV) + self.oTF.AddPoint(self.OpacB-.01, self.OpacV) + self.oTF.AddPoint(self.OpacB+.01, 0) + self.oTF.AddPoint(self.scalar_range[1], 0) + self.rW.Render() + + +def instantiateVTKviewer(fileLoc): + # Create the reader for the data + reader = vtk.vtkXMLImageDataReader() + reader.SetFileName(fileLoc) + reader.Update() + flowField = reader.GetOutput() + scalar_range = flowField.GetScalarRange() + + # Create a custom lut, it's used both as a mapper and at the scalar_bar + lut = vtk.vtkLookupTable() + lut.SetTableRange(scalar_range) + lut.Build() + + # create the scalar_bar + scalar_bar = vtk.vtkScalarBarActor() + scalar_bar.SetLookupTable(lut) + scalar_bar.SetNumberOfLabels(6) + scalar_bar.GetLabelTextProperty().SetFontFamilyToCourier() + scalar_bar.GetLabelTextProperty().SetJustificationToRight() + scalar_bar.GetLabelTextProperty().SetVerticalJustificationToCentered() + scalar_bar.GetLabelTextProperty().BoldOff() + scalar_bar.GetLabelTextProperty().ItalicOff() + scalar_bar.GetLabelTextProperty().ShadowOff() + scalar_bar.GetLabelTextProperty().SetColor(0, 0, 0) + + # Setup rendering + # Create transfer mapping scalar value to opacity + opacityTransferFunction = vtk.vtkPiecewiseFunction() + opacityTransferFunction.AddPoint(0, .1) + opacityTransferFunction.AddPoint(2.95, .1) + opacityTransferFunction.AddPoint(3.05, 0) + opacityTransferFunction.AddPoint(7, 0) + + # Create transfer mapping scalar value to color + colorTransferFunction = vtk.vtkColorTransferFunction() + for s in np.linspace(scalar_range[0], scalar_range[1], 200): + col = [0, 0, 0] + lut.GetColor(s, col) + colorTransferFunction.AddRGBPoint(s, col[0], col[1], col[2]) + + # The property describes how the data will look + volumeProperty = vtk.vtkVolumeProperty() + volumeProperty.SetColor(colorTransferFunction) + volumeProperty.SetScalarOpacity(opacityTransferFunction) + volumeProperty.SetInterpolationTypeToLinear() + + volumeProperty.SetComponentWeight(0, 1.0) + volumeProperty.SetComponentWeight(1, 0.0) + volumeProperty.SetComponentWeight(2, 0.0) + + # The mapper / ray cast function know how to render the data + volumeMapper = vtk.vtkGPUVolumeRayCastMapper() + volumeMapper.SetInputConnection(reader.GetOutputPort()) + + # The volume holds the mapper and the property and + # can be used to position/orient the volume + volume = vtk.vtkVolume() + volume.SetMapper(volumeMapper) + volume.SetProperty(volumeProperty) + + # Setup camera + camera = vtk.vtkCamera() + camera.SetPosition(-800, -400, 300) + camera.SetFocalPoint(0, 0, 0) + camera.SetViewUp(0, 0, 1) + + renderer = vtk.vtkRenderer() + renderer.AddVolume(volume) + renderer.AddActor(scalar_bar) + renderer.SetBackground(240/255, 240/255, 240/255) + renderer.SetActiveCamera(camera) + renderer.ResetCamera() + + return renderer, opacityTransferFunction, scalar_range, scalar_bar diff --git a/wind_tools/VTKbased/flowfieldWriter.py b/wind_tools/VTKbased/flowfieldWriter.py new file mode 100644 index 0000000..e459121 --- /dev/null +++ b/wind_tools/VTKbased/flowfieldWriter.py @@ -0,0 +1,221 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Feb 12 16:07:47 2018 + +@author: roald +""" + + +import vtk +import numpy as np +from pandas import DataFrame + + +def create_vti_from_df(fileName, df, spacing, dimensions, origin): + """Create a .vti fileformat for use with the vtk toolbox + + input: + fileName: The fileName of the .vti file to be created + df: a pandas table with the columns, x,y,z,u,v,w of all relevant flow info + spacing: The spacing in the x, y and z direction between points + dimensions: a tuple or list of integers with the number of points in xyz + origin: the origin of the flow field, for reconstructing turbine coords + + output: - + + Roald Storm, 2018 """ + imageData = vtk.vtkImageData() + imageData.SetDimensions(list(int(x) for x in dimensions)) + imageData.AllocateScalars(vtk.VTK_DOUBLE, 3) + imageData.SetSpacing(spacing) + imageData.SetOrigin(origin) + + # Fill every entry of the image data with "2.0" + i = 0 + for z in range(int(dimensions[2])): + for y in range(int(dimensions[1])): + for x in range(int(dimensions[0])): + imageData.SetScalarComponentFromDouble( + x, y, z, 0, df.u[i]) + imageData.SetScalarComponentFromDouble( + x, y, z, 1, df.v[i]) + imageData.SetScalarComponentFromDouble( + x, y, z, 2, df.w[i]) + i += 1 + + writer = vtk.vtkXMLImageDataWriter() + writer.SetFileName(fileName) + writer.SetInputData(imageData) + writer.Write() + + +def create_df_from_vti(fileName): + """Read flow array output from SOWFA + + + input: + filename: name of .vti file to open + + output: + df: a pandas table with the columns, x,y,z,u,v,w of all relevant flow info + spacing: The spacing in the x, y and z direction between points + dimensions: a tuple or list of integers with the number of points in xyz + origin: the origin of the flow field, for reconstructing turbine coords + + Roald Storm, 2018 """ + + reader = vtk.vtkXMLImageDataReader() + reader.SetFileName(fileName) + reader.Update() + + readerOutput = reader.GetOutput() + dimensions = readerOutput.GetDimensions() + Nx, Ny, Nz = dimensions + spacing = readerOutput.GetSpacing() + origin = readerOutput.GetOrigin() + + xRange = np.arange(0, Nx*spacing[0], spacing[0]) + yRange = np.arange(0, Ny*spacing[1], spacing[1]) + zRange = np.arange(0, Nz*spacing[2], spacing[2]) + + pts = np.array([(x, y, z) for z in zRange for y in yRange for x in xRange]) + vals = np.zeros(pts.shape) + + for z in range(0, Nz): + for y in range(0, Ny): + for x in range(0, Nx): + vals[z*Ny*Nx + y*Nx + x] = [ + readerOutput.GetScalarComponentAsDouble(x, y, z, 0), + readerOutput.GetScalarComponentAsDouble(x, y, z, 1), + readerOutput.GetScalarComponentAsDouble(x, y, z, 2)] + + df = DataFrame(columns=('x', 'y', 'z', 'u', 'v', 'w')) + + df['u'] = vals[:, 0] + df['v'] = vals[:, 1] + df['w'] = vals[:, 2] + df['x'] = pts[:, 0] + df['y'] = pts[:, 1] + df['z'] = pts[:, 2] + + return df, spacing, dimensions, origin + + +def create_df_from_vtk_poly(fileName): + """Read flow array output from SOWFA + + + input: + filename: name of .vtk file to open with a 2D flowfield slice specified + as PolyData + + output: + df: a pandas table with the columns, x,y,z,u,v,w of all relevant flow info + spacing: The spacing in the x, y and z direction between points + dimensions: a tuple or list of integers with the number of points in xyz + origin: the origin of the flow field, for reconstructing turbine coords + + Roald Storm, 2018 """ + + reader = vtk.vtkGenericDataObjectReader() + reader.SetFileName(fileName) + reader.Update() + readerOutput = reader.GetOutput() + + if not reader.IsFilePolyData(): + raise ValueError('This .vtk file does not have PolyData') + + if readerOutput.GetCell(0).GetPoints().GetData().GetNumberOfTuples() != 4: + raise ValueError('This .vtk file is not a 2d flowfield slice') + + ncells = readerOutput.GetNumberOfCells() + cellData = readerOutput.GetCellData() + + pts = np.zeros([ncells, 3]) + vals = np.zeros([ncells, 3]) + data = cellData.GetArray(0) + + for i in range(data.GetNumberOfTuples()): + vals[i] = data.GetTuple(i) + pts[i] = [min(i) for i in zip(*[readerOutput.GetCell(i).GetPoints().GetData().GetTuple(x) for x in range(4)])] + + df = DataFrame(columns=('x', 'y', 'z', 'u', 'v', 'w')) + + df['u'] = vals[:, 0] + df['v'] = vals[:, 1] + df['w'] = vals[:, 2] + df['x'] = pts[:, 0] + df['y'] = pts[:, 1] + df['z'] = pts[:, 2] + + df = df.sort_values(by=['z', 'y', 'x']).reset_index(drop=True) + + dimensions = (len(df[(df['y'] == 0)].index), + len(df[(df['x'] == 0)].index), + 1) + spacing = (df[df['x'] != 0].iloc[0].loc['x'], + df[df['y'] != 0].iloc[0].loc['y'], + 0) + origin = (0.0, 0.0, 0.0) + + return df, spacing, dimensions, origin + + +def create_df_from_vtk_unstructuredgrid(fileName): + """Read flow array output from SOWFA + + + input: + filename: name of .vtk file to open with 3D flowfield specified as an + unstructured grid + + output: + df: a pandas table with the columns, x,y,z,u,v,w of all relevant flow info + spacing: The spacing in the x, y and z direction between points + dimensions: a tuple or list of integers with the number of points in xyz + origin: the origin of the flow field, for reconstructing turbine coords + + Roald Storm, 2018 """ + + reader = vtk.vtkGenericDataObjectReader() + reader.SetFileName(fileName) + reader.Update() + readerOutput = reader.GetOutput() + + if not reader.IsFileUnstructuredGrid(): + raise ValueError('This .vtk file does not have Unstructured Grid Data') + + if readerOutput.GetCell(0).GetPoints().GetData().GetNumberOfTuples() != 8: + raise ValueError('This .vtk file is not a 3d flowfield slice') + + ncells = readerOutput.GetNumberOfCells() + cellData = readerOutput.GetCellData() + + pts = np.zeros([ncells, 3]) + vals = np.zeros([ncells, 3]) + data = cellData.GetArray(1) + + for i in range(data.GetNumberOfTuples()): + vals[i] = data.GetTuple(i) + pts[i] = [min(i) for i in zip(*[readerOutput.GetCell(i).GetPoints().GetData().GetTuple(x) for x in range(8)])] + + df = DataFrame(columns=('x', 'y', 'z', 'u', 'v', 'w')) + + df['u'] = vals[:, 0] + df['v'] = vals[:, 1] + df['w'] = vals[:, 2] + df['x'] = pts[:, 0] + df['y'] = pts[:, 1] + df['z'] = pts[:, 2] + + df = df.sort_values(by=['z', 'y', 'x']).reset_index(drop=True) + + dimensions = (len(df[(df['y'] == 0) & (df['z'] == 0)].index), + len(df[(df['x'] == 0) & (df['z'] == 0)].index), + len(df[(df['x'] == 0) & (df['y'] == 0)].index)) + spacing = (df[df['x'] != 0].iloc[0].loc['x'], + df[df['y'] != 0].iloc[0].loc['y'], + df[df['z'] != 0].iloc[0].loc['z']) + origin = (0.0, 0.0, 0.0) + + return df, spacing, dimensions, origin diff --git a/wind_tools/VTKbased/slicerTool.py b/wind_tools/VTKbased/slicerTool.py new file mode 100644 index 0000000..2010319 --- /dev/null +++ b/wind_tools/VTKbased/slicerTool.py @@ -0,0 +1,174 @@ +#!/usr/bin/python3 +# -*- coding: utf-8 -*- + +import vtk +from PyQt5.QtWidgets import QWidget, QSlider, QRadioButton, QGridLayout +from PyQt5.QtCore import Qt +from vtk.qt.QVTKRenderWindowInteractor import QVTKRenderWindowInteractor + + +class slicerInterface(QWidget): + def __init__(self, fileLoc): + super().__init__() + + # Start the VTK viewer and the user interface + (self.ren, self.scalar_bar, self.dims, self.cellSizes, self.reader, + self.lut, self.scalar_range) = instantiateVTKviewer2(fileLoc) + + self.initUI() + self.rW = self.vtkWidget.GetRenderWindow() + self.iren = self.rW.GetInteractor() + self.rW.AddRenderer(self.ren) + + # Instantiate the dynamical slice that will be changed on user input + self.dynActor = makePlaneAt((1, 0, 0), (0, 0, 0), + self.reader, self.lut, self.scalar_range) + self.ren.AddActor(self.dynActor) + self.rButtonMap = {'X': 0, 'Y': 1, 'Z': 2} + self.sliceVal = 0 + self.ax = 0 + + # Change the VTK element interaction style and show the GUI + self.iren.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera()) + self.iren.Initialize() + self.show() + + def initUI(self): + Grid = QGridLayout() + self.setLayout(Grid) + + rButtonX = QRadioButton('X') + rButtonX.setChecked(True) + rButtonX.toggled.connect(self.on_radio_button_toggled) + Grid.addWidget(rButtonX, 0, 0) + + rButtonY = QRadioButton('Y') + rButtonY.toggled.connect(self.on_radio_button_toggled) + Grid.addWidget(rButtonY, 0, 1) + + rButtonZ = QRadioButton('Z') + rButtonZ.toggled.connect(self.on_radio_button_toggled) + Grid.addWidget(rButtonZ, 0, 2) + + sld = QSlider(Qt.Horizontal) + sld.setSingleStep(.1) + sld.setGeometry(30, 40, 100, 30) + sld.valueChanged[int].connect(self.changeSliderValue) + Grid.addWidget(sld, 1, 0, 1, 3) + + self.vtkWidget = QVTKRenderWindowInteractor() + Grid.addWidget(self.vtkWidget, 2, 0, 1, 3) + + self.setGeometry(300, 100, 1000, 800) + self.setWindowTitle('Slicer Interface') + + def changeSliderValue(self, value): + self.sliceVal = value/100 + self.changeSlice() + + def on_radio_button_toggled(self): + self.ax = self.rButtonMap[self.sender().text()] + self.changeSlice() + + def changeSlice(self): + self.ren.RemoveActor(self.dynActor) + Nor = [0, 0, 0] + Nor[self.ax] = 1 + Ori = [0, 0, 0] + Ori[self.ax] = self.dims[self.ax]*self.cellSizes[self.ax]*self.sliceVal + self.dynActor = makePlaneAt(Nor, Ori, self.reader, + self.lut, self.scalar_range) + self.ren.AddActor(self.dynActor) + self.rW.Render() + + +def instantiateVTKviewer2(fileLoc): + # Read the source file. + reader = vtk.vtkXMLImageDataReader() + reader.SetFileName(fileLoc) + reader.Update() + flowField = reader.GetOutput() + scalar_range = flowField.GetScalarRange() + + # Create a custom LookUpTable. + # The lut is used both at the mapper and at the scalar_bar + lut = vtk.vtkLookupTable() + lut.SetTableRange(scalar_range) + lut.Build() + + dims = flowField.GetDimensions() + cellSizes = flowField.GetSpacing() + + # create the scalar_bar + scalar_bar = vtk.vtkScalarBarActor() + scalar_bar.SetLookupTable(lut) + scalar_bar.SetNumberOfLabels(8) + scalar_bar.GetLabelTextProperty().SetFontFamilyToCourier() + scalar_bar.GetLabelTextProperty().SetJustificationToRight() + scalar_bar.GetLabelTextProperty().SetVerticalJustificationToCentered() + scalar_bar.GetLabelTextProperty().BoldOff() + scalar_bar.GetLabelTextProperty().ItalicOff() + scalar_bar.GetLabelTextProperty().ShadowOff() + scalar_bar.GetLabelTextProperty().SetColor(0, 0, 0) + + # Create the three planes at the back sides of the volume + origin = flowField.GetOrigin() + box = [o + d*s for o, d, s in zip(origin, dims, cellSizes)] + + cutActor1 = makePlaneAt((1, 0, 0), (box[0]*.99, 0, 0), + reader, lut, scalar_range) + cutActor2 = makePlaneAt((0, 1, 0), (0, box[1]*.99, 0), + reader, lut, scalar_range) + cutActor3 = makePlaneAt((0, 0, 1), (0, 0, box[2]*.01), + reader, lut, scalar_range) + + # Setup camera + camera = vtk.vtkCamera() + camera.SetPosition(-800, -400, 300) + camera.SetFocalPoint(0, 0, 0) + camera.SetViewUp(0, 0, 1) + + # Setup lights, One for each plane to prevent any kind of shadow artefacts + light1 = vtk.vtkLight() + light1.SetPosition(0, 0, 0) + light1.SetFocalPoint(1, 0, 0) + light2 = vtk.vtkLight() + light2.SetPosition(0, 0, 0) + light2.SetFocalPoint(0, 1, 0) + light3 = vtk.vtkLight() + light3.SetPosition(0, 0, 1000) + light3.SetFocalPoint(0, 0, 1) + + # Setup rendering + renderer = vtk.vtkRenderer() + renderer.AddActor(cutActor1) + renderer.AddActor(cutActor2) + renderer.AddActor(cutActor3) + renderer.AddActor(scalar_bar) + renderer.SetBackground(240/255, 240/255, 240/255) + renderer.SetActiveCamera(camera) + renderer.ResetCamera() + renderer.AddLight(light1) + renderer.AddLight(light2) + renderer.AddLight(light3) + + return (renderer, scalar_bar, dims, cellSizes, reader, lut, scalar_range) + + +def makePlaneAt(Nor, Ori, reader, lut, scalar_range): + plane = vtk.vtkPlane() + plane.SetOrigin(Ori) + plane.SetNormal(Nor) + + planeCut = vtk.vtkCutter() + planeCut.SetInputConnection(reader.GetOutputPort()) + planeCut.SetCutFunction(plane) + + cutMapper = vtk.vtkPolyDataMapper() + cutMapper.SetInputConnection(planeCut.GetOutputPort()) + cutMapper.SetScalarRange(scalar_range) + cutMapper.SetLookupTable(lut) + + Actor = vtk.vtkActor() + Actor.SetMapper(cutMapper) + return Actor diff --git a/wind_tools/flow/flow_field.py b/wind_tools/flow/flow_field.py index 02e5fbe..032b759 100644 --- a/wind_tools/flow/flow_field.py +++ b/wind_tools/flow/flow_field.py @@ -16,19 +16,20 @@ import pandas as pd -def get_flow_file(case_folder): +def get_flow_file(case_folder, fileType = ".vtk"): """Given a case folder, find the flow file - input: case_folder: name of case folder + input + case_folder: name of case folder output: - flow_file: full path name of flow file""" + flow_file: full path name of flow file""" array_folder = os.path.join(case_folder,'array.mean') time_folder = os.path.join(array_folder,os.listdir(array_folder)[0]) - flow_file = os.path.join(time_folder,os.listdir(time_folder)[0]) - flow_file = os.path.join(time_folder,os.listdir(time_folder)[0]) + fileNames = [fi for fi in os.listdir(time_folder) if fi.endswith(fileType)] + flow_file = os.path.join(time_folder,fileNames[0]) return flow_file @@ -40,7 +41,9 @@ def read_flow_frame_SOWFA(filename): input: filename: name of flow array to open output: - df: a pandas table with the columns, x,y,z,u,v,w of all relavent flow info + df: a pandas table with the columns, x,y,z,u,v,w of all relevant flow info + spacing: The spacing in the x, y and z direction between points + dimensions: a tuple or list of integers with the number of points in xyz origin: the origin of the flow field, for reconstructing turbine coords Paul Fleming, 2018 """ @@ -53,7 +56,8 @@ def read_flow_frame_SOWFA(filename): if 'SPACING' in read_data: spacing = tuple([float(d) for d in read_data.rstrip().split(' ')[1:]]) if 'DIMENSIONS' in read_data: - dimensions = tuple([float(d) for d in read_data.rstrip().split(' ')[1:]]) + dims = tuple([float(d) for d in read_data.rstrip().split(' ')[1:]]) + dimensions = tuple(int(x) for x in dims) if 'ORIGIN' in read_data: origin = tuple([float(d) for d in read_data.rstrip().split(' ')[1:]]) @@ -75,11 +79,11 @@ def get_flow_frame_FLORIS(floris): """Read flow array output from SOWFA - input: floris: a floris model which has already been run to extract flow from + input: + floris: a floris model which has already been run to extract flow from output: - df: a pandas table with the columns, x,y,z,u,v,w of all relavent flow info - origin: the origin of the flow field, for reconstructing turbine coords + df: a pandas table with the columns, x,y,z,u,v,w of all relevant flow info Paul Fleming, 2018 """