Global principal component analysis (PCA) has been successfully introduced for modeling distributed parameter systems (DPSs). In spite of the merits, this method is not feasible due to parameter variations and multiple operating domains. A novel multimode spatiotemporal modeling method based on the locally weighted PCA (LW-PCA) method is developed for large-scale highly nonlinear DPSs with parameter variations, by separating the original dataset into tractable subsets. This method implements the decomposition by making full use of the dependence among subset densities. First, the spatiotemporal snapshots are divided into multiple different Gaussian components by using a finite Gaussian mixture model (FGMM). Once the components are derived, a Bayesian inference strategy is then applied to calculate the posterior probabilities of each spatiotemporal snapshot belonging to each component, which will be regarded as the local weights of the LW-PCA method. Second, LW-PCA is adopted to calculate each locally weighted snapshot matrix, and the corresponding local spatial basis functions (SBFs) can be generated by the PCA method. Third, all the local temporal models are estimated using the extreme learning machine (ELM). Thus, the local spatiotemporal models can be produced with local SBFs and corresponding temporal model. Finally, the original system can be approximated using the sum form of each local spatiotemporal model. Unlike global PCA, which uses global SBFs to construct a global spatiotemporal model, LW-PCA approximates the original system by multiple local reduced SBFs. Numerical simulations verify the effectiveness of the developed multimode spatiotemporal model.