Determining covariance of irregularly spaced spatial data I'm comparing concentration $C$ of a contaminant in the same spatial region at two time point 2000 and 2010 with sample size of $N_{2000}$ = 51 and $N_{2010}$ = 26 (not all the samples are from the same location), mean of $\mu(C)_{2000}$ = 47 and $\mu(C)_{2010}$ = 27 (determined by block kriging of all point observations) and variance of $V(C)_{2000}$ = 89 and $V(C)_{2010}$ = 68 (kriging variance). To determine if there has been any significant change over the last 10 years, we first need to determine the variance of change in the area:
$V(\Delta C) = V(C)_{2000} + V(C)_{2010} - V(C)_{2000,2010}$
where, $V(\Delta C)$ is the variance of change over time; and $V(C)_{2000,2010}$ is the covariance between the two temporal samples. Does anyone know how to determine the $V(C)_{2000,2010}$ term in the above equation?
 A: If

not all the samples are from the same location

is equivalent to

almost all the samples are from the same location

or put differently that $2000$ and $2010$ spatial supports broadly intersect (as in $\rm{Fig.1}$), the approach below is applyable.

What "Broadly" means ? Whether focusing only on the intersection does or does not change your question of reasearch. E.g. if the intersection is restricted to a too small, say, urban area, and that you were originally interested in the overall metropolitan area, your question of reasearch would be changed and what follows will thus not suit you.

Focusing only on the intersection of the two years spatial support
You can build a $2D$ kriging-interpolated continuum on which you can then project a grid whose peremiter consists of the convex envelope of the set  of points belonging to the intersection described above (and shown in $\rm{Fig.1}$). The nodes of the so-projected grid are going to be used as "individuals" within the bootstrapping proccess. Like so (see $\rm{Fig.2}$), each node will have year-$2000$ and year-$2010$ concentrations attached.

The steps are:


*

*Get, for each boostrap resample $b=1, ..., k$ and for all of its individuals (nodes) $i=1, ..., n$ (duplicated positions $\iff$ not all positions from the original sample, since the resample size is still $n$): $C_{b,i,2000}$, $C_{b,i,2010}$ and  $C_{b,i,2000} \times C_{b,i,2010}$

*Compute the three sample means of each resample $b=1, ..., k$ as follows
$\forall b, \overline{C_{b,2000}} = \frac{1}{n}\sum_{i=1}^n C_{b,i,2000}$
$\forall b, \overline{C_{b,2010}} = \frac{1}{n}\sum_{i=1}^n C_{b,i,2010}$
$\forall b, \overline{C_{b,2000} \times C_{b,2010}}= \frac{1}{n}\sum_{i=1}^n C_{b,i,2000} \times C_{b,i,2010}$
Once you are provided with these three empirical bootstrap distributions of sample means, you may want to compute the corresponding three empirical bootstrap means


*Compute $\hat{E}(C_{2000})$, $\hat{E}(C_{2010})$ and $\hat{E}(C_{2000} \times C_{2010})$, as follows:
$\hat{E}(C_{2000}) = \frac{1}{k} \sum_{b=1}^k \overline{C_{b,2000}}$
$\hat{E}(C_{2010}) = \frac{1}{k} \sum_{b=1}^k \overline{C_{b,2010}}$
$\hat{E}(C_{2000} \times C_{2010}) = \frac{1}{k} \sum_{b=1}^k \overline{C_{b,2000} \times C_{b,2010}}$
And finally (also reusing your notation for covariance)


*Compute $V(C)_{2000,2010} = cov(C_{2000},C_{2010}) = \hat{E}(C_{2000} \times C_{2010}) - \hat{E}(C_{2000})\hat{E}(C_{2010})$


And do not forget to set $k$ as close as possible to $\frac{(2n-1)!}{n!(n-1)!}$ (which is very likely to be computationally expensive).

Naturally, I could have skipped step 2, computing directly the three empirical bootstrap means, e.g. $\hat{E}(C_{2000}) = \frac{1}{kn} \sum_{b=1}^k \sum_{i=1}^nC_{b,i,2000}$

The critical points in this approach are (i) how such it is possible for the practician to bootstrap over his grid nodes and (ii) how such she can access data related to each grid node and process them (computing and storing results)

My sources are theoretical: using the two definitions of bootstrapping and covariance.
