It is a direct consequence of the usual (weak) spatial Markov property enjoyed by the Ising model. Namely, if $\Lambda$ is a (deterministic) finite subset of $\mathbb{Z}^d$, and $f$ a local function with support inside $\Lambda$, then the expected value of $f$ in the box $\Lambda$, with a given frozen configuration $\omega$ outside $\Lambda$ only depends on the value taken by $\omega$ along the boundary of $\Lambda$ (assuming nearest-neighbor interactions, as in the linked paper). This follows immediately from the DLR equation.
The strong Markov property is the analogous property when $\Lambda$ is a random finite subset of $\mathbb{Z}^d$. In order for the property to be satisfied in this case, it is necessary that knowing $\Lambda$ only imposes conditions on the spins outside $\Lambda$ (otherwise that would have an effect on the expected value of $f$); this is the reason for asking that $\Lambda$ be determined from outside.
An example would be: "Let $\Lambda$ be the outermost circuit of $+$ spins surrounding the support of $f$ in a region $\Delta$ of $\mathbb{Z}^d$". This imposes constraints outside $\Lambda$ (since, e.g., there cannot be a circuit of $+$ spins surrounding $\Lambda$ inside the region $\Delta$), but not inside.
The proof of the strong version is then a corollary of the (weak) Markov property: just decompose the expectation of $f$ according to the realization of the random set $\Lambda$, and, for each realization, apply the weak Markov property.
PS: Forgive the plug, but I'd like to point out two recent works that provide alternative approaches to the Aizenman-Higuchi theorem, and yield considerably stronger results: arXiv:1003.6034 (Ising model) and arXiv:1205.4659 (Potts model).
EDIT: I've looked at your slightly more detailed question on math.stackexchange, but I still don't see precisely where your difficulty lies. So, let me try to give an explicit example. To keep things easy, let us only consider the nearest-neighbor model, and work in the finite box $\Lambda=\{-n,\ldots,n\}^2$ with $+$ boundary conditions. Let us consider the random variable $\sigma_0$ giving the spin at the center of the box. Given a contour $\gamma$ (I assume you know what this is) surrounding the origin, let $\mathcal{E}_\gamma$ be the event that $\gamma$ is the outermost contour surrounding the origin, and $\mathcal{E}_\varnothing$ be the event that no such contour exist. Then the expectation of $\sigma_0$ can be decomposed according to $\gamma$:
$$
\mu^+_\Lambda(\sigma_0) = \sum_\gamma \mu^+_\Lambda(\mathcal{E}_\gamma) \mu^+_\Lambda(\sigma_0 \,|\, \mathcal{E}_\gamma),
$$
where the sum is over all contours $\gamma$ in $\Lambda_n$ surrounding the origin, including the degenerate case $\gamma=\emptyset$.
The point, now, is that, when $\gamma\neq\emptyset$,
$$
\mu^+_\Lambda(\sigma_0 \,|\, \mathcal{E}_\gamma) = \mu^-_{\rm int(\gamma)}(\sigma_0),
$$
since, conditionnally on $\mathcal{E}_\gamma$, there is a path of $-$ spins along $\gamma$ (on the "inner" side of $\gamma$), which decouples the configuration inside $\gamma$ from the configuration outside. Note, and this is crucial, that the fact that we chose the outermost contour implies that the configuration inside $\rm int(\gamma)$ is unconstrained (while the configuration outside $\gamma$ is very much constrained, since it must not destroy the fact that $\gamma$ was outermost).
The last identity is the strong Markov property for the finite-volume nearest-neighbor Ising model. The property in your question is the analogous one for the infinite volume Gibbs measure.