1.) It's clear that any rotation or reflection of a lattice produces an isometric torus. So we can assume that our lattices have bases of the form $\{(t,0),(x,y)\}$ where $t>0$ and $y>0$. We can of course also change the basis so that $0 < x \le t$. I then claim that the classes of metrics on the torus are in bijection with the choices of $t,x,y$ satisfying these conditions. You can prove this by calculating the distances of points in the lattice closest to the origin. These correspond to the lengths of the shortest closed geodesics on your flat torus, which are clearly isometry invariants.
2.) Topologically, the only way a torus can decompose as a (nontrivial) product is as the product of circles. So a metric torus, if it were a product, would have to be given as $S^1_a \times S^1_b$ where $a,b>0$ denote the lengths of the circle factors. It follows that the two shortest closed geodesics on the torus are just the circle factors intersecting orthogonally. But for a lattice as you've given, the two shortest closed geodesics come from $u$ and $v$ and intersect at angle $\arctan(y/x) \not= \pi/2$.
The core of the issue is whether the definition of the cross-product is invariant. The definition that Wikipedia provides uses coordinates, and is hence manifestly chart dependent. This is mostly fine for ordinary uses of the cross product, where we are mainly sticking to $\mathbb{R}^3$, but the context here is a bit different. You have a vector space, in your case $T_xM$, and you want to define a cross-product on this vector space in a coordinate invariant manner. It is not at all clear that Wikipedia's definition achieves this.
The question then is really this: given a $3$-dimensional inner product space $V$. How do we define a cross product on $V$ that does not depend on a particular choice of isomorphism with $\mathbb{R}^3$? There is actually a rather natural way of defining an invariant cross product, which is to use the Hodge star.
Definition: Let $(V,g)$ be a $3$-dimensional oriented inner product space. Given $\mathbf{u},\mathbf{v}\in V$, we define their cross-product to be
$$\mathbf{u} \times \mathbf{v} \equiv \star(\mathbf{u}\wedge \mathbf{v}),$$
where $\star$ is the Hodge star. This is a manifestly coordinate independent definition, which for $V = \mathbb{R}^3$ reduces to the ordinary cross-product.
Let $(\mathbf{v}_1,\mathbf{v}_2,\mathbf{v}_3)$ be a basis for $V$. Then in components, we have
$$\mathbf{u} \wedge \mathbf{v} = \frac{1}{2}(u^iv^j-u^jv^i)\,\mathbf{v}_i\wedge \mathbf{v}_j,$$
where we use Einstein summation notation for repeated indices. Avoiding a rather lengthy calculation, you can show that the Hodge star is given in coordinates by
$$(\star (\mathbf{u} \wedge \mathbf{v}))_k = \sqrt{g}\,u^iv^j\epsilon_{ijk},$$
where $\sqrt{g}$ is the square root of the determinant of $g$ in whatever coordinate system you're in. These are components of a covector, so if we raise the index (using the inverse metric, not the metric itself, as you did in your OP) we get
$$(\mathbf{u} \times \mathbf{v})^\ell = \sqrt{g}\,u^iv^jg^{k\ell}\epsilon_{ijk},$$
which is very similar to the expression from Wikipedia, albeit with the inclusion of the determinant factor.
Ultimately, that determinant is the source of all your troubles. When you write down an expression in coordinate form, you need to make sure that it is tensorial so that they define genuinely coordinate invariant quantities. The easiest way to make sure your expressions are tensorial (besides remaining manifestly coordinate invariant) is to make sure the constituent objects you use are tensors themselves. This is fine for the (inverse) metric and your vectors, but the Levi-Civita symbol is not a tensor.
As the name suggests, the Levi-Civita symbol $\epsilon_{ijk}$ is a symbol whose components are defined to be completely alternating. The Levi-Civita symbol exists independently of any choice of coordinate system and is not designed to be a coordinate invariant object that transforms appropriately (for example, it doesn't really make sense to raise and lower indices of the Levi-Civita symbol). However, a closely related concept is the Levi-Civita tensor (the Levi-Civita symbol is distinct from the Levi-Civita tensor, although people are unfortunately quite sloppy about the distinction), which I will denote with a tilde as $\tilde{\epsilon}$. This is the object defined as
$$\tilde{\epsilon}_{ijk} = \sqrt{g}\,\epsilon_{ijk}.$$
This is a bonafide tensor (in fact, the Riemannian volume form, which underlies the fact that cross-products give volume) which transforms appropriately. You can see then that the cross-product is given by
$$(\mathbf{u}\times \mathbf{v})^\ell = u^iv^jg^{k\ell}\tilde{\epsilon}_{ijk},$$
and everything works out correctly if the $\epsilon$ in the Wikipedia definition was intended to be the Levi-Civita tensor and not symbol (again, people are unfortunately quite sloppy about this). Now you can see that all components of the formula above are genuine tensors, and so the expression is guaranteed to be well-defined in any coordinate system. Try your calculation again, but this time include the determinant factor (and use the inverse metric), and hopefully you'll find that everything works out.
Best Answer
This is not an actual answer, but this is too long for a comment.
If you have already shown that all sectional curvatures of your torus $(S^1\times S^1,g)$ vanish (which I haven't checked), then the metric $g$ is flat by definition. From the classification of flat tori, there exist a lattice $\Lambda \subset \Bbb R^2$ and an isometry $f \colon \Bbb R^2/\Lambda \to (S^1\times S^1,g)$. The metric on the left is the one obtained from the Euclidean metric on $\Bbb R^2$.
Here is how to construct the corresponding lattice. Let $\gamma_1$ (resp. $\gamma_2$) be the shortest (resp. second shortest) closed unit speed geodesic in $(S^1\times S^1,g)$. Let $p_0$ be the point where they intersect. Up to translating the time parameter, $\gamma_1$ and $\gamma_2$, we can assume that $\gamma_1(0) = \gamma_2(0) = p_0$. Now, up to replacing $\gamma_2(t)$ with $\gamma_2(-t)$, we can assume that $\langle \gamma'_1(0),\gamma'_2(0) \rangle \geqslant 0$. Let $u = \gamma'(0)$ and $v = \gamma'(0)$. Let $\Lambda = \Bbb Z u \oplus \Bbb Z v \subset T_{p_0}(S^1\times S^1) = \Bbb R^2$ (this identification is canonical). The exponential map $\exp_{p_0} \colon \Bbb R^2 \to (S^1\times S^1,g)$ is a local isometry, and descends as an isometry $f \colon \Bbb R^2 / \Lambda \to (S^1\times S^1,g)$.
This will be possible if and only if $|u|=|v| = 1$, and $\langle u,v\rangle = 0$, i.e. if $\{u,v\}$ is an orthonormal basis of $T_{p_0}(S^1\times S^1)$. This is because two flat tori $\Bbb R^2 / \Lambda_1$ and $\Bbb R^2 / \Lambda_2$ are isometric if and only if there exists a linear isometry of $\Bbb R^2$ that sends $\Lambda_1$ to $\Lambda_2$.
Yes, this is true for any diffeomorphism from $\Bbb R^2 / \Bbb Z^2$ to itself (so that the comment you added in edit does not make sense).
Finally, let me comment that I computed the area of your torus thanks to WolframAlpha, which seems to be close to 2.6. The usual torus has unit area. I might have made an error somewhere, but it seems to me that you torus isn't isometric to the standard flat torus.