A general spatial interpolation method for tidal properties has been developed by solving a partial differential equation with a combination of different orders of harmonic operators using a mixed finite element method. Numerically, the equation is solved implicitly without iteration on an unstructured triangular mesh grid. The paper demonstrates the performance of the method for tidal property fields with different characteristics, boundary complexity, number of input data points, and data point distribution. The method has been successfully applied under several different tidal environments, including an idealized distribution in a square basin, coamplitude and cophase lines in the Taylor semi-infiite rotating channel, and tide coamplitude and cophase lines in the Bohai Sea and Chesapeake Bay. Compared to Laplace’s equation that NOAA/NOS currently uses for interpolation in hydrographic and oceanographic applications, the multiple-order harmonic equation method eliminates the problem of singularities at data points, and produces interpolation results with better accuracy and precision.