Q:

2d convolution gives not the desired output

I want to use the 2D convolution in the same way I did here in 1D. Unfortunately the output in the former case does not have the desired shape. Let n = 5, then
h_0 = (1 / 4) * np.array([1, 2, 1])
x = np.random.rand(n)
np.convolve(h_0, x, 'same')

>>> array([0.65498075, 0.72729356, 0.51417706, 0.34597679, 0.1793755])

but
h_00 = np.kron(h_0, h_0)
h_00 = np.reshape(h_00, (3, 3))
x = np.random.rand(n, n)
scipy.signal.convolve2d(h_00, x, 'same', boundary='symm')

>>> array([[1.90147294, 1.6541233 , 1.82704077],
           [1.55228912, 1.3641027 , 1.55536069],
           [1.61190909, 1.45159935, 1.58266083]])

I would have expected a (5, 5) output array.

A:

The docs for scipy.signal.convolve2d regarding the mode parameter clearly state

mode
...
same
    The output is the same size as in1, centered with respect to the ‘full’ output.

So, given that you pass the kernel first, your output will be the same size as the kernel, not the array you are filtering. To fix, swap the first two inputs:
scipy.signal.convolve2d(x, h_00, 'same', boundary='symm')

Confusion likely arises from the behavior of numpy.convolve, which does the following:

mode : {‘full’, ‘valid’, ‘same’}, optional
...
‘same’:
    Mode ‘same’ returns output of length max(M, N). Boundary effects are still visible.

Numpy interprets the larger array as the kernel regardless of argument order. This is possible because with a single dimension, there is always an unambiguous winner.

