Remote-sensing instruments have enabled the collection of big spatial data over large domains such as entire continents or the globe. This trend is particularly apparent in large-scale natural environmental monitoring, which has commonly multivariate data. Statistical inference and prediction is generally prohibitive when the data size goes very large. Basis-function representations are well suited to such big spatial data, as they can enable fast computations for large datasets and they provide flexibility to deal with the complicated dependence structures. We will discuss a multi-resolution approximation (MRA) that uses basis functions at multiple resolutions to achieve fast inference and that can (approximately) represent any covariance structure. We present the block MRA for multivariate data, which is based on a multi-resolution partitioning of the spatial domain and a latent dimension for processes. It can deal with massive datasets with multiple variables and provides great flexibility in capturing the within and cross location dependence among processes.