The answer to this is yes! Let $g\colon\mathbb{R}\to\mathbb{R}$ be essentially unbounded on every nontrivial open interval and set $f(x,y)=g(x+y)$. Then, f is essentially unbounded on every measurable rectangle $A\times B$ of positive Lebesgue measure.
To construct g you can take $g_0(x)=1_{\{1>x>0\}}x^{-1/2}$ which is integrable but not square integrable on any neighborhood of 0. Then, letting $q_1,q_2,\ldots$ be an enumeration of the rationals, $g(x)=\sum_n2^{-n}g_0(x-q_n)$ is integrable (so, is finite almost everywhere) but not square integrable on any open interval, so is not essentially bounded.
To see that this construction of $f$ works, pick any two sets $A,B \subseteq \mathbb{R}$ with positive Lebesgue measure. By intersecting with a bounded interval, we can restrict to sets with finite measure. I'll write $I_A(x)$ for the indicator function of a set A and $(I_A*I_B)(x)=\int I_A(x-y)I_B(y)\,dy$ for the convolution, which is a continuous function. For any $K>0$ let $C=\{x\in\mathbb{R}\colon\vert g(x)\vert\ge K\}$, which has positive measure in any nontrivial open interval. Then,
$$
\begin{align}
\mu((A\times B)\cap \{x\in\mathbb{R}^2\colon \vert f(x)\vert\ge K\})&=\int\int I_A(x)I_B(y)I_C(x+y)\,dx\,dy\\
&=\int (I_A*I_B)(x)I_C(x)\,dx.
\end{align}
$$
The convolution $I_A*I_B$ is a nonnegative function with positive integral $\int (I_A*I_B)(x)\,dx=\mu(A)\mu(B)>0$, so it must be strictly positive at some point. By continuity, it is bounded below by a positive number on some open interval, so the integral displayed above is positive. Therefore, $f$ is essentially unbounded on $A\times B$.