Visualize Covariance when only probability mass and marginal functions are given

Question

I am trying to intuitively understand Covariance like here. So if a random sample set given, I could draw rectangles with them, one of the cornes being fixated on mean $(\overline{x},\overline{y})$. Something like below.

I saw an example for Covariance where just the joint pmf was given and Covariance was calculated as below.

This was weird, because if I try to visualize the points location as per formula used above, it would be perfect rectangle like below, suggesting the X and Y are independent. It does not make sense to directly draw them from joint pmf, to observe any trend, then how come they could be used in teh Covariance formula? How do I visualize this particular situation when only pmf is given? I understand inequal probability would have had given a scatter if samples are drawn from this distribution, but then we should have used sample data in Covariance formula, not these right?

The above example was taken from Probability and Statistics for engineering and the sciences by Devore link here, page 208.

Update 1:
I have constructed 3d bar as below, to illustrate that probability also adding weight (here height), to the measure, which would not be evident from 2D perspective (which is inclusive in this case). Looking at 3D, one could say, oh ok green is more (with higher bars, larger volume), so positive covariance. Will this visual inference would be correct? If this is so, only issue here is the 3d bars are not rightly rendered issue created

Answer 1

In an Answer to this question @Glen_b shows how to generate a bivariate normal sample with a given population correlation $\rho.$ Below, I use this method to show samples of size $n = 2000$ from each of three bivariate normal distributions with respective correlations $\rho = .1, -.5, .9.$ In each case $\mu_x = \mu_y = 0$ and $\sigma_x = \sigma_y = 1,$ so that all marginal distributions are standard normal.

set.seed(1103);  n=2000
rho = .1; x1 = rnorm(n);  z = rnorm(n);  y1 = rho*x1 + sqrt(1-rho^2)*z
rho =-.5; x2 = rnorm(n);  z = rnorm(n);  y2 = rho*x2 + sqrt(1-rho^2)*z
rho = .9; x3 = rnorm(n);  z = rnorm(n);  y3 = rho*x3 + sqrt(1-rho^2)*z
r1=round(cor(x1,y1),3); r2=round(cor(x2,y2),3); r3=round(cor(x3,y3),3)
par(mfrow=c(1,3))
   plot(x1,y1, pch=".", main=paste("r = ",r1))
     abline(h=0, col="green2"); abline(v=0, col="green2")
   plot(x2,y2, pch=".", main=paste("r = ",r2))
     abline(h=0, col="green2"); abline(v=0, col="green2")
   plot(x3,y3, pch=".", main=paste("r = ",r3))
    abline(h=0, col="green2"); abline(v=0, col="green2")
par(mfrow=c(1,1))

You are correct that concentrations of points in the first and third quadrants add positive terms to the sample covariance and hence tend to make the sample correlation positive; similarly, points in quadrants two and four add negative terms to the sample covariance.

However, do not confuse zero correlation with independence: if $X$ and $Y$ are independent, then $\rho_{x,y} \approx 0.$ But correlation near $0$ is not necessarily an indication of independence. (Perhaps see the Wikipedia article on correlation and dependence for examples.)