This is a cute problem! I toyed with it and didn't really get anywhere - I got the strong impression that it requires fields of mathematics that I am not expert in.

Indeed, given that the problem seems related to that of counting integer solutions to the equation $f(x,y) = c$, one may need to use arithmetic geometry tools (e.g. Faltings' theorem). In particular if we could reduce to the case when the genus is just 0 or 1 then presumably one could kill off the problem. (One appealing feature of this approach is that arithmetic geometry quantities such as the genus are automatically invariant (I think) with respect to invertible polynomial changes of variable such as $(x,y) \mapsto (x,y+P(x))$ or $(x,y) \mapsto (x+Q(y),y)$ and so seem to be well adapted to the problem at hand, whereas arguments based on the raw degree of the polynomial might not be.)

Of course, Faltings' theorem is ineffective, and so might not be directly usable, but perhaps some variant of it (particularly concerning the dependence on c) could be helpful. [Also, it is overkill - it controls rational solutions, and we only care here about integer ones.] This is far outside of my own area of expertise, though...

The other thing that occurred to me is that for fixed c and large x, y, one can invert the equation $f(x,y) = c$ to obtain a Puiseux series expansion for y in terms of x or vice versa (this seems related to resolution of singularities at infinity, though again I am not an expert on that topic; certainly Newton polytopes seem to be involved). In some cases (if the exponents in this series expansion are favourable) one could then use Archimedean counting arguments to show that f cannot cover all the natural numbers (this is a generalisation of the easy counting argument that shows that a 1D polynomial of degree 2 or more cannot cover a positive density set of integers), but this does not seem to work in all cases, and one may also have to use some p-adic machinery to handle the other cases. One argument against this approach though is that it does not seem to behave well with respect to invertible polynomial changes of variable, unless one works a lot with geometrical invariants.

Anyway, to summarise, it seems to me that one has to break out the arithmetic geometry and algebraic geometry tools. (Real algebraic geometry may also be needed, in order to fully exploit the positivity, though it is also possible that positivity is largely a red herring, needed to finish off the low genus case, but not necessary for high genus, except perhaps to ensure that certain key exponents are even.)

EDIT: It occurred to me that the polynomial $f(x,y)-c$ might not be irreducible, so there may be multiple components to the associated algebraic curve, each with a different genus, but presumably this is something one can deal with. Also, the geometry of this curve may degenerate for special c, but is presumably stable for "generic" c (or maybe even all but finitely many c).

It also occurs to me that one use of real algebraic geometry here is to try to express f as something like a sum of squares. If there are at least two nontrivial squares in such a representation, then f is only small when both of the square factors are small, which is a 0-dimensional set and so one may then be able to use counting arguments to conclude that one does not have enough space to cover all the natural numbers (provided that the factors are sufficiently "nonlinear"; if for instance $f(x,y)=x^2+y^2$ then the counting arguments barely fail to provide an obstruction, one has to use mod p arguments or something to finish it off...)