Call a $2\times 2$ matrix with coefficients in $\{0,1,2,3,\ldots\}$ positively reduced if any row or column reduction (given by replacing a row/column by itself minus the other row/column) produces at least one strictly negative entry.
A matrix $\left(\begin{array}{cc}a&b\\c&d\end{array}\right)$ (with coefficients $a,b,c,d\in\{0,1,2\ldots\}$) of strictly positive determinant $n=ad-bc>0$ is positively reduced if and only if $0\leq \max(b,c)<\min(a,d)$. One shows easily that $\max(a,d)\leq n$. The number of positively reduced matrices of given strictly positive determinant is thus finite.
(Observation: Positive reduction can be seen as a generalization of Euclids algorithm which is positive reduction on (columns of) a $1\times 2$ matrix.)
We denote by $\mathcal R(n)$ the finite set of all positively reduced matrices of determinant $n>0$ and by $r(n)=\sharp(\mathcal R(n))$ the number of elements in $\mathcal R(n)$. (The set $\mathcal R(0)$ of positively reduced matrices with determinant $0$ consists of all matrices with at most one strictly positive entry and is infinite. The set $\mathcal R(-n)$ is isomorphic to $\mathcal R(n)$ by row-exchange.)
The first $25$ values of the sequence $r(1),r(2),\ldots$ are given by $$1,2,3,5,5,8,7,11,10,14,11,19,13,20,18,24,17,30,19,31,26,26,32,23,44,26$$ and are not recognized by https://oeis.org/.
The numbers $r(n)$ and the set $\mathcal R(n)$ seem to have interesting arithmetical properties:
We have seemingly $r(n)\geq n$ with equality if and only if $n=1$ or $n$ is prime.
More precisely, for a prime $p$, (reductions modulo $p$ of) elements in $\mathcal R(p)$ have all rank $1$ over the finite field $\mathbb F_p$. The two row-vectors of a matrix in $\mathcal R(p)$ are thus collinear over $\mathbb F_p$ and represent a unique element $[a:b]$ of the projective line $\mathbb P^1(\mathbb F_p)$ over $\mathbb F_p$. This induces a bijection between $\mathcal R(p)$ and $\mathbb P^1(\mathbb F_p)\setminus\{[1:1]\}$. (See update below for a proof.)
The two points $[1:0]$ and $[0:1]$ for example correspond to the two diagonal matrices $\left(\begin{array}{cc}1&0\\0&p\end{array}\right)$ and $\left(\begin{array}{cc}p&0\\0&1\end{array}\right)$.
Update: The definition of positive reducedness can be reformulated as follows. The two rows $v,w$ of a positively reduced matrix $M$ of determinant $n$ generate a sublattice $\Lambda(M)=\mathbb Zv+\mathbb Zw$ of index $n$ in $\mathbb Z^2$. The map $\mathcal R(n)\ni M\longmapsto \Lambda(M)$ is injective but never surjective for $n>1$. More precisely, a sublattice $\Lambda$ of index $n$ in $\mathbb Z^2$ has an associated (unbounded) "Newton polygon" $P(\Lambda)$ by considering the convex hull of $\Lambda\cap(\mathbb N^2\setminus\{(0,0)\})$.
A sublattice $\Lambda$ of index $n$ in $\mathbb Z^2$ is the of the form $\Lambda(M)$ for $M$ in $\mathcal R(n)$ if and only if the boundary $\partial P(\Lambda)$ of the polygone $P(\Lambda)$ intersects the diagonal $\mathbb R(1,1)$ in a point $Q$ which is not integral. The two rows of $M$ are then given by the two consecutive integral points on $\partial P(\Lambda)$ separated by $Q$.
This implies the result easily for $n=p$ prime: The only sublattice of prime index $p$ in $\mathbb Z^2$ with $\partial P(\Lambda)$ intersecting $\mathbb N(1,1)$ is the lattice $\mathbb Z(p,0)+\mathbb Z(1,1)$.
This should give a fast algorithm by removing the number of "bad" lattices (with $\partial P(\Lambda)$ intersecting $\mathbb Z(1,1)$) from the total number of sublattices of index $n$ in $\mathbb Z^2$. (The number of sublattices of given index $n$ in $\mathbb Z^d$ is given by $$\prod_{p\vert n}{e_p+d-1\choose d-1}_p$$ where $n$ has prime-decomposition $n=\prod_{p\vert N}p^{e_p}$ and where ${a\choose b}_q$ is the $q$-binomial, see Y. M. Zou, “Gaussian binomials and the number of sublattices”, Acta Crystallogr. 62 (2006), no. 5, p. 409--410). (End of update)
More generally, the function $n\longmapsto r(n)$ behaves almost multiplicatively:
We have experimentally for small prime-powers: \begin{align*} r(p^2)&=p^2+1,\\ r(p^3)&=p^3+p^2-p+1,\\ r(p^4)&=p^4+p^3-p+2,\\ r(p^5)&=p^5+p^4+p^3-p^2-p+2,\\ r(p^6)&=p^6+p^5+p^4-p^2-p+3. \end{align*}
$r(pq)=(p+1)(q-1)+2$ for primes $p<q$,
$r(p^kq)=\frac{p^{k+1}-1}{p-1}(q-1)+k+1$ for primes $q>p^k$ (checked for small cases),
for three distinct primes $p_1,p_2,p_3>p_1p_2$ we have $r(p_1p_2p_3)=(p_1+1)(p_2+1)(p_3-1)+4$,
Question: Does someone have any ideas for proofs? The function $n\longmapsto r(n)$ seems to have interesting arithmetic properties: Is there something similar to a modular form lurking somewhere? Ideas for fast algorithms (for computing $r(n)$ or $\mathcal R(n)$) would also be wellcome. (I think my algorithm below is worse than $n^2$.)
(Update on questions: The above update transforms the questions into the determination of the set (and its number of elements) of "bad" index $n$-sublattices of $\mathbb Z^2$.)
Computations: We set $r(N)=r_1(N)+r_2(N)+2r_3(N)+2r_4(N)+4r_5(N)$
where
$$r_1(N)=\sum_{d\vert N}\left(2\min(d,N/d)-1\right)$$
counts all triangular matrices in $\mathcal R(N)$,
$$r_2(N)=\sum_{d\vert N,d^2>N,d\equiv N/d\pmod 2}1$$
counts all matrices of the form $\left(\begin{array}{cc}a&b\\b&a\end{array}\right)$ with $0<b<a$ in $\mathcal R(N)$,
$$r_3=\sum_{1\leq \delta< \alpha,\alpha\delta<N,(\alpha+\delta)\vert N-\alpha\delta}1$$
counts all matrices of the form $\left(\begin{array}{cc}a&b\\b&d\end{array}
\right)$
with $0<b<d<a$ in $\mathcal R(N)$ (the factor $2$ before $r_3$ accounts
for the symmetry exchanging the role of $a$ and $d$),
$$r_4(N)=\sum_{1\leq \beta<\alpha,\alpha^2<N,2\alpha-\beta\vert N-\alpha^2}1$$
counts all matrices of the form $\left(\begin{array}{cc}a&b\\c&a\end{array}\right)$ with $0<c<b<a$ in $\mathcal R(N)$ (the factor $2$ before $r_4$ accounts
for the symmetry exchanging the role of $b$ and $c$) and
$$r_5(N)=\sum_{2\leq \delta<\alpha,\alpha\delta<N}\sum_{k,k\vert N-\alpha\delta,\alpha<k<\alpha+\delta}1$$
counts the generic type $\left(\begin{array}{cc}a&b\\c&d\end{array}\right)$
with $0<c<b<d<a$ in $\mathcal R(N)$ (the factor $4$ accounts
for the symmetric roles of elements in the pairs $(a,d)$ and $(b,c)$).
Formulae for $r_3,r_4,r_5$ are based on the identity
$$\det\left(\begin{array}{cc} a&b\\c&d\end{array}\right)=ad-bc
=(a-c)(d-c)+c(a+d-b-c)=\alpha\delta+c(\alpha+\delta-\beta)$$
with $\alpha=a-c,\beta=b-c,\delta=d-c$ entries of the positively reduced upper triangular matrix
$\left(\begin{array}{cc}\alpha&\beta\\0&\gamma\end{array}\right)$
obtained by subtracting the minimal entry $c$ from
all entries of a positively reduced matrix $\left(\begin{array}{cc} a&b\\c&d\end{array}\right)$
with $0\leq c\leq b<d\leq a$.
Final Remark: There are two possible generalizations to higher dimensions of positive reducedness. Both are problematic.
The matrix generalization: Positive reducedness can of course be defined for matrices of any size. Counting positively reduced matrices accordingly to non-zero determinants is however no longer interesting in dimension larger than $2$: For any $x$ in $\{0,1,2,\ldots\}$ the matrix $\left(\begin{array}{ccc}4+x&2+x&1+x\\x&1+x&3+x\\1+x&1+x&2+x\end{array}\right)$ has determinant $1$ and is positively reduced. (Counting positively reduced matrices accordingly to entry-sums is of course always possible and not very difficult in principle.) See comments by Dima Pasechnik for possibles remedies using total positivity.
The point of view of subgroups of finite index in $\mathbb Z^d$. One can of course define a "Newton polytope" and look at the intersection of its boundary with $\mathbb N(1,1,\ldots,1)$. This intersection can be integral (bad, I guess) or lie on a face of codimension larger than $1$ (acceptable or not?). Even in the case of a codimension $1$ face, there is a problem: It can happen that integral points of this face generate a sublattice of index $kd$ in $\mathbb Z^d$ for some integer $k>1$.
In any case, there is no longer an obvious way to associate a matrix to such a sublattice (there is of course an obvious sublattice associated to an integral matrix by considering $\mathbb Z$-linear combinations of rows.)
