The discovery of underlying surface partial differential equation (PDE) from observational data has significant implications across various fields, bridging the gap between theory and observation, enhancing our understanding of complex systems, and providing valuable tools and insights for applications. In this paper, we propose a novel approach, termed physical-informed sparse optimization (PIS), for learning surface PDEs. Our approach incorporates both $L_2$ physical-informed model loss and $L_1$ regularization penalty terms in the loss function, enabling the identification of specific physical terms within the surface PDEs. The unknown function and the differential operators on surfaces are approximated by some extrinsic meshless methods. We provide practical demonstrations of the algorithms including linear and nonlinear systems. The numerical experiments on spheres and various other surfaces demonstrate the effectiveness of the proposed approach in simultaneously achieving precise solution prediction and identification of unknown PDEs.