The existing methods for retrieving high-spatial resolution zenith tropospheric delay (ZTD), such as interpolation and fitting methods, are not suitable for large areas and large height differences. A high-precision ZTD interpolation (HPZI) method is developed by considering the influence of height difference. The fifth-generation reanalysis dataset (ERA5) of the European Center for Medium-Range Weather Forecasts (ECMWF) data is introduced to calculate the proportion between zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD), and assumed it is equal at the collocated ERA5 and GNSS station. The ZHD and ZWD at GNSS stations are processed by considering different height correction models. Two height planes (2500 and 4400 m) are determined by averaging the heights of ERA5 and GNSS stations in Qinghai Tibet region. 96 GNSS stations in this area over the whole year of 2019 were selected, and the result shows that: the proposed HPZI method has good performance at different height ranges and different seasons, respectively, when compared with the existing polynomial interpolation and spherical harmonic interpolation methods. The averaged root mean square and mean absolute error of ZTD difference at 96 GNSS stations are 15.9/15.6 mm and 13/12.7 mm at two height planes, while those values of the interpolation method are 18.3/17.9 mm and 14.7/14.3 mm, respectively. Compared with the previous studies, the HPZI method has the highest accuracy (less than 16 mm) and considers the season influence, which was not investigated before. Such results verify the highest accuracy and robustness of the HPZI method for retrieving ZTD in large areas and large height differences.