Correlation function at zero distance

I'm confused about the definition of the correlation function (at equal time). I know it is defined from the thermal average of the scalar product of two random variables (for example the spins of a magnetic system) as $$G(r)=\langle s_1(0) s_2(r)\rangle - \langle s_1(0)\rangle \langle s_2(r)\rangle$$ where I have used the notation of connected correlator, substracting the uncorrelated variables from the main correlator. Morover we know that the correlation function is given asymptotically by the formula $$G(r) \sim \frac{e^{-r/\xi}}{r^{d-2+\eta}}$$ where $$\xi$$ is the correlation length of the model and $$\eta$$ is a critical exponent. From this expression and from the analogy with correlators in QFT (which diverge when the fields they involve are situated at the same position), I would say that the correlation function close to zero distance diverges $$G(r)\to \infty$$ for $$r\to 0$$. But actually all the plots I have found - like this one

(By User:Mgibby5 - Own work, CC BY-SA 4.0, Link)

show that $$G(r)=1$$. I don't understand why is it possible and what's the physical interpretation for that.

The correlation function measures, well, the correlation between two random events. Imagine $$s_1$$ and $$s_2$$ are random variables measuring the outcome of two coin flips - $$+1$$ if the coin shows up heads, and $$-1$$ if the coin shows up tails. If the two coins are completely uncorrelated - the outcome of one doesn't affect the outcome of the other - then you can convince yourself that the correlation function $$G = \langle s_1 s_2 \rangle - \langle s_1 \rangle \langle s_2 \rangle = 0$$. Likewise, if the two coins are perfectly correlated - if one shows up as heads, the other will always show up as heads, and vice versa - then $$G = 1$$. Random variables are perfectly correlated with themselves (if coin #1 is heads, then coin #2 will also be heads if they're the same coin!), so we would expect $$G=1$$ in this case.
As an aside, you'll notice that the fact that $$G=1$$ depends on my choice to represent heads/tails as $$\pm 1$$, respectively. If I had chosen $$\pm 2$$, or if I had chosen tails to be represented by $$0$$, then perfect correlation would not cause $$G=1$$. However, in the spin case, up/down spins are represented by $$\pm 1$$, so there isn't an ambiguity here. It can be useful to normalize your variables such that perfect correlation implies $$G=1$$.
Back to your question. In the derivation of the correlation function at the critical point, it's implicitly assumed during the derivation that $$r$$ is well-separated from $$0$$. This is because the derivation takes the distance between two spins, $$r$$, to be a continuous variable. Normally, one defines an Ising model on a discrete lattice, with a countable number of sites. If one only cares about long distance behavior, they can ignore the discreteness of the lattice, and "pretend" that the Ising spin is a continuous variable defined at every point in space. Upon making this assumption, they can derive the expression for $$G(r)$$ cited above. If instead one was to not make this assumption, and really take it seriously that there are only a countable number of Ising spins on a lattice, one can instead derive a form for the correlation function that approaches $$1$$ as the separation between the two spins goes to zero (you can look up things like the "transfer matrix" method for examples of this. My computer is about to run out of battery, but I'll try to add a link later if I get a chance.)