A practical problem of interest in remote sensing is to increase the spatial resolution of a coarse spatial resolution image by fusing the information of that image with another fine spatial resolution image (from the same sensor or from sensors on different satellites). Thus, the problem is how to introduce spatial ‘detail’ into a coarse spatial resolution image (decrease the pixel size) such that it is coherent with the spectral information of the image. Cokriging provides a geostatistical solution to the problem and has several interesting advantages: it is a sound statistical method by being unbiased and minimizing a prediction variance (c.f. ad hoc procedures), it takes into account the effect of pixel size, and also autocorrelation in each image as well as the cross-correlation between images, it may be extended to incorporate extra information from other sources and it provides an estimation of the uncertainty of the final predictions. When formulating the cokriging system, semivariograms and cross-semivariograms (or covariances and cross-covariances) appear, some of which cannot be estimated from data directly. Cross-variograms between different variables as well as cross-semivariograms between different supports for the same variable are required. The problem is solved by using linear systems theory in which any variable for any pixel size is seen as the output of a linear system when the input is the same variable on a point support. In remote-sensing applications, the linear system is specified by the point-spread function (or impulse response) of the sensor. Linear systems theory provides the theoretical relations between the different semivariograms and cross-semivariograms. Overall, one must ensure that the whole set of covariances and cross-covariances is positive-definite and models must be estimated for non-observed semivariograms and cross-semivariograms. The models must also be realistic, taking into account, for example, the parabolic behaviour close to the origin presented in regularized semivariograms and cross-semivariograms. The solution proposed is to find by numerical deconvolution a positive-definite set of point covariances and cross-covariances and then any required model may be obtained by numerical convolution of the corresponding point model. The first step implies several numerical deconvolutions where some model parameters are fixed, while others are estimated using the available experimental semivariograms and cross-semivariograms, and some goodness-of-fit measure. The details of the proposed procedure are presented and illustrated with an example from remote sensing.