Suppose you have a 1D array, you plot the position of its values on the x axis, they are so dense that you can't tell the spatial distribution, you use a 1D histogram to show the distribution by count of boxes along the x axis. Problem solved.

Then you have two 1D arrays, a list of 2D dots in (x, y) axes. You plot their positions on the x-y plane, again they are so dense and overlap with each other. You want to view the distribution better by count of boxes in the plane, so you try a 2D diagram. Problem solved.

Here is an example

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

# prepare 2D random dots centered at (0, 0)
n = 100000
x = np.random.randn(n)
y = x + np.random.randn(n)

# plot data
fig1 = plt.figure()
plt.plot(x,y,'.r')
plt.xlabel('x')
plt.ylabel('y')

gives

# plot 2D histogram using pcolor
fig2 = plt.figure()
plt.hist2d(x, y, bins=100)
plt.xlabel('x')
plt.ylabel('y')
cbar = plt.colorbar()
cbar.ax.set_ylabel('Counts')

gives

If you have the raw data from the counts, you could use plt.hexbin to create the plots for you (IMHO this is better than a square lattice): Adapted from the example of hexbin:

import numpy as np
import matplotlib.pyplot as plt

n = 100000
x = np.random.standard_normal(n)
y = 2.0 + 3.0 * x + 4.0 * np.random.standard_normal(n)
plt.hexbin(x,y)

plt.show()

If you already have the Z-values in a matrix as you mention, just use plt.imshow or plt.matshow:

XB = np.linspace(-1,1,20)
YB = np.linspace(-1,1,20)
X,Y = np.meshgrid(XB,YB)
Z = np.exp(-(X**2+Y**2))
plt.imshow(Z,interpolation='none')