<p>Ocean observations are inherently characterized by irregular temporal and spatial distributions, as well as heterogeneous spatial resolutions and error characteristics arising from the use of diverse observational platforms and techniques. To enable their application across a broad range of scientific and practical problems, it is essential to map these heterogeneous datasets into temporally and spatially consistent gridded products. Optimal Interpolation remains the most widely adopted algorithm for the mapping of oceanographic data. Two principal implementations of the optimal interpolation algorithm are commonly employed. The first, known as the basic optimal interpolation, is derived from the theory of optimal estimation and involves computationally intensive matrix operations, posing significant challenges when applied to high-dimensional problems. The second, referred to as the point-wise optimal interpolation, reduces computational complexity through point-wise estimation, thereby circumventing high-dimensional operations; however, this approach results in a substantially higher overall computational cost. In this study, a novel optimal interpolation algorithm is proposed that utilizes the Kronecker product to approximate the background error covariance matrix. This formulation enables the decomposition of high-dimensional matrix operations into smaller, computationally tractable sub-problems, thereby improving the scalability of optimal interpolation for large spatial domains with dense observational coverage. Building upon this framework, a multi-scale optimal interpolation method is further developed to enhance the integration of observational datasets with widely varying spatial resolutions, thereby improving the accuracy and applicability of the resulting gridded products.</p>