PythonScript to load data as colored point-cloud Background

I want to render 3D points with color information for scientific visualization. I have several hundred thousands of points up to 2 million. For each point I have [xyz] coordinates and two scalars [t, tau]. The scalars should be shown as color and are currently mapped to uv-coordinates.

Wanted result

I want to achieve something like this but with (lots) more points/spheres:



Current solution


Here is my script: visualisePointCloud.py
Here is some example data: point_data.txt


My workflow


load the data 
create/copy a sphere mesh
move the copied mesh to [xyz] coordinates
define the uv coordinates for the whole mesh with respect to [t, tau]
repeat 2. to 4. for each data entry
join all created meshes


Problem

While loading the points and creating the scene I get a slowdown. The first 1000 points work fine, but the overall process gets quickly slower. I am sure it has something to do with the commands I use from the API and the amount of objects I create. 

I had two ideas to solve the problem - neither worked:


Create a sphere mesh for each data entry and join them in the end. This leads to a huge amount of created objects, which slows down blender (compare this question).
Create the meshes for a subset of the data entries and join these to a single object/mesh. Do this until all points are loaded. This leads to the problem, that the joining process gets more complicated and increasingly takes more time.




Questions


Are there better commands I can use to speed up the process?
Is there a better solution than copy a sphere mesh for each data entry and adjust its UV-coordinates for coloring. E.g. create one mesh and add vertices to it - how can I alter the UV-coordinates this way.
Is there generally a better way to render a colored point
cloud, e.g. this one is missing the colors.
Updated:

This script will add the pointcloud data as mesh vertices, then create a low-poly icosphere as an instanced object on each vert (using mesh.from_pydata instead of bmesh as advised by BatFINGER).

Currently this does not use the UV data in the latter two columns in the data file. I've also scaled the coordiantes by a factor of 25 (can be changed in the script) to make everything more visible and easier to manipulate in the viewport. This script takes ~0.075 seconds on my laptop for the 10,001 verts in this file.

import bpy
import numpy as np
from mathutils import Vector
from time import time

t = time()
C = bpy.context

scale = 25

# Open textfile to numpy array and scale
pointcloudfilepath = '/Users/username/Downloads/point_data.txt'
with open( pointcloudfilepath ) as fh:
    data = np.loadtxt( fh ) * scale

print( f'vertcount = {len(data)}')

# Create and arrange mesh data
verts = [ Vector( data[i,:3] ) for i in range(data.shape[0]) ]
m     = bpy.data.meshes.new('pc')
m.from_pydata(verts, [], [])

# Create mesh object and link to scene collection
o = bpy.data.objects.new( 'pc', m )
C.scene.collection.objects.link( o )

# Add minimal icosphere
bpy.ops.mesh.primitive_ico_sphere_add( subdivisions = 1, radius = 0.05 )
isobj = bpy.data.objects[ C.object.name ]

# Set instancing props
for ob in [ isobj, o ]:
    ob.instance_type               = 'VERTS'
    ob.show_instancer_for_viewport = False
    ob.show_instancer_for_render   = False

# Set instance parenting (parent icosphere to verts)
o.select_set(True)
C.view_layer.objects.active = o

bpy.ops.object.parent_set( type = 'VERTEX', keep_transform = True )

print( f'Total time = {time() - t} seconds' )


I've used each sphere's location as material color in a simple shader, you can apply any other logic that fits your requirements better. Here's the node setup:



And final result: