I have an algebraic system of two equations in two unknowns that I would like to solve with NSolve within a specified region.

However, while Mathematica solves the system if I omit the domain specification, it fails to do so when I include it.

Here is an example:

eq1 = (-0.9781476007338057` + Cos[x] Cos[y])^2 + 
    Cos[y]^2 Sin[x]^2 + (-0.20791169081775931` + Sin[y])^2 == 
    0.04370479853238872`;

eq2 = -0.058944236842231254` (-0.9781476007338057` + Cos[x] Cos[y]) - 
    0.20040259242104053` Cos[y] Sin[x] - 
    0.008317236704697833` (-0.20791169081775931` + Sin[y]) == 0;

NSolve returns four solutions:

NSolve[{eq1, eq2}, {x, y}]

{{x -> -3.13939, y -> 3.14158}, {x -> -3.12983, 
  y -> 2.72301}, {x -> 0.0022017, y -> 0.0000114031}, {x -> 0.0117594,
   y -> 0.418582}}

I am only interested in the one that has positive $x$ and is sufficiently different from $(0,0)$. Hence I try to solve

NSolve[{eq1, eq2, x \[Element] Interval[{0.01, 0.5}], 
  y \[Element] Interval[{-0.5, 0.5}]}, {x, y}]

which Mathematica returns unevaluated.

Can anybody help?