Python vtk.util.numpy_support.numpy_to_vtk() Examples

def create_surface_points(self, vertices):
        """
        Method to create the points that form the surfaces
        Args:
            vertices (numpy.array): 2D array (XYZ) with the coordinates of the points

        Returns:
            vtk.vtkPoints: with the coordinates of the points
        """
        Points = vtk.vtkPoints()
        if self.ve != 1:
            vertices[:, 2] = vertices[:, 2] * self.ve
        #     raise NotImplementedError('Vertical exageration for surfaces not implemented yet.')
        # for v in vertices:
        #     v[-1] = self.ve * v[-1]
        #     Points.InsertNextPoint(v)
        Points.SetData(numpy_to_vtk(vertices))

        return Points 

def __create_structured_grid(ptsMat, dims, models=None):
        """An internal helper to build out structured grids"""
        # Adjust if result was 2D:
        if ptsMat.shape[1] == 2:
            # Figure out which dim is null
            nullDim = dims.index(None)
            ptsMat = np.insert(ptsMat, nullDim, np.zeros(ptsMat.shape[0]), axis=1)
        if ptsMat.shape[1] != 3:
            raise RuntimeError('Points of the mesh are improperly defined.')
        # Convert the points
        vtkPts = _vtk.vtkPoints()
        vtkPts.SetData(_nps.numpy_to_vtk(ptsMat, deep=True))
        # Uncover hidden dimension
        for i, d in enumerate(dims):
            if d is None:
                dims[i] = 0
            dims[i] = dims[i] + 1
        output = _vtk.vtkStructuredGrid()
        output.SetDimensions(dims[0], dims[1], dims[2]) # note this subtracts 1
        output.SetPoints(vtkPts)
        # Assign the model('s) to the object
        return assign_cell_data(output, models=models) 

def CreateTriangles(self, pc, triangles):
        "Create a mesh from triangles"

        self.verts = vtk.vtkPoints()
        self.points_npy = pc[:, :3].copy()
        self.verts.SetData(numpy_support.numpy_to_vtk(self.points_npy))

        nTri = len(triangles)
        self.cells = vtk.vtkCellArray()
        self.pd = vtk.vtkPolyData()
        # - Note that the cell array looks like this: [3 vtx0 vtx1 vtx2 3 vtx3 ... ]
        self.cells_npy = np.column_stack(
            [np.full(nTri, 3, dtype=np.int64), triangles.astype(np.int64)]
        ).ravel()
        self.cells.SetCells(nTri, numpy_support.numpy_to_vtkIdTypeArray(self.cells_npy))

        self.pd.SetPoints(self.verts)
        self.pd.SetPolys(self.cells)

        self.SetupPipelineMesh() 

def volume_to_vtk(volelement):
    """Convert the volume element to a VTK data object.

    Args:
        volelement (:class:`omf.volume.VolumeElement`): The volume element to
            convert

    """
    output = volume_grid_geom_to_vtk(volelement.geometry)
    shp = get_volume_shape(volelement.geometry)
    # Add data to output
    for data in volelement.data:
        arr = data.array.array
        arr = np.reshape(arr, shp).flatten(order='F')
        c = nps.numpy_to_vtk(num_array=arr, deep=True)
        c.SetName(data.name)
        loc = data.location
        if loc == 'vertices':
            output.GetPointData().AddArray(c)
        else:
            output.GetCellData().AddArray(c)
    return pyvista.wrap(output)


# Now set up the display names for the docs 

def assign_cell_data(vtkDS, models=None):
    """Assign the model(s) to the VTK dataset as CellData

    Parameters
    ----------

    vtkDS : pyvista.Common
        Any given VTK data object that has cell data

    models : dict(numpy.ndarray)
        Name('s) and array('s). Match number of cells

    """
    nc = vtkDS.GetNumberOfCells()
    if models is not None:
        for name, mod in models.items():
            # Convert numpy array
            if mod.shape[0] != nc:
                raise RuntimeError('Number of model cells ({}) (first axis of model array) for "{}" does not match number of mesh cells ({}).'.format(mod.shape[0], name, nc))
            vtkDoubleArr = _nps.numpy_to_vtk(mod, deep=1)
            vtkDoubleArr.SetName(name)
            vtkDS.GetCellData().AddArray(vtkDoubleArr)
    return vtkDS 

def convert_cell_conn(cell_connectivity):
    """Converts cell connectivity arrays to a cell matrix array that makes sense
    for VTK cell arrays.
    """
    cellsMat = np.concatenate(
            (
                np.ones((cell_connectivity.shape[0], 1), dtype=np.int64)*cell_connectivity.shape[1],
                cell_connectivity
            ),
            axis=1).ravel()
    return nps.numpy_to_vtk(cellsMat, deep=True, array_type=vtk.VTK_ID_TYPE) 

def AddColors(self, cdata):
        """"Add colors (Nx3 NumPy array) to the current point cloud visualization object        
        NumPy array should be of type uint8 with R, G and B values between 0 and 255
        """

        nColors = cdata.shape[0]
        nDim = cdata.shape[1]

        if self.pd.GetNumberOfPoints() == 0:
            raise Exception("No points to add color for")

        if nColors != self.pd.GetNumberOfPoints():
            raise Exception(
                "Supplied number of colors incompatible with number of points"
            )

        if nDim != 3:
            raise Exception("Expected Nx3 array of colors")

        # Optimized version of creating the scalars array
        # - We need to be sure to keep the NumPy array,
        #   because the numpy_support function creates a pointer into it
        # - We need to take a copy to be sure the array is contigous
        self.colors_npy = cdata.copy().astype(np.uint8)
        self.colors = numpy_support.numpy_to_vtk(self.colors_npy)

        self.pd.GetPointData().SetScalars(self.colors)
        self.pd.Modified() 

def AddNormals(self, ndata):
        "Add surface normals (Nx3 NumPy array) to the current point cloud visualization object"

        nNormals = ndata.shape[0]
        nDim = ndata.shape[1]

        if nDim != 3:
            raise Exception("Expected Nx3 array of surface normals")

        if self.pd.GetNumberOfPoints() == 0:
            raise Exception("No points to add normals for")

        if nNormals != self.pd.GetNumberOfPoints():
            raise Exception(
                "Supplied number of normals incompatible with number of points"
            )

        # Optimized version of creating the scalars array
        # - We need to be sure to keep the NumPy array,
        #   because the numpy_support function creates a pointer into it
        # - We need to take a copy to be sure the array is contigous
        self.normals_npy = ndata.copy()
        self.normals = numpy_support.numpy_to_vtk(self.normals_npy)

        self.pd.GetPointData().SetNormals(self.normals)
        self.pd.Modified() 

def add_data(output, data):
    """Adds data arrays to an output VTK data object"""
    for d in data:
        arr = d.array.array
        c = nps.numpy_to_vtk(num_array=arr)
        c.SetName(d.name)
        loc = d.location
        if loc == 'vertices':
            output.GetPointData().AddArray(c)
        else:
            output.GetCellData().AddArray(c)
    return output 

def point_set_to_vtk(pse):
    """Convert the point set to a :class:`pyvista.PolyData` data object.

    Args:
        pse (:class:`omf.pointset.PointSetElement`): The point set to convert

    Return:
        :class:`pyvista.PolyData`
    """

    points = pse.geometry.vertices
    npoints = pse.geometry.num_nodes

    # Make VTK cells array
    cells = np.hstack((np.ones((npoints, 1)),
                       np.arange(npoints).reshape(-1, 1)))
    cells = np.ascontiguousarray(cells, dtype=np.int64)
    cells = np.reshape(cells, (2*npoints))
    vtkcells = vtk.vtkCellArray()
    vtkcells.SetCells(npoints, nps.numpy_to_vtk(cells, deep=True, array_type=vtk.VTK_ID_TYPE))

    # Convert points to vtk object
    pts = vtk.vtkPoints()
    pts.SetNumberOfPoints(pse.geometry.num_nodes)
    pts.SetData(nps.numpy_to_vtk(points))

    # Create polydata
    output = vtk.vtkPolyData()
    output.SetPoints(pts)
    output.SetVerts(vtkcells)

    # Now add point data:
    add_data(output, pse.data)

    add_textures(output, pse.textures, pse.name)

    return pyvista.wrap(output) 

def surface_geom_to_vtk(surfgeom):
    """Convert the triangulated surface to a :class:`pyvista.UnstructuredGrid`
    object

    Args:
        surfgeom (:class:`omf.surface.SurfaceGeometry`): the surface geomotry to
            convert
    """

    output = vtk.vtkUnstructuredGrid()
    pts = vtk.vtkPoints()
    cells = vtk.vtkCellArray()

    # Generate the points
    pts.SetNumberOfPoints(surfgeom.num_nodes)
    pts.SetData(nps.numpy_to_vtk(surfgeom.vertices))

    # Generate the triangle cells
    cellConn = surfgeom.triangles.array
    cellsMat = np.concatenate(
        (np.ones((cellConn.shape[0], 1), dtype=np.int64)*cellConn.shape[1], cellConn),
        axis=1).ravel()
    cells = vtk.vtkCellArray()
    cells.SetNumberOfCells(cellConn.shape[0])
    cells.SetCells(cellConn.shape[0],
            nps.numpy_to_vtk(cellsMat, deep=True, array_type=vtk.VTK_ID_TYPE))

    # Add to output
    output.SetPoints(pts)
    output.SetCells(vtk.VTK_TRIANGLE, cells)
    return pyvista.wrap(output) 

def set_scalar_field_cmap(self, cmap: Union[str, dict] = 'viridis',
                              regular_grid_actor = None) -> None:
        """
        Args:
            cmap:
            regular_grid_actor (Union[None, vtkRenderingOpenGL2Python.vtkOpenGLActor):
        """
        if regular_grid_actor is None:
            regular_grid_actor = self.regular_grid_actor

        if type(cmap) is dict:
            self._cmaps = {**self._cmaps, **cmap}
            cmap = cmap.keys()
        elif type(cmap) is str:
            if cmap == 'lith':
                hex_colors = list(self._get_color_lot(is_faults=False))
                n_colors = len(hex_colors)
                cmap_ = mcolors.ListedColormap(hex_colors)
                col = cmap_(np.linspace(0, 1, n_colors)) * 255
                self._cmaps[cmap] = numpy_to_vtk(col, array_type=3)
            if cmap not in self._cmaps.keys():
                col = matplotlib.cm.get_cmap(cmap)(np.linspace(0, 1, 250)) * 255
                nv = numpy_to_vtk(col, array_type=3)
                self._cmaps[cmap] = nv
        else:
            raise AttributeError('cmap must be either a name of a matplotlib string or a dictionary containing the '
                                 'rgb values')
        # Set the scalar field color map
        regular_grid_actor.GetMapper().GetLookupTable().SetTable(self._cmaps[cmap]) 

def __init__(self, nparray):
        super(VTKActorWrapper, self).__init__()

        self.nparray = nparray

        nCoords = nparray.shape[0]
        nElem = nparray.shape[1]

        self.verts = vtk.vtkPoints()
        self.cells = vtk.vtkCellArray()
        self.scalars = None

        self.pd = vtk.vtkPolyData()
        self.verts.SetData(vtk_np.numpy_to_vtk(nparray))
        self.cells_npy = np.vstack([np.ones(nCoords,dtype=np.int64),
                               np.arange(nCoords,dtype=np.int64)]).T.flatten()
        self.cells.SetCells(nCoords,vtk_np.numpy_to_vtkIdTypeArray(self.cells_npy))
        self.pd.SetPoints(self.verts)
        self.pd.SetVerts(self.cells)

        self.mapper = vtk.vtkPolyDataMapper()
        self.mapper.SetInputDataObject(self.pd)

        self.actor = vtk.vtkActor()
        self.actor.SetMapper(self.mapper)
        self.actor.GetProperty().SetRepresentationToPoints()
        self.actor.GetProperty().SetColor(0.0,1.0,0.0) 

def getVtkFromNumpy(numpyArray):

    def MakeCallback(numpyArray):
        def Closure(caller, event):
            closureArray = numpyArray
        return Closure

    vtkArray = numpy_support.numpy_to_vtk(numpyArray)
    vtkArray.AddObserver('DeleteEvent', MakeCallback(numpyArray))
    return vtkArray 

def save_vtk_labels(filename, tracts, scalars, lines_indices=None):

	lengths = [len(p) for p in tracts]
	line_starts = ns.numpy.r_[0, ns.numpy.cumsum(lengths)]
	if lines_indices is None:
		lines_indices = [ns.numpy.arange(length) + line_start for length, line_start in izip(lengths, line_starts)]

	ids = ns.numpy.hstack([ns.numpy.r_[c[0], c[1]] for c in izip(lengths, lines_indices)])
	vtk_ids = ns.numpy_to_vtkIdTypeArray(ids, deep=True)

	cell_array = vtk.vtkCellArray()
	cell_array.SetCells(len(tracts), vtk_ids)
	points = ns.numpy.vstack(tracts).astype(ns.get_vtk_to_numpy_typemap()[vtk.VTK_DOUBLE])
	points_array = ns.numpy_to_vtk(points, deep=True)

	poly_data = vtk.vtkPolyData()
	vtk_points = vtk.vtkPoints()
	vtk_points.SetData(points_array)
	poly_data.SetPoints(vtk_points)
	poly_data.SetLines(cell_array)
	poly_data.GetPointData().SetScalars(ns.numpy_to_vtk(scalars))
	poly_data.BuildCells()
#    poly_data.SetScalars(scalars)

	if filename.endswith('.xml') or filename.endswith('.vtp'):
		writer = vtk.vtkXMLPolyDataWriter()
		writer.SetDataModeToBinary()
	else:
		writer = vtk.vtkPolyDataWriter()
		writer.SetFileTypeToBinary()


	writer.SetFileName(filename)
	if hasattr(vtk, 'VTK_MAJOR_VERSION') and vtk.VTK_MAJOR_VERSION > 5:
		writer.SetInputData(poly_data)
	else:
		writer.SetInput(poly_data)
	writer.Write() 

def render(self, **kwargs):
        """ Plots the volume and the control points. """
        # Calling parent function
        super(VisVolume, self).render(**kwargs)

        # Initialize a list to store VTK actors
        vtk_actors = []

        # Start plotting
        for plot in self._plots:
            # Plot control points
            if plot['type'] == 'ctrlpts' and self.vconf.display_ctrlpts:
                # Points as spheres
                pts = np.array(plot['ptsarr'], dtype=np.float)
                vtkpts = numpy_to_vtk(pts, deep=False, array_type=VTK_FLOAT)
                vtkpts.SetName(plot['name'])
                temp_actor = vtkh.create_actor_pts(pts=vtkpts, color=vtkh.create_color(plot['color']),
                                                   name=plot['name'], index=plot['idx'])
                vtk_actors.append(temp_actor)

            # Plot evaluated points
            if plot['type'] == 'evalpts' and self.vconf.display_evalpts:
                pts = np.array(plot['ptsarr'], dtype=np.float)
                vtkpts = numpy_to_vtk(pts, deep=False, array_type=VTK_FLOAT)
                vtkpts.SetName(plot['name'])
                temp_actor = vtkh.create_actor_pts(pts=vtkpts, color=vtkh.create_color(plot['color']),
                                                   name=plot['name'], index=plot['idx'])
                vtk_actors.append(temp_actor)

        # Render actors
        return vtkh.create_render_window(vtk_actors, dict(KeyPressEvent=(self.vconf.keypress_callback, 1.0)),
                                         figure_size=self.vconf.figure_size) 

def render(self, **kwargs):
        """ Plots the volume and the control points. """
        # Calling parent function
        super(VisVoxel, self).render(**kwargs)

        # Initialize a list to store VTK actors
        vtk_actors = []

        # Start plotting
        for plot in self._plots:
            # Plot control points
            if plot['type'] == 'ctrlpts' and self.vconf.display_ctrlpts:
                # Points as spheres
                pts = np.array(plot['ptsarr'], dtype=np.float)
                vtkpts = numpy_to_vtk(pts, deep=False, array_type=VTK_FLOAT)
                vtkpts.SetName(plot['name'])
                temp_actor = vtkh.create_actor_pts(pts=vtkpts, color=vtkh.create_color(plot['color']),
                                                   name=plot['name'], index=plot['idx'])
                vtk_actors.append(temp_actor)

            # Plot evaluated points
            if plot['type'] == 'evalpts' and self.vconf.display_evalpts:
                faces = np.array(plot['ptsarr'][1], dtype=np.float)
                filled = np.array(plot['ptsarr'][2], dtype=np.int)
                grid_filled = faces[filled == 1]
                temp_actor = vtkh.create_actor_hexahedron(grid=grid_filled, color=vtkh.create_color(plot['color']),
                                                          name=plot['name'], index=plot['idx'])
                vtk_actors.append(temp_actor)

        # Render actors
        return vtkh.create_render_window(vtk_actors, dict(KeyPressEvent=(self.vconf.keypress_callback, 1.0)),
                                         figure_size=self.vconf.figure_size) 

def numpy_to_image(numpy_array):
    """
    @brief Convert a numpy 2D or 3D array to a vtkImageData object
    @param numpy_array 2D or 3D numpy array containing image data
    @return vtkImageData with the numpy_array content
    """

    shape = numpy_array.shape
    if len(shape) < 2:
        raise Exception('numpy array must have dimensionality of at least 2')

    h, w = shape[0], shape[1]
    c = 1
    if len(shape) == 3:
        c = shape[2]

    # Reshape 2D image to 1D array suitable for conversion to a
    # vtkArray with numpy_support.numpy_to_vtk()
    linear_array = np.reshape(numpy_array, (w*h, c))
    vtk_array = numpy_to_vtk(linear_array)

    image = vtkImageData()
    image.SetDimensions(w, h, 1)
    image.AllocateScalars()
    image.GetPointData().GetScalars().DeepCopy(vtk_array)

    return image 

def vtk_points(points, deep=True):
    """Convert numpy points to a vtkPoints object."""
    if not points.flags['C_CONTIGUOUS']:
        points = np.ascontiguousarray(points)
    vtkpts = vtk.vtkPoints()
    vtkpts.SetData(nps.numpy_to_vtk(points, deep=deep))
    return vtkpts 

def _from_arrays(self, x, y, z):
        """Create VTK rectilinear grid directly from numpy arrays.

        Each array gives the uniques coordinates of the mesh along each axial
        direction. To help ensure you are using this correctly, we take the unique
        values of each argument.

        Parameters
        ----------
        x : np.ndarray
            Coordinates of the nodes in x direction.

        y : np.ndarray
            Coordinates of the nodes in y direction.

        z : np.ndarray
            Coordinates of the nodes in z direction.

        """
        # Set the coordinates along each axial direction
        # Must at least be an x array
        x = np.unique(x.ravel())
        self.SetXCoordinates(numpy_to_vtk(x))
        if y is not None:
            y = np.unique(y.ravel())
            self.SetYCoordinates(numpy_to_vtk(y))
        if z is not None:
            z = np.unique(z.ravel())
            self.SetZCoordinates(numpy_to_vtk(z))
        # Ensure dimensions are properly set
        self._update_dimensions() 

def x(self, coords):
        """Set the coordinates along the X-direction."""
        self.SetXCoordinates(numpy_to_vtk(coords))
        self._update_dimensions()
        self.Modified() 

def z(self, coords):
        """Set the coordinates along the Z-direction."""
        self.SetZCoordinates(numpy_to_vtk(coords))
        self._update_dimensions()
        self.Modified() 

def save_vtk(filename, tracts, lines_indices=None, scalars = None):

	lengths = [len(p) for p in tracts]
	line_starts = ns.numpy.r_[0, ns.numpy.cumsum(lengths)]
	if lines_indices is None:
		lines_indices = [ns.numpy.arange(length) + line_start for length, line_start in izip(lengths, line_starts)]

	ids = ns.numpy.hstack([ns.numpy.r_[c[0], c[1]] for c in izip(lengths, lines_indices)])
	vtk_ids = ns.numpy_to_vtkIdTypeArray(ids, deep=True)

	cell_array = vtk.vtkCellArray()
	cell_array.SetCells(len(tracts), vtk_ids)
	points = ns.numpy.vstack(tracts).astype(ns.get_vtk_to_numpy_typemap()[vtk.VTK_DOUBLE])
	points_array = ns.numpy_to_vtk(points, deep=True)

	poly_data = vtk.vtkPolyData()
	vtk_points = vtk.vtkPoints()
	vtk_points.SetData(points_array)
	poly_data.SetPoints(vtk_points)
	poly_data.SetLines(cell_array)
	poly_data.BuildCells()

	if filename.endswith('.xml') or filename.endswith('.vtp'):
		writer = vtk.vtkXMLPolyDataWriter()
		writer.SetDataModeToBinary()
	else:
		writer = vtk.vtkPolyDataWriter()
		writer.SetFileTypeToBinary()


	writer.SetFileName(filename)
	if hasattr(vtk, 'VTK_MAJOR_VERSION') and vtk.VTK_MAJOR_VERSION > 5:
		writer.SetInputData(poly_data)
	else:
		writer.SetInput(poly_data)
	writer.Write() 

def add_scalar(var, timestamp):

    vtkFile = vtk.vtkXMLUnstructuredGridReader()
    vtkFile.SetFileName('cosipy.vtu')
    vtkFile.Update()
    
    # Find cellId by coordinates
    pointLocator =  vtk.vtkPointLocator()
    pointLocator.SetDataSet(vtkFile.GetOutput())
    pointLocator.BuildLocator()
    
    ds = xr.open_dataset('../data/output/Peru_20160601-20180530_comp4.nc')
    ds = ds.sel(time=timestamp)
    
    ds_sub = ds[var].stack(x=['south_north','west_east']) 
    ds_sub = ds_sub.dropna(dim='x')
    lats = ds_sub.x.lat.values
    lons = ds_sub.x.lon.values
    data = ds_sub.values
    print(lats)

    numPoints = vtkFile.GetOutput().GetNumberOfPoints()
    scalar = np.empty(numPoints)
    scalar[:] = np.nan

    interpField = numpy_support.numpy_to_vtk(scalar)
    interpField.SetName(var)
    vtkFile.GetOutput().GetPointData().AddArray(interpField)
    vtkFile.Update()

    print('Write points \n')
    for i in np.arange(len(data)):
        # Get height
        alt = ds.HGT.sel(lat=lats[i],lon=lons[i]).values/6370000.0
        
        pointId = vtk.mutable(0) 
        Id = vtk.vtkIdList()
        pointId = pointLocator.FindClosestPoint([lons[i],lats[i],alt])
        vtkFile.GetOutput().GetPointData().GetArray(var).InsertTuple1(pointId,data[i])

    writer = vtk.vtkXMLUnstructuredGridWriter()
    writer.SetFileName('cosipy.vtu')
    writer.SetInputData(vtkFile.GetOutput())
    writer.Write()

    #plotSurface(vtkFile) 

def plot_topography(self, topography = None, scalars='geo_map', **kwargs):
        """

        Args:
            topography: gp Topography object, np.array or None
            scalars:
            **kwargs:

        Returns:

        """
        if topography is None:
            topography = self.model._grid.topography.values
        rgb = False

        # Create vtk object
        cloud = pv.PolyData(topography)

        # Set scalar values
        if scalars == 'geo_map':
            arr_ = np.empty((0, 3), dtype='int')

            # Convert hex colors to rgb
            for val in list(self._color_lot):
                rgb = (255 * np.array(mcolors.hex2color(val)))
                arr_ = np.vstack((arr_, rgb))

            sel = np.round(self.model.solutions.geological_map[0]).astype(int)[0]
          #  print(arr_)
          #  print(sel)

            scalars_val = numpy_to_vtk(arr_[sel-1], array_type=3)
            cm = None
            rgb = True

        elif scalars == 'topography':
            scalars_val = topography[:, 2]
            cm = 'terrain'

        elif type(scalars) is np.ndarray:
            scalars_val = scalars
            scalars = 'custom'
            cm = 'terrain'
        else:
            raise AttributeError()

        topo_actor = self.p.add_mesh(cloud.delaunay_2d(), scalars=scalars_val, cmap=cm, rgb=rgb, **kwargs)
        self.vista_topo_actors[scalars] = topo_actor
        return topo_actor 

def plot_topography(
            self,
            topography = None,
            scalars="geomap",
            **kwargs
    ):
        if not topography:
            try:
                topography = self.model._grid.topography.values
            except AttributeError:
                print("Unable to plot topography: Given geomodel instance "
                      "does not contain topography grid.")
                return

        polydata = pv.PolyData(topography)

        rgb = False
        if scalars == "geomap":
            arr_ = np.empty((0, 3), dtype=int)
            # convert hex colors to rgb
            for val in list(self._get_color_lot(faults=False)):
                rgb = (255 * np.array(mcolors.hex2color(val)))
                arr_ = np.vstack((arr_, rgb))

            sel = np.round(self.model.solutions.geological_map[0]).astype(int)[0]

            scalars_val = numpy_to_vtk(arr_[sel - 1], array_type=3)
            cm = None
            rgb = True
        elif scalars == "topography":
            scalars_val = topography[:, 2]
            cm = 'terrain'
        elif type(scalars) is np.ndarray:
            scalars_val = scalars
            scalars = 'custom'
            cm = 'terrain'
        else:
            raise AttributeError("Parameter scalars needs to be either \
                'geomap', 'topography' or a np.ndarray with scalar values")

        topography_actor = self.p.add_mesh(
            polydata.delaunay_2d(),
            scalars=scalars_val,
            cmap=cm,
            rgb=rgb,
            **kwargs
        )
        self._surface_actors["topography"] = topography_actor
        return topography_actor 

def set_topography(self):
        # Create points on an XY grid with random Z coordinate
        vertices = copy.copy(self.geo_model._grid.topography.values)

        points = vtk.vtkPoints()
        # for v in vertices:
        #     v[-1] = v[-1]
        #     points.InsertNextPoint(v)
        if self.ve !=1:
            vertices[:, 2]= vertices[:, 2]*self.ve
        points.SetData(numpy_to_vtk(vertices))

        # Add the grid points to a polydata object
        polydata = vtk.vtkPolyData()
        polydata.SetPoints(points)
        #
        # glyphFilter = vtk.vtkVertexGlyphFilter()
        # glyphFilter.SetInputData(polydata)
        # glyphFilter.Update()
        #
        # # Create a mapper and actor
        # pointsMapper = vtk.vtkPolyDataMapper()
        # pointsMapper.SetInputConnection(glyphFilter.GetOutputPort())

        #
        # pointsActor = vtk.vtkActor()
        # pointsActor.SetMapper(pointsMapper)
        # pointsActor.GetProperty().SetPointSize(3)
        # pointsActor.GetProperty().SetColor(colors.GetColor3d("Red"))

        # Triangulate the grid points
        delaunay = vtk.vtkDelaunay2D()
        delaunay.SetInputData(polydata)
        delaunay.Update()

        # Create a mapper and actor
        triangulatedMapper = vtk.vtkPolyDataMapper()
        triangulatedMapper.SetInputConnection(delaunay.GetOutputPort())

        triangulatedActor = vtk.vtkActor()
        triangulatedActor.SetMapper(triangulatedMapper)

        self.topography_surface = triangulatedActor
        self._topography_polydata = polydata
        self._topography_delauny = delaunay
        self.ren_list[0].AddActor(triangulatedActor)
        self.ren_list[1].AddActor(triangulatedActor)
        self.ren_list[2].AddActor(triangulatedActor)
        self.ren_list[3].AddActor(triangulatedActor)
        try:
            if self.geo_model.solutions.geological_map is not None:
                self.set_geological_map()
        except AttributeError as ae:
            warnings.warn(str(ae)) 

def vector_poly_data(orig, vec):
    """Create a vtkPolyData object composed of vectors."""
    # shape, dimension checking
    if not isinstance(orig, np.ndarray):
        orig = np.asarray(orig)

    if not isinstance(vec, np.ndarray):
        vec = np.asarray(vec)

    if orig.ndim != 2:
        orig = orig.reshape((-1, 3))
    elif orig.shape[1] != 3:
        raise ValueError('orig array must be 3D')

    if vec.ndim != 2:
        vec = vec.reshape((-1, 3))
    elif vec.shape[1] != 3:
        raise ValueError('vec array must be 3D')

    # Create vtk points and cells objects
    vpts = vtk.vtkPoints()
    vpts.SetData(nps.numpy_to_vtk(np.ascontiguousarray(orig), deep=True))

    npts = orig.shape[0]
    cells = np.empty((npts, 2), dtype=pyvista.ID_TYPE)
    cells[:, 0] = 1
    cells[:, 1] = np.arange(npts, dtype=pyvista.ID_TYPE)
    vcells = pyvista.utilities.cells.CellArray(cells, npts)

    # Create vtkPolyData object
    pdata = vtk.vtkPolyData()
    pdata.SetPoints(vpts)
    pdata.SetVerts(vcells)

    # Add vectors to polydata
    name = 'vectors'
    vtkfloat = nps.numpy_to_vtk(np.ascontiguousarray(vec), deep=True)
    vtkfloat.SetName(name)
    pdata.GetPointData().AddArray(vtkfloat)
    pdata.GetPointData().SetActiveVectors(name)

    # Add magnitude of vectors to polydata
    name = 'mag'
    scalars = (vec * vec).sum(1)**0.5
    vtkfloat = nps.numpy_to_vtk(np.ascontiguousarray(scalars), deep=True)
    vtkfloat.SetName(name)
    pdata.GetPointData().AddArray(vtkfloat)
    pdata.GetPointData().SetActiveScalars(name)

    return pyvista.PolyData(pdata) 

def surface_grid_geom_to_vtk(surfgridgeom):
    """Convert the 2D grid to a :class:`pyvista.StructuredGrid` object.

    Args:
        surfgridgeom (:class:`omf.surface.SurfaceGridGeometry`): the surface
            grid geometry to convert

    """
    surfgridgeom._validate_mesh()

    output = vtk.vtkStructuredGrid()

    axis_u = np.array(surfgridgeom.axis_u)
    axis_v = np.array(surfgridgeom.axis_v)
    axis_w = np.cross(axis_u, axis_v)
    if not check_orthogonal(axis_u, axis_v, axis_w):
        raise ValueError('axis_u, axis_v, and axis_w must be orthogonal')
    rotation_mtx = np.array([axis_u, axis_v, axis_w])
    ox, oy, oz = surfgridgeom.origin

    # Make coordinates along each axis
    x = ox + np.cumsum(surfgridgeom.tensor_u)
    x = np.insert(x, 0, ox)
    y = oy + np.cumsum(surfgridgeom.tensor_v)
    y = np.insert(y, 0, oy)

    z = np.array([oz])

    output.SetDimensions(len(x), len(y), len(z))

    # Build out all nodes in the mesh
    xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')
    xx, yy, zz, = xx.ravel('F'), yy.ravel('F'), zz.ravel('F')
    zz += surfgridgeom.offset_w
    points = np.c_[xx, yy, zz]

    # Rotate the points based on the axis orientations
    points = points.dot(rotation_mtx)

    # Convert points to vtk object
    pts = vtk.vtkPoints()
    pts.SetNumberOfPoints(len(points))
    pts.SetData(nps.numpy_to_vtk(points))
    # Now build the output
    output.SetPoints(pts)

    return pyvista.wrap(output) 

def volume_grid_geom_to_vtk(volgridgeom):
    """Convert the 3D gridded volume to a :class:`pyvista.StructuredGrid`
    (or a :class:`pyvista.RectilinearGrid` when apprropriate) object contatining
    the 2D surface.

    Args:
        volgridgeom (:class:`omf.volume.VolumeGridGeometry`): the grid geometry
            to convert
    """
    volgridgeom._validate_mesh()

    ox, oy, oz = volgridgeom.origin

    # Make coordinates along each axis
    x = ox + np.cumsum(volgridgeom.tensor_u)
    x = np.insert(x, 0, ox)
    y = oy + np.cumsum(volgridgeom.tensor_v)
    y = np.insert(y, 0, oy)
    z = oz + np.cumsum(volgridgeom.tensor_w)
    z = np.insert(z, 0, oz)

    # If axis orientations are standard then use a vtkRectilinearGrid
    if check_orientation(volgridgeom.axis_u, volgridgeom.axis_v, volgridgeom.axis_w):
        output = vtk.vtkRectilinearGrid()
        output.SetDimensions(len(x), len(y), len(z)) # note this subtracts 1
        output.SetXCoordinates(nps.numpy_to_vtk(num_array=x))
        output.SetYCoordinates(nps.numpy_to_vtk(num_array=y))
        output.SetZCoordinates(nps.numpy_to_vtk(num_array=z))
        return pyvista.wrap(output)

    # Otherwise use a vtkStructuredGrid
    output = vtk.vtkStructuredGrid()
    output.SetDimensions(len(x), len(y), len(z)) # note this subtracts 1

    # Build out all nodes in the mesh
    xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')
    points = np.c_[xx.ravel('F'), yy.ravel('F'), zz.ravel('F')]

    # Rotate the points based on the axis orientations
    rotation_mtx = np.array([volgridgeom.axis_u, volgridgeom.axis_v, volgridgeom.axis_w])
    points = points.dot(rotation_mtx)

    # Convert points to vtk object
    pts = vtk.vtkPoints()
    pts.SetNumberOfPoints(len(points))
    pts.SetData(nps.numpy_to_vtk(points))
    # Now build the output
    output.SetPoints(pts)

    return pyvista.wrap(output)