You can just use the 3 argument form of np.where for this:

>>> import numpy as np
>>> x = np.array([1,2,31,32,4,0,3,0,0,0])
>>> z = np.array([99] * len(x))
>>> y = np.where(x != 0, x, z)
>>> y
array([ 1,  2, 31, 32,  4, 99,  3, 99, 99, 99])

Alpha blending is usually done using the Porter & Duff equations:

where src and dst would correspond to your foreground and background images, and the A and RGB pixel values are assumed to be floating point, in the range [0, 1].

For your specific example:

src_rgb = fg[..., :3].astype(np.float32) / 255.0
src_a = fg[..., 3].astype(np.float32) / 255.0
dst_rgb = bg[..., :3].astype(np.float32) / 255.0
dst_a = bg[..., 3].astype(np.float32) / 255.0

out_a = src_a + dst_a*(1.0-src_a)
out_rgb = (src_rgb*src_a[..., None]
           + dst_rgb*dst_a[..., None]*(1.0-src_a[..., None])) / out_a[..., None]

out = np.zeros_like(bg)
out[..., :3] = out_rgb * 255
out[..., 3] = out_a * 255

Output:

You can use slicing to replace first len(smaller_array) elements in larger_array with smaller_array as follows:

arr[:len(arr2)] = arr2

Why don't you use imshow instead?

You can plot a 2D image by doing:

plt.imshow(Image1, cmap='gray') # I would add interpolation='none'

Afterwards, you can easily overlay the segmentation by doing:

plt.imshow(Image2_mask, cmap='jet', alpha=0.5) # interpolation='none'

Changing the alpha will change the opacity of the overlay.

Additionaly, why do you create 2 masks? Only one should be enough, you can do:

Image2_mask = ma.masked_array(Image2 > 0, Image2)

Practical example:

import numpy as np
mask = np.zeros((10,10))
mask[3:-3, 3:-3] = 1 # white square in black background
im = mask + np.random.randn(10,10) * 0.01 # random image
masked = np.ma.masked_where(mask == 0, mask)

import matplotlib.pyplot as plt
plt.figure()
plt.subplot(1,2,1)
plt.imshow(im, 'gray', interpolation='none')
plt.subplot(1,2,2)
plt.imshow(im, 'gray', interpolation='none')
plt.imshow(masked, 'jet', interpolation='none', alpha=0.7)
plt.show()

Because you are limited to writing contiguous blocks of data to GDAL rasters, the best way around this is to create the output array first, then write the array to the output raster. Using this method, you can use numpy's boolean indexing to dictate what is written. Your last for loop would work as follows using this logic:

import numpy

no_data = -2147483647
output_array = numpy.full((newheight, newwidth), no_data, 'int32')

for f in files:
    raster = gdal.Open(f, GA_ReadOnly)
    f_xmin, f_pwidth, f_xskew, f_ymax, f_yskew, f_pheight = raster.GetGeoTransform()
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    xoffset = int((f_xmin - xmin) / f_pwidth)
    yoffset = int((f_ymax - ymax) / f_pheight)  # Assumes grids align
    band = raster.GetRasterBand(1)
    data = band.ReadAsArray()
    # Create a boolean array that allows data where they exist
    mask = (output_array == no_data) & (data != band.GetNoDataValue())  # Note: this assumes you prefer data in order based on the files list
    # Allocate data to output array
    output_array[yoffset:yoffset + rows, xoffset:xoffset + cols][mask] = data[mask]

# Write data to the output raster
outband.WriteArray(output_array)
outband.FlushCache()

However, you will need to fully load the raster in to memory prior to flushing it to the disk. If memory-management is an issue, you can either use numpy.memmap, or read the overlapping block from the output raster for each iteration, which would look like this:

for f in files:
    raster = gdal.Open(f, GA_ReadOnly)
    f_xmin, f_pwidth, f_xskew, f_ymax, f_yskew, f_pheight = raster.GetGeoTransform() #geotransform for the old files being iterated over
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    xoffset = int((f_xmin - xmin) / f_pwidth)
    yoffset = int((f_ymax - ymax) / f_pheight)
    band = raster.GetRasterBand(1)
    data = band.ReadAsArray(0, 0, cols, rows)
    data_update = outband.ReadAsArray(xoffset, yoffset, cols, rows)
    mask = (data_update == no_data) & (data != band.GetNoDataValue())
    data_update[mask] = data[mask]
    outband.WriteArray(data_update, xoffset, yoffset)
    outband.FlushCache()

One last general note of caution: your method assumes that raster grids and data types align.

Your issue is trying to work with lists or object arrays.

Numpy is simply not designed for that. Use a flat array instead.

# convert subixs to a flat array
flat_subixs = np.concatenate(subixs)
# compute the lengths once
lengths = np.array([len(l) for l in subixs])
arr = np.repeat(overlay, lengths)
print(arr)

msk = (overlay == 1) | (overlay == 3)
arr[flat_subixs[np.repeat(msk, lengths)]] = np.repeat([44, 48, 47], lengths[msk])
print(arr)

Output:

[0 0 0 0 1 1 1 4 3]
[ 0  0  0  0 44 44 48  4 47]

This way all the slow computations (concatenation and lengths) are performed only once, and all other operations are vectorized.

If you want to mutate overlay this is also possible, just create arr once in the end:

msk = (overlay == 1) | (overlay == 3)
overlay[msk] = [44, 48, 47]
# eventually add other transforms here

# finally create "arr" once
arr = np.repeat(overlay, lengths)
print(arr)
# [ 0  0  0  0 44 44 48  4 47]

You need to make your arrays numpy arrays so that np.where works correctly:

M = np.array([
    [bg,bg,bg,bg,bg,bg],
    [bg,bb,bb,bb,bg,bg],
    [bg,bb,bb,bb,bg,bg],
    [bg,bg,bg,bg,bg,bg],
    [bg,bg,bg,bg,bg,bg],
    [bg,bg,bg,bg,bg,bg]
    ])

L = np.array([
    [tr,tr,tr,tr,tr,tr],
    [tr,tr,tr,tr,tr,tr],
    [tr,tr,tr,tr,tr,tr],
    [tr,tr,tr,tl,tr,tr],
    [tr,tr,tr,tl,tr,tr],
    [tr,tr,tr,tr,tr,tr]
    ])