Usually to render a large quantity of gird aligned cubes is to use a method known as "Marching cubes".
Note: the main alternative is to use ray casting and voxels (point cloud data), which is less common in video games, however it is the standard in medical imaging. See also: A Fast Voxel Traversal Algorithm for Ray Tracing.
You will see "Marching cubes" come up often when people talk about improving the performance of "voxels". The idea is that you will be iterating over the array of cubes and merging them into a single mesh in such way that hidden faces are removed, and also merges co-planar adjacent faces.
However, what I have described so far is simply Meshing.
Marching cubes (no quotes this time) goes a step further. In a similar way to how you can pick a tile from a tilemap based on its neighbors (Autotiling), Marching cubes will pick polygons depending on the neighbor cubes.
See also Combinations for tiling two textures together.
The original paper on Marching Cubes describes 15 configurations:
For each cube, you will pick a mesh from these configurations using a lookup table based on what vertices are hidden or visible. Since each cube has eight vertices, we have 256 entries on the table, which are reduced to the ones above by considering symmetries.
See Marching Cubes Animation | Algorithms Visualized.
See also Nvidia GPU Gems 3 - Chapter 1.
You do not have to use the same set of cases as described here. Hopefully you can see how you could edit this to match what you want.
Finally, if we do Marching cubes minus the Meshing, you would not be merging the cubes into a single mesh... instead you would simply identify what vertices are hidden and what vertices are visible and pick/modify the meshes depending on that.
Wait, that is what you are doing. Huh.