I know the pdf of a multivariate normal distribution can be expressed as
$$ {\displaystyle (2\pi )^{-{\frac {k}{2}}}\det({\boldsymbol {\Sigma }})^{-{\frac {1}{2}}}\,e^{-{\frac {1}{2}}(\mathbf {x} -{\boldsymbol {\mu }})^{\!{\mathsf {T}}}{\boldsymbol {\Sigma }}^{-1}(\mathbf {x} -{\boldsymbol {\mu }})}} $$
when $\boldsymbol{\Sigma}$ is positive-definite.
Does it make sense to think of a bivariate normal variable as a point in the x-y Cartesian plane, where the x coordinate comes from a Gaussian distribution ${\mathcal {N}}(\mu_1 ,\sigma_1 ^{2})$ and the y coordinate comes from a Gaussian distribution ${\mathcal {N}}(\mu_2 ,\sigma_1 ^{2})$?
I'm trying to simulate sampling 1000 times on this bivariate normal variable with the following Python code though, I'm not sure whether my understanding. Could someone please give a review?
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
x2 = np.random.normal(loc=3.0, scale=.8, size=1000)
y2 = np.random.normal(loc=2.0, scale=.9, size=1000)
plt.scatter(x2, y2)
plt.show()

You can think points sampled from a bivariate normal distribution as points in the plane but the way they are located depends on the mean vector and covariance matrix of the distribution you are considering.
In your case the procedure you have suggested works only if the two components of the random vector are independent (say uncorrelated, in the gaussian world the two are equivalent).
Suppose you have a bivariate random vector $$ \begin{bmatrix} X \\ Y \end{bmatrix} \sim \mathcal{N}\biggl{(} \begin{bmatrix} \mu_1 \\ \mu_2\end{bmatrix}, \begin{bmatrix} \sigma_{11} && \sigma_{12} \\ \sigma_{21} && \sigma_{22}\end{bmatrix} \biggr{)} $$ then if the correlation between $X$ and $Y$ is zero, both $\sigma_{12}$ and $\sigma_{21}$ are zero. In this special case we get $\Sigma$ diagonal matrix, and you can actually think to sample indpendently from $X \sim \mathcal{N}(\mu_1, \sigma_{11})$ and $Y \sim \mathcal{N}(\mu_2, \sigma_{22})$ and see the resulting point $(x,y)$ as a point whose coordinates are independently given by such sample. Geometrically speaking, if you perform such a sampling (where $\sigma_{12} \neq \sigma_{22}$) you will see the points distributed more or less as an ellipse having axis aligned with the axis of the reference plane. Having the axis aligned with the coordinate axis translates the independence between the variables of the random vector (or lack of correlation, that is, knowing something about X doesn’t tell anything about Y and viceversa). This is the case pictured in your sampling. In the even more special case of equal variance between $X$ and $Y$, meaning that $\Sigma = I$, where $I$ is the identity matrix, you will have points distributed as a circle with an opportune radius.
The reason for this is that the quadratic form $(\boldsymbol{x} - \boldsymbol{\mu})^T \Sigma (\boldsymbol{x} - \boldsymbol{\mu})$ also known as Mahalanobis distance, is the equation of an ellipse, whose shape depends on the mean vector $\boldsymbol{\mu}$ and covariance matrix $\Sigma$.
In the most general case $\Sigma$ is not diagonal, and you will have correlation between the variables $X$ and $Y$. Still you will have points placed according to an ellipse, but this time the axis of the ellipse will not be aligned with the coordinate axis. Then, if one variable increase, say X, you know that Y will increase or decrease depending on the type of correlation you have between X and Y.
In this case is better to not think to sample indipendenlty from $X$ and $Y$ but to sample from the joint distribution of the random vector $\begin{bmatrix} X \\ Y \end{bmatrix}$, otherwise you will end up discarding the possible correlation between the variables (for example if X is the weight of a person and Y the height, indipendent sampling would mean to discard the fact that taller person are probably heavier than smaller ones)