We consider the use of kernel methods for the extraction of spatiotemporal patterns from PDE simulation data. We begin with a brief overview of the formulation of pattern extraction as the eigenvalue problem for a kernel integral operator. We then consider spatiotemporal correlation kernels that act on the product space of temporal snapshots and spatial gridpoints. The introduced family of kernels leads to kernel eigenfunctions that are invariant under the spatial symmetries of the dynamics and inherently spatiotemporal. As a result, each one of them is capable of encoding complex and dynamically relevant patterns. We pay specific attention to the scalable implementation of the presented method via low rank approximation and distributed GPU computing. We conclude with numerical results for the modified Hasegawa-Wakatani equations of plasma dynamics.