With technological advances, more functional data with high dimensionality are available in various fields such as neuroimaging analysis. Due to infinite dimensionality, functional principal component analysis is an important tool for dimension reduction, which however is scarcely researched in high dimensions. We propose sparse principal component analysis for high dimensional functional data based on the relationship between orthonormal basis expansions and multivariate K-L representations. Two sparsity regimes of interest are investigated with theoretical guarantees for the resulting estimators. Simulation and real data examples are provided to lend empirical support to the proposed method, which also performs well in subsequent analysis such as classification.