Estimation of covariance matrix function is fundamentally important for many problems. In these kinds of problems, the number of available observations at each area is often small relative to the high dimensional nature of the covariance matrix function, and not large enough to give an non-degenerate estimate. Here, we propose a local shrinkage covariance method to estimate covariance matrix function. The performance of the new proposal is compared to some non-shrinkage nonparametric methods by simulations under a wide range of scenarios. We also apply the proposal to a real MEG neuroimaging dataset, showing that the method is promising.