If I understand correctly, you have a multivariate normal in $N$ dimensions, and you are given $N – q$ coordinates. You want the multivariate normal distribution in $q$ dimensions conditioned on those coordinates, and you want to evaluate the CDF of the latter at “other” coordinates for which you know only the upper bounds.

There is no function in SciPy for this purpose, and I’m not aware of another library that has this functionality, so let’s write our own code to do the job as you suggest (get conditional MVN, evaluate CDF).

The Wikipedia article on the multivariate normal distribution has a section entitled “Conditional Distribution” (https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Conditional_distributions) the contains relevant theory. Here is code that generates an example MVN and coordinates, implements the equations to calculate the conditional probability, and performs a sanity check (assuming $N – q = 2$) using numerical integration.

import numpy as np
from scipy import stats
from scipy.integrate import dblquad
rng = np.random.default_rng(238492432)

n = 6  # dimensionality  
qc = 4  # number of given coordinates
q = n - qc  # number of other coordinates (must be 2 if you want check to work)
x = rng.random(n)  # generate values for all axes
# the first `q` are the "other" coordinates for which you want the CDF
# the rest are "given"

A = rng.random(size=(n, n))  # generate covariance matrix 
A = A + A.T + np.eye(n)*n
mu = rng.random(n)  # generate mean
dist0 = stats.multivariate_normal(mean=mu, cov=A)

# Generate MVN conditioned on x[q:] 
s11 = A[:q, :q]  # partition covariance matrix
s12 = A[:q, q:]
s21 = A[q:, :q]
s22 = A[q:, q:]
mu1 = mu[:q]  # partition mean
mu2 = mu[q:]
x1 = x[:q]  # "other" values
x2 = x[q:]  # given values

a = x2
inv_s22 = np.linalg.inv(s22)
mu_c = mu1 + s12 @ inv_s22 @ (a - mu2)
A_c = s11 - s12 @ inv_s22 @ s21
dist = stats.multivariate_normal(mean=mu_c, cov=A_c)

# Check (assumes q = 2)
def pdf(y, x):
    return dist0.pdf(np.concatenate(([x, y], x2)))

p1 = dblquad(pdf, -np.inf, x[0], -np.inf, x[1])[0]  # joint probability
p2 = dblquad(pdf, -np.inf, np.inf, -np.inf, np.inf)[0]  # marginal probability

# These should match (approximately)
dist.cdf(x1), p1/p2
# (0.25772255281364065, 0.25772256555864476)

For the other part of your question, it sounds like you want to “marginalize out” the coordinates for which you have upper limits, and you want to evaluate the PDF of the resulting distribution. However, that is a separate issue, and I think you should post a separate question for that in which you describe what you need in more detail. In the meantime, consider reviewing the “Marginal distributions” section (https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Marginal_distributions) of the Wikipedia article.