Many approaches for principal component analysis (PCA) with sparse entries are available. While sparsity facilitates interpretation, it also facilitates visualization of the complex result of multiple PCAs. Still, none of these approaches is applicable for PCA of a set of covariance matrices simultaneously. This becomes essential for spatial data, where covariance matrices for different locations can be calculated. Thus, we develop an algorithm to solve the complex optimization problem and apply it to simulations and real data examples.