What is the difference between scale invariance and self-similarity?

From the point of view of non-linear dynamics where self-similarity plays an important role if the attractor is a fractal I would say that the difference is one between continuous and discrete transformations.

A self-similar transformation like the one producing the Cantor set or the Sierpinski triangle proceeds by discrete stages. The fractal which is the limit when the number of stages $N$ tends to infinity shows self-similarity (i.e. is identical to itself) only for a discrete number of stages.

For instance when zooming on the Sierpinski triangle, one may not zoom anywhere and by any zooming factor. One has to zoom only with a factor $1/3$ and center the zoom on the symmetry axis of the triangle. So basically the number of self-similar objects is an integer and has for characteristic the self-similarity dimension which is a number $D$ such as $N = L^D$ where $N$ is the number of copies produced by changing the size by $L$.

As for scale invariance which is not so largely used, it is a statement that $f(\mu x) = \mu^D \, f(x)$ with some constant $D$. The property is continuous and true for every $x$. Fractal attractors are generally not exactly scale invariant - they have often 2 or several different scalings.

Hence from this point of view the self-similarity and scale invariance may only be identical in a discrete number of points for simple fractals which have a unique scaling factor. (I am aware that this does not adress spin lattices but it answers the question in the frame of the chaos theory.)

The discussion p. 24 on renormalons by Rosten does not look to me related to the issue of scale invariance at all. As for the one on pp. 8-9, I think this is just a remark in passing about typical spin configurations at the critical point. Let me give precise definitions to explain what is going on in the language of probability theory.

Take say the nearest neighbor Ising model in dimension $d\ge 2$. This is a collection of probability measures $\mu_{\beta,h}$ on the space $\Omega=\{-1,1\}^{\mathbb{Z}^d}$ which depends on the two parameters given by the inverse temperature $\beta\in [0,\infty)$ and the magnetic field $h\in\mathbb{R}$. The critical case to which the "false picture" remark by Rosten applies to is the particular case $\beta=\beta_{\rm c}$ and $h=0$. The remark pertains to some geometric property $P(\sigma)$ of a spin configuration $\sigma=(\sigma_{\mathbf{x}})_{\mathbf{x}\in\mathbb{Z}^{d}}\in\Omega$. Here the property $P(\sigma)$ is something like "the configuration $\sigma$ has all these nested oceans/islands of up and down spins". The belief that Rosten disputes (following Section 1.4.2 of Delamotte as Kostya rightly pointed out) is: $P(\sigma)$ is true $\mu_{\beta_{\rm c},0}$-almost surely. This belongs to the chapter on the Strong Law of Large Numbers in probability theory.

Now the property $P(\sigma)$ definitely has a "self-similar" flavor to it but I would not at all call that scale invariance. The latter has to do with a different chapter, namely, the one on the Central Limit Theorem.

Let me just work in the single phase region $\beta\le \beta_{\rm c}$. Then the general picture is that for any $\beta, h$ you can find a quantity $[\phi]_{\beta,h}$ for which the following happens. Pick some fixed integer $L>1$. For each integer $n\ge0$ define a probability measure $\nu_{\beta,h,L,n}$ on the space of temperate distributions $S'(\mathbb{R}^d)$ as the law of the random distribution $\phi_n(x)$ given by $$ \phi_{n}(x)=L^{-n(d-[\phi]_{\beta,h})}\sum_{\mathbf{y}\in\mathbb{Z}^d} (\sigma_{\mathbb{y}}-\langle \sigma_{\mathbb{y}}\rangle)\ \delta^d(x-L^{-n}\mathbf{y}) $$ with $\sigma$ sampled according to the lattice measure $\mu_{\beta,h}$, and where the magnetization $\langle \sigma_{\mathbb{y}}\rangle$ is the expectation of a single spin for that measure. Then the $\phi_n$ should converge in (probability) distribution to a random (Schwartz) distribution with law $\nu_{\beta,h,L,\infty}$. Just from the existence and uniqueness of this limit, it follows that the limit law or scaling limit $\nu_{\beta,h,L,\infty}$ is scale invariant with respect to the group $L^{\mathbb{Z}}$, i.e., rescaling by powers of $L$ (also the linear size of a block in a block-spin RG procedure). This is a discrete scale invariance property (related to Stan's answer). In fact here one expects this to hold for any $L$ and thus $\nu_{\beta,h,L,\infty}$ should be invariant by the full multiplicative group $(0,\infty)$. In other words it should have the continuous scale invariance property. BTW, the difference between discrete and continuous scale invariances is discussed for instance in this review by Sornette. In our situation continuous scale invariance means in particular that say the two point function satisfies $$ \langle\phi(\lambda x)\phi(\lambda y)\rangle=\lambda^{-2[\phi]_{\beta,h}}\langle\phi(x)\phi(y)\rangle $$ for all $\lambda\in (0,\infty)$.

Now the key question is: how to pick the crucial quantity $[\phi]_{\beta,h}$ (the scaling dimension of the spin field) needed to see the scale invariance?

For $d=2$, one has $[\phi]_{\beta_{\rm c},0}=\frac{1}{8}$ and the limit $\nu$ is non-Gaussian: the $m=3$ unitary minimal conformal field theory.

For $d\ge 4$, one has $[\phi]_{\beta_{\rm c},0}=\frac{d-2}{2}$ and the limit is the massless Gaussian field.

The most interesting situation is $d=3$ where according to best present estimates $[\phi]_{\beta_{\rm c},0}\simeq 0.5181489$ (see this article).

Now in any dimension and when $(\beta,h)\neq (\beta_{\rm c},0)$ one has $[\phi]_{\beta_{\rm c},0}=\frac{d}{2}$. In this case the limit is called Gaussian white noise. The two point function should be homogeneous of degree $-d$. This might suggest that it should be given by $$ \langle\phi(x)\phi(y)\rangle\sim \frac{1}{|x-y|^d} $$ but that is wrong. The unique rotationally invariant element of $S'(\mathbb{R}^d)$ (up to multiplicative constant) with homogeneity of degree $-d$ is $\delta^d(x)$. One is really in a central limit regime where one divides a sum of $N=L^{nd}$ real random variables by $\sqrt{N}$. Finally, if instead of the right guess for $[\phi]_{\beta,h}$ one uses some number $\alpha>0$ and writes $$ \phi_{n}(x)=L^{-n(d-\alpha)}\sum_{\mathbf{y}\in\mathbb{Z}^d} (\sigma_{\mathbb{y}}-\langle \sigma_{\mathbb{y}}\rangle)\ \delta^d(x-L^{-n}\mathbf{y}) $$ then the following happens. If $\alpha<[\phi]_{\beta,h}$, there is too much suppression and the $\phi_n$ converge in distribution to the (non random) identically zero field. If $\alpha>[\phi]_{\beta,h}$ there is loss of tightness, i.e., probability mass escaping to infinity and the limit of the $\phi_n$ does not exist. Again, this is like the central limit theorem if you don't divide by the correct power of $N$, namely, $\sqrt{N}$.

BTW, pretty much everything I said above is substantiated by rigorous mathematical proofs. For instance, the case $(\beta,h)\neq(\beta_{\rm c},0)$ with convergence to white noise follows from a general result by Newman for FKG systems. For $d=2$ and $(\beta,h)=(\beta_{\rm c},0)$, this is in the recent work of Chelkak, Hongler and Izyurov, Dubédat, as well as Garban and Newman. Of course, the main open case is $d=3$ and $(\beta,h)=(\beta_{\rm c},0)$, i.e., the infamous 3D Ising CFT. As far as I know, the best rigorous result so far (see these recent lectures by Duminil-Copin) is that at $(\beta,h)=(\beta_{\rm c},0)$ the two-point function satisfies $$ \frac{c_2}{|\mathbf{x}-\mathbf{y}|^{2}}\le \langle \sigma_{\mathbf{x}}\sigma_{\mathbf{y}}\rangle\le \frac{c_1}{|\mathbf{x}-\mathbf{y}|} $$ for come positive constants $c_1,c_2$, when the actual conjecture is more like $$ \langle \sigma_{\mathbf{x}}\sigma_{\mathbf{y}}\rangle \sim \frac{1}{|\mathbf{x}-\mathbf{y}|^{1.0362978}}\ . $$