We study a stochastic heat equation with piecewise constant diffusivity $θ$ having a jump at a hypersurface $Γ$ that splits the underlying space $[0,1]^d$, $d\geq2,$ into two disjoint sets $Λ_-\cupΛ_+.$ Based on multiple spatially localized measurement observations on a regular $δ$-grid of $[0,1]^d$, we propose a joint M-estimator for the diffusivity values and the set $Λ_+$ that is inspired by statistical image reconstruction methods. We study convergence of the domain estimator $\hatΛ_+$ in the vanishing resolution level regime $δ\to 0$ and with respect to the expected symmetric difference pseudometric. As a first main finding we give a characterization of the convergence rate for $\hatΛ_+$ in terms of the complexity of $Γ$ measured by the number of intersecting hypercubes from the regular $δ$-grid. Furthermore, for the special case of domains $Λ_+$ that are built from hypercubes from the $δ$-grid, we demonstrate that perfect identification with overwhelming probability is possible with a slight modification of the estimation approach. Implications of our general results are discussed under two specific structural assumptions on $Λ_+$. For a $β$-Hölder smooth boundary fragment $Γ$, the set $Λ_+$ is estimated with rate $δ^β$. If we assume $Λ_+$ to be convex, we obtain a $δ$-rate. While our approach only aims at optimal domain estimation rates, we also demonstrate consistency of our diffusivity estimators, which is strengthened to a CLT at minimax optimal rate for sets $Λ_+$ anchored on the $δ$-grid.