Key challenges in the analysis of highly multivariate large-scale spatial stochastic processes, where both the number of components (p) and spatial locations (n) can be large, include achieving maximal sparsity in the joint precision matrix, ensuring efficient computational cost for its generation, accommodating asymmetric cross-covariance in the joint covariance matrix, and delivering scientific interpretability. We propose a cross-MRF model class, consisting of a mixed spatial graphical model framework and cross-MRF theory, to collectively address these challenges in one unified framework across two modelling stages. The first stage exploits scientifically informed conditional independence (CI) among p component fields and allows for a step-wise parallel generation of joint covariance and precision matrix, enabling a simultaneous accommodation of asymmetric cross-covariance in joint covariance matrix and sparsity in joint precision matrix. The second stage extends the first-stage CI to doubly CI among both p and n and unearths the cross-MRF via an extended Hammersley-Clifford theorem for multivariate spatial stochastic processes. This results in the sparsest possible representation of the joint precision matrix and ensures its lowest generation complexity. We demonstrate with 1D simulated comparative studies and 2D real-world data.