The presence of clouds and their shadows in remotely sensed images limits their potential uses for extracting information. The commonly used methods for replacing clouded pixels by land cover reflection estimates usually yield poor results if the images being combined exhibit radical differences in target radiance due, for example, to large date separation and high temporal variability. This study focuses on introducing geostatistical techniques for interpolating the DN values of clouded pixels in multispectral remotely sensed images using traditional ordinary cokriging and standardized ordinary cokriging. Two case studies were conducted in this study. The first case study shows that the methods work well for the small clouds in a heterogeneous landscape even when the images being combined show high temporal variability. Although the basic spatial structure in large size clouds can be captured, image interpolation-related artefacts such as smoothing effects are visually apparent in a heterogeneous landscape. The second case study indicates that the cokriging methods work better in homogenous regions such as the dominantly agricultural areas in United States Midwest. Various statistics including both global statistics and local statistics are employed to confirm the reliability of the methods.