Here's a suggestion:

Define a non-negative decreasing function $w(r)$ measuring interaction strength. Given each object its own independent $N(0,1)$ random variable $N_i$. Now set $$ Y_i=\frac{\sum_{j}w(\|x_i-x_j\|)N_j}{\sqrt{\sum_j w(\|x_i-x_j\|)^2}}, $$ where $x_i$ denotes the location of the $i$th object.

Then the $Y_i$ are correlated $N(0,1)$ random variables. If two objects are co-located the normal random variables agree.

Finally set $t_i=\Phi^{-1}(p_i)$ (i.e. $\mathbb P(N < t_i)=p_i$) and set $X_i=1$ if $Y_i < p_i$ and 0 otherwise.

With this setup you can write down the covariance of $Y_i$ and $Y_k$ explicitly: it's just $$ \text{Cov}(Y_i,Y_k)=\frac{\sum_j w(\|x_i-x_j\|)w(\|x_k-x_j\|)} {\sqrt{\sum_j w(\|x_i-x_j\|)^2\sum_j w(\|x_k-x_j\|)^2}}. $$

If you write this as $\cos\theta_{ik}$ then you can write the covariance of $X_i$ and $X_k$ as an integral: $$ 1/(2\pi)\int_{x < t_1,\,\cos\theta_{ik}x+\sin\theta_{ik}y < t_2}e^{-r^2/2}\,dxdy-p_ip_k. $$

Here's a suggestion:

Define a non-negative decreasing function $w(r)$ measuring interaction strength. Given each object its own independent $N(0,1)$ random variable $N_i$. Now set $$ Y_i=\frac{\sum_{j}w(\|x_i-x_j\|)N_j}{\sqrt{\sum_j w(\|x_i-x_j\|)^2}}, $$ where $x_i$ denotes the location of the $i$th object.

Then the $Y_i$ are correlated $N(0,1)$ random variables. If two objects are co-located the normal random variables agree.

Finally set $t_i=\Phi^{-1}(p_i)$ (i.e. $\mathbb P(N < t_i)=p_i$) and set $X_i=1$ if $Y_i < p_i$ and 0 otherwise.

With this setup you can write down the covariance of $Y_i$ and $Y_k$ explicitly: it's just $$ \text{Cov}(Y_i,Y_k)=\frac{\sum_j w(\|x_i-x_j\|)w(\|x_k-x_j\|)} {\sqrt{\sum_j w(\|x_i-x_j\|)^2\sum_j w(\|x_k-x_j\|)^2}}. $$

If you write this as $\cos\theta_{ik}$ then you can write the covariance of $X_i$ and $X_k$ as an integral: $$ 1/(2\pi)\int_{ x < t_1\;,\; cos\theta_{ik}x+\sin\theta_{ik}y < t_2} e^{-(x^2+y^2)/2}\,dxdy-p_ip_k. $$