Existing potential evapotranspiration (PET) retrieval methods can obtain high-precision PET values at specific stations but not well used, which becomes the focus of this paper. A multi-source PET fusion (MPF) method using station- and grid-based PET is proposed to obtain surface-based high-precision PET. The grid-based empirical PET (GEP) model is initially established by analyzing the relationship between precipitable water vapor/temperature and the PET difference between Penman–Monteith and Thornthwaite model. Subsequently, the MPF model is developed by fusing the station- and grid-based PET. In addition, the optimal weights of station- and grid-based PET for MPF model are determined by introducing the Helmert variance component estimation method. Loess Plateau (LP) area is selected to perform the experiment. The corresponding data of 33 global navigation satellite system stations, 84 meteorological stations, and grid-based points with spatial resolution of 0.25° × 0.25° provided by ERA-Interim are used over the period of 2012–2018. The statistical result shows that the average root mean square (RMS) of MPF-derived PET at 84 stations is 6.53 mm/month in LP area. In addition, the RMS improvement rate (IR) of MPF-derived PET is 80.52% when compared with that of the TH model. Comparisons of long-time-series MPF-derived PET, standardized precipitation evapotranspiration index (SPEI), precipitation smoothing index (PSI), standardized precipitation index (SPI), and standardized precipitation conversion index (SPCI) also show the good performance of proposed MPF method. The RMS IRs of SPEI calculated using the MPF-derived PET are 64.3%, 69.0%, 68.1%, and 8.2%, respectively, when compared with the TH-derived SPEI at 1-, 3-, 6-, and 12-month scales, respectively. Such results verified the effectiveness and robustness of the proposed MPF method for obtaining surface-based PET and SPEI values.