scipy.stats.multivariate_normal.cdf works on rectangular regions, giving the probability that components of a random vector would be less than components of the given vector. If we instead want the probability that a sample lies inside the ellipsoid determined by Mahalanobis distance, this can be done using chi2.cdf (see this article):

from scipy.stats import chi2
import numpy as np

mean = np.array([1, 2])
covariance = np.array([[1, 0.8],[0.8, 1]])
x = np.array([2, 3])
y = x - mean
r2 = y @ np.linalg.inv(covariance) @ y
print(chi2.cdf(r2, len(x))) # 0.4262465792625671