High-dimensional multivariate spatiotemporal data arise more and more frequently in a wide range of applications; however, there are relatively few statistical methods that can be readily used to deal with large dataset while incorporating spatial, temporal and variable dependencies simultaneously. In this paper, we propose a new approach to utilize the dependence structures in the variable, space and time dimensions to achieve dimension reduction and to facilitate spatial/temporal predictions in a high-dimensional multivariate spatiotemporal setting. Simultaneously dealing with variate, spatial and temporal covariances is made possible by modeling the multivariate spatiotemporal data as a time series of matrices whose rows and columns correspond to sampling sites and variables, respectively. We proposed methods of spatial and temporal prediction based on the estimated latent structure. Asymptotic properties of the proposed methods are established. Performance of our method is investigated on both synthetic and real datasets.