Is the following not a counterexample, if we follow the (IMHO) standard agreement that $\deg x^iy^j=i+j$?
Let $n=4$, $q=5$. This time $\sqrt{n}-1=2-1=1$ so we are looking for a linear polynomial. So we have $f(x,y)=ax+by+c$. If we specify $f(1,1)=f(1,2)=f(2,1)=0$, then we have $ 0=a+b+c=a+2b+c=2a+b+c, $ so the only solution to these three equations is $a=b=c=0$. So if we also insist that $f(2,2)=1$, then $f$ cannot be linear.
===============================
OTOH, if we declare that $\deg x^iy^j=\max\{i,j\}$, then we can do Lagrange interpolation in two steps. Let us for each $j\in [\sqrt n]$ find a polynomial $f_j(y)\in F_q[y]$ such that $\deg f_j(y)=\sqrt n-1$, $f_j(j)=1$ and $f_j(i)=0$ for all $i\in[\sqrt n], i\neq j$. This is just usual Lagrange interpolation. We can let $ f_j(y)=K_j\,\frac{(y-1)(y-2)\cdots (y-\sqrt n)}{y-j} $ for an appropriate constant $K_j\in F_q^*$
Then as a next step let's for all $j\in[\sqrt n]$ find a polynomial $g_j(x)\in F_q[x]$ such that for all $i\in[\sqrt n]$ we have $g_j(i)=f(i,j)$. Again, Lagrange interpolation tells that such a polynomial $g_j(x)$ exists and is of degree $\le \sqrt n -1$, because its value was specified at $\sqrt n$ points.
Putting all this together we get a bivariate polynomial $ h(x,y)=\sum_{k\in[\sqrt n]} g_k(x) f_k(y)\in F_q[x,y] $ that has the prescribed value at all the points $(i,j)\in[\sqrt n]\times [\sqrt n]$. Furthermore, all the monomials $a_{ij}x^iy^j$ of $h(x,y)$ have $i,j\le\sqrt n-1$.
=============================
Edit: At some point I had a suspicion that a required polynomial $h(x,y)$ extending the given function $f$ at the given set $S$ of $n$ points on the plane could be found within the span of any $n$ monomials $x^iy^j, i,j. A moment's reflection shows that as the $x$-coordinates of the points of $S$ are constrained to a set of size $\sqrt n$, we can write the restriction of any higher power $x^i$ function to the set $S$ to lower degree terms. The recipe is simply to replace $x^i$ with its remainder when divided by the polynomial $(x-1)(x-2)\cdots (x-\sqrt n)$ that vanishes on all of $S$. Similarly for the $y$ variable.
OTOH we clearly need an $n$-dimensional space of polynomials to extend all the functions with prescribed values on a set of $n$ points. Therefore to be able to extend all the functions we are forced to use all the polynomials in the linear span of $\{x^iy^j\mid 0\le i,j <\sqrt n\}$. The restriction to functions with the property $f(S)\subseteq \{0,1\}$ does not change this fact, because within that set we can find an obvious set of $n$ linearly independent functions.