I have a discretized region $\mathscr{R}$ and two functions $x(u,v)$ and $y(u,v)$. I wish to find (numerically of course) the values of $(u,v)$ such that

$$
\begin{align}
\left(x(u,v),y(u,v)\right)\in \mathscr{R}
\end{align}
$$

To make it clearer let me take the functions $ \left(\frac{1}{3} \left(u^2+v^2\right),\frac{1}{3} \left(u^2-v^2\right)\right)$ and a geographical region, with

   scaled = TransformedRegion[RegionResize[DiscretizeRegion[GeoGridPosition[Entity["Country","Portugal"]["Polygon"], "Mercator"]], 1/2], TranslationTransform[{9.5, -43.4}]]; 
   Show[ParametricPlot[{(u^2 + v^2)/3, (u^2 - v^2)/3}, {u, 0, 2}, {v, 0, 2}], scaled]

so I have the following, and the idea is to get the set of $(u,v)$s that produce the blue part.

Now, my initial idea was to use NSolve but it didn’t work. Also, I haven’t seen any good implementation of inverse functions on two variables. Any ideas ?