I can prove the conjecture, although I can't give a bijection. I've left up a previous answer which does the case $m=3$, since it serves as an example of many of the ideas here.
Here is the key observation:
Theorem: There is a unique array $f(m,n)$, for $m$ and $n$ nonnegative integers, satisfying:
$f(m,0) = f(0,n)=1$
If we fix $m \geq 1$, then $f(m,n)$ is of the form $(m+1)^n p_m(n)$ where $p_m(x)$ is a polynomial of degree $\leq m-1$.
If we fix $n \geq 1$, then $f(m,n)$ is of the form $(n+1)^m p_n(m)$ where $p_n(x)$ is a polynomial of degree $\leq n-1$.
Uniqueness: Let's suppose that we have already found unique values of $f(m,n)$ whenever $\min(m,n) < k$. In particular, we know the values of $f(k,x)$ and $f(x,k)$ for $0 \leq x \leq k-1$. This is enough values to uniquely determine a degree $k-1$ polynomial, so the values we already know determine the values of $f(k,x)$ and $f(x,k)$ for all $x$.
Existence: The proof of uniqueness gives an algorithm to compute $f$. The only potential issue is that $f(k,k)$ is determined twice, as a polynomial in its first variable and in its second variable. However, the algorithm is symmetric in $m$ and $n$, so the two interpolating polynomials coincide. $\square$
For example, we have $f(1,0) = 1 = 2^0 \cdot 1$ so we must have $p_1(n)=1$ and $f(1,n) = 2^n$. Similarly, $f(2,0)=1=2^0 \cdot 1$ and $f(2,1) = 4=3^1 \cdot (4/3)$, so we must have $p_2(n) = 1+n/3$ and $f(2,n) = 3^n(1+n/3) = 3^{n-1} (n+3)$. Using this algorithm, I computed $f(m,n)$ for $0 \leq m,n \leq 10$. Here is the data, and some Mathematica code, if you'd like to play with it.
f[m_, n_] := f[m, n] =
If[m == 0 || n == 0, 1,
If[m <= n,
(m + 1)^n * InterpolatingPolynomial[Table[{k, f[m, k]/(m + 1)^k}, {k, 0, m - 1}], x] /. x -> n,
(n + 1)^m * InterpolatingPolynomial[Table[{k, f[k, n]/(n + 1)^k}, {k, 0, n - 1}], x] /. x -> m]]
Values of $f(m,n)$ for $0 \leq m,n \leq 10$.
1,1,1,1,1,1,1,1,1,1,1
1,2,4,8,16,32,64,128,256,512,1024
1,4,15,54,189,648,2187,7290,24057,78732,255879
1,8,54,328,1856,9984,51712,260096,1277952,6160384,29229056
1,16,189,1856,16145,129000,968125,6925000,47690625,318500000,2073828125
1,32,648,9984,129000,1475856,15450912,151201728,1403288064,12479546880,107152782336
1,64,2187,51712,968125,15450912,219682183,2862173104,34828543449,401200569280,4418300077219
1,128,7290,260096,6925000,151201728,2862173104,48658878080,760774053888,11122973573120,153936323805184
1,256,24057,1277952,47690625,1403288064,34828543449,760774053888,15047189968833,274908280855680,4707029151392121
1,512,78732,6160384,318500000,12479546880,401200569280,11122973573120,274908280855680,6199170628499200,129708290461760000
1,1024,255879,29229056,2073828125,107152782336,4418300077219,153936323805184,4707029151392121,129708290461760000,3283463201858585471
So, we can prove the conjecture by showing that both $F_{m,n}$ and $S_{m,n}$ obey conditions 1, 2 and 3.
For condition 1, we can just define $S_{m,0}$, $S_{0,n}$, $F_{m,0}$, and $F_{0,n}$ to be singletons (and I claim this is actually the most natural definition); then we just need to check that our verification of polynomiality goes all the way down to the zero case.
Since conditions 2 and 3 are symmetric, and the definitions of $F_{m,n}$ and $S_{m,n}$ are symmetric, we just check condition $2$.
Verification of condition 2 for forests: Let $A_{m,b}$ be the set of isomorphism classes of bicolored forests whose white vertices are labeled $\{1,2,\ldots,m \}$ and which have $b$ black vertices, each of degree $\geq 2$. We claim that
$$|F_{m,n}| = \sum_b |A_{m,b}| n(n-1)(n-2) \cdots (n-b+1) (m+1)^{n-b} . \qquad (1) $$
Proof: Take a forest in $F_{m,n}$, with the $m$ vertices colored white and the $n$ vertices colored black.
Delete all black vertices of degree $0$ and $1$ and forget the numbering of the remaining black vertices. This gives a forest in $A_{m,b}$. To undo this process, we first must take the $b$ black vertices and decide which of the $n$ vertices of $K_{m,n}$ they will be; there are $n(n-1)(n-2) \cdots (n-b+1)$ ways to do this.
(We used that trees in $A_{m,b}$ have no automorphisms.) Then we must take the $n-b$ remaining black vertices and either join them to one of the $m$ white vertices or leave them unconnected; there $(m+1)^{n-b}$ ways to do this.
Formula (1) is clearly $(m+1)^n$ times a polynomial, it remains to check that the polynomial has degree at most $m-1$. We must check that $A_{m,b}$ is empty if $b \geq m$. Indeed, since every black vertex has degree at least $2$, a forest in $A_{m,b}$ has at least $2b$ edges. But, since it is a forest, it also has at most $m+b-1$ vertices. So $2b \leq m+b-1$ and $b \leq m-1$. $\square$
Verification of condition 2 for margins of $(0,1)$ matrices: Consider a vector $C = (C_1, \ldots, C_n) \in \{ 0,1,\ldots, m \}^n$ of columns sums, and consider how many row sums $(R_1, \ldots, R_m)$ are compatible with it. Let $c_j = \# \{ k : C_k = j \}$, so $c_0+c_1+\cdots + c_m=n$.
By the Gale-Ryser theorem, $(R_1, \ldots, R_m)$ is compatible with $(c_0, \cdots, c_m)$ if and only if $\sum R_i = \sum j c_j$ and, for all index sets $1 \leq i_1 < i_2 < \cdots < i_k \leq m$, we have
$$R_{i_1} + R_{i_2} + \cdots + R_{i_k} \leq \sum_j \min(j,k) c_j . $$
This is the defining list of inequalities of the permutahedron, so this condition can alternately be stated as saying that $(R_1, \ldots, R_m)$ is in the convex hull of the $m!$ permutations of $(c_1+c_2+\cdots+c_m, c_2+\cdots + c_m, \cdots, c_{m-1}+c_m, c_m)$.
By Theorem 11.3 of Postnikov's Permutohedra, associahedra, and beyond, the number of such $(R_1, \ldots, R_m)$ is a polynomial $g_m(c_0,c_1, \ldots, c_m)$ in the $c_j$, of degree $m-1$. So
$$|S_{m,n}| = \sum_{c_0+\cdots + c_m=n} \frac{n!}{c_0! c_1! \cdots c_m!} \ g_m(c_0,c_1, \ldots, c_m). \qquad (2)$$
Let $\partial_j$ be $\tfrac{\partial}{\partial x_j}$.
There is some polynomial $h_m$ in the $\partial_j$, of degree $m-1$, such that
$$\left. h_m(\partial_0, \ldots, \partial_m) \left( x_0^{c_0} \cdots x_m^{c_m} \right) \right|_{(1,1,\ldots,1)} = g_m(c_1, \ldots, c_m). \qquad (3)$$
Combining (2), (3) and the binomial expansion of $(x_0+x_1+\cdots+x_m)^n$, we obtain
$$|S_{m,n}| = \left. h_m(\partial_0, \ldots, \partial_m) \left( x_0+x_1+\cdots+x_m \right)^n \right|_{(1,1,\ldots,1)} . \qquad (4) $$
Let $h_m(t,t,\ldots,t) = \sum_{k=0}^{m-1} h_{m,k} t^k$. Then $(4)$ is $\sum_{k=0}^{m-1} h_{m,k} n(n-1)(n-2) \cdots (n-k+1) (m+1)^{n-k}$, which is a polynomial of the desired form. $\square$
Best Answer
It seems from your description that you don't need to actually compute the matrix, but just insure that the necessary and sufficient conditions exist for one to be computed. For that, I think that all you'll need is J. H. Ryser's paper, "Combinatorial properties of matrices of zeros and ones." You may also find Richard A. Brualdi's paper "Matrices of zeros and ones with fixed row and column vectors" very readable.
Applying those, we first establish some terminology. As you've already said, the matrix ${\bf A}$ is an $m \times n$ (where $m$ and $n$ are positive integers) matrix of ones and zeros. Further, we have row and column sum vectors ${\bf R} = (r_1,\dots,r_m$) and ${\bf S}=(s_1,\dots,s_n)$ which are nonnegative integral vectors. First, it's clear that if there are $\tau$ ones in the matrix ${\bf A}$, then $\displaystyle \tau = \sum_{i=0}^{m}r_i = \sum_{j=0}^{n}s_j$. If that condition fails, then obviously the answer is "no."
If that succeeds, then following Brualdi, we can define a function to help with the evaluation. First, let $I \subseteq \left\{ 1,\dots,m \right\}$ and $J \subseteq \left\{1,\dots, n \right\}$. Define $\displaystyle t(I,J) =|I||J|+\sum_{i \notin I}r_i - \sum_{j \notin J}s_j$. All that remains is to demonstrate that $t(I,J) \ge 0$ for all $I\subseteq\left\{1,\dots,m\right\}$ and $J\subseteq\left\{1,\dots,n\right\}$.