Much effort has been devoted to generalizing the Calder\'on--Zygmund theory from Euclidean spaces to metric measure spaces, or spaces of homogeneous type. Here the underlying space $\mathbb{R}^n$ with Euclidean metric and Lebesgue measure is replaced by a set $X$ with a general metric or quasi-metric and a doubling measure. Further, one can replace the Laplacian operator that underpins the Calder\'on--Zygmund theory by more general operators~$L$ satisfying heat kernel estimates. I will present recent joint work with P.~Chen, X.T.~Duong, J.~Li and L.X.~Yan along these lines. We develop the theory of product Hardy spaces $H^p_{L_1,L_2}(X_1 \times X_2)$, for $1 \leq p < \infty$, defined on products of spaces of homogeneous type, and associated to operators $L_1$, $L_2$ satisfying Davies--Gaffney estimates. This theory includes definitions of Hardy spaces via appropriate square functions, an atomic Hardy space, a Calder\'on--Zygmund decomposition, interpolation theorems, and the boundedness of a class of Marcinkiewicz-type spectral multiplier operators.