Placing cubes in a grid is very slow (using python script) I'm trying to generate a scene comprising of voxels generated from a 3D data matrix.
However the script that I'm using is super slow (~1h per image), and most of the time is consumed while building the scene (pyhton loop) instead of the rendering.
Does anyone see where the problem comes from?

import bpy
import bmesh
import numpy as np
import numpy
import scipy.io
import sys

# Data loading
inputFile = sys.argv[-2]
outputFile = sys.argv[-1]
data = scipy.io.loadmat(inputFile)['data']  
voxel_dims = np.array([256, 256, 32])
data = np.reshape(data, (voxel_dims))

# Create blender scene
context = bpy.context
scene = context.scene

def clear_scene(scene):
    #Clears the current scene
    for obj in scene.objects:
        scene.objects.unlink(obj)
        bpy.data.objects.remove(obj)
    return

def base_object():
    #Creates base cube
    bpy.ops.mesh.primitive_cube_add(radius=0.4, location=(0,0,0))
    return scene.objects.active

colors = {  
          '_1':(1.0,0.0,0.0),
          '_2':(0.0,1.0,0.0),,
          '_3':(0.0,0.0,1.0),
         }
color_map = {
             1: '_1',
             2: '_2',
             3: '_3'
            }        

def creatematerial(name, color):
    red = color[0]
    green = color[1]
    blue = color[2]
    alpha = 1.0
    mat=(bpy.data.materials.new(name=name))    
    mat.use_nodes = True
    Diffuse_BSDF = mat.node_tree.nodes['Diffuse BSDF']
    Diffuse_BSDF.inputs[0].default_value = [red, green, blue, alpha]
    mat.diffuse_color = [red, green, blue]
    return mat

clear_scene(scene) 
bpy.context.scene.render.engine = 'CYCLES'
bpy.ops.object.empty_add(type='PLAIN_AXES', location=(0,0,0))
empty = bpy.context.active_object
materials=dict()    
for k,v, in colors.items():
    materials[k] = creatematerial(k,v)
bo = base_object()
mesh = bo.data

# Scene creation (This is the loop that takes forever!)
for idx,label in enumerate(color_map):
    col = colors[color_map[label]]
    xx, yy, zz = np.where(data == label)
    bo.data.materials.clear()
    bo.data.materials.append(materials[color_map[label]])
    for i in range(len(xx)):
        (x,y,z) = (xx[i], yy[i], zz[i])
        thiscube = bo.copy()
        thiscube.data = thiscube.data.copy()
        scene.objects.link(thiscube)
        thiscube.location=(x,y,z)
        thiscube.parent=empty


scene.objects.unlink(bo)
bpy.data.objects.remove(bo)

# Place camera and render scene 
bpy.ops.object.camera_add(location=(130.409,114.09,46.750),rotation=(1.201,-0.0000464,3.111))
bpy.context.object.data.clip_end = 200
point=bpy.ops.object.lamp_add(type='POINT',location=(15,15,60))
bpy.ops.object.lamp_add(type='SUN',location=(15,15,80))
bpy.context.scene.cycles.film_transparent = True
bpy.context.scene.render.resolution_percentage = 100
bpy.context.scene.cycles.samples = 256
bpy.context.scene.render.resolution_x = 512
bpy.context.scene.render.resolution_y = 480
bpy.context.scene.cycles.device = 'GPU'
bpy.context.scene.camera = bpy.data.objects['Camera']
bpy.context.scene.render.filepath = outputFile
bpy.ops.render.render(write_still=True)


Every entry of the 3D matrix has an id which corresponds to its color. The dimension of the matrix is 256x256x32, however only a subset of the voxels is occupied.
As far as I can see the code creates for each data point a seperate Blender Mesh Object. For large datasets this will make everything super slow.
I would recommend you to create only one mesh per label-class so that with 3 label classes you would have only 3 Mesh Objects.

Since you are actually using my code, I took your question as feedback and updated it accordingly on github.
https://github.com/cagcoach/BlenderPlot