Kernel methods are an efficient tool to build surrogate models from scattered data. In this talk, we consider scattered data approximation by kernel methods on product regions of equal and different dimensionality. On each of these regions, we assume quasi-uniform but unstructured data sites. We introduce multilevel hierarchies and construct optimal sparse grids for scattered data approximation on the product region. For this, we derive new improved error estimates for the respective kernel interpolation error by invoking duality arguments. An efficient algorithm to solve the underlying linear system of equations is proposed. The algorithm is based on the sparse grid combination technique, where a sparse direct solver is used for the elementary anisotropic tensor product kernel interpolation problems. The application of the sparse direct solver is facilitated by applying a samplet matrix compression to each univariate kernel matrix, resulting in an essentially sparse representation of the latter. In this way, we obtain a method that is able to deal with large problems up to billions of interpolation points, especially in case of reproducing kernels of nonlocal nature. Numerical results are presented to qualify and quantify the approach.