Here are a couple fully generic ways with region functions. I'm not sure which is best yet.

First, we might want to use scaled, or we might want to convert to a BoundaryMeshRegion via

Second, we might want to use InverseTransformedRegion[$reg,f,n$], which gives the region $\{p \in \mathbb{R}^n \mid f(p) \in reg\}$—in other words, exactly what you want, minus the bounds restrictions!

Or, with surface being your function of $u$ and $v$, we might want to use ImplicitRegion[surface[u,v] \[Element] reg, {{u, 0, 2}, {v, 0, 2}}].

Finally, if we discretize the result, we might want to discretize it via DsicretizeRegion or via BoundaryDiscretizeRegion. I'm not sure which is best.

Though, that part might depend on your application: if you don't need to discretize the resulting region, don't. Discretization of these regions takes a lot of time! The first one to return (DiscretizeRegion[uvImplicit]; see below) takes between 3 and 4 minutes on my laptop, depending on whether I'm connected to power. A lot of functions will work directly (and numerically) with the returned InverseTransformedRegion or ImplicitRegion anyway, so you might not need to discretize at all.

I wrote a test to time all different combinations. Here's that test, plus all the code you might want to use anyway.

ClearAll[scaled, bscaled, method, reg, surface, uvImplicit, uvInverseTransformed, inputregiontype]

scaled = 
  TransformedRegion[
   RegionResize[
    DiscretizeRegion[
     GeoGridPosition[Entity["Country", "Portugal"]["Polygon"], 
      "Mercator"]], 1/2], TranslationTransform[{9.5, -43.4}]];

bscaled = BoundaryDiscretizeRegion[scaled];

surface[u_, v_] := {(u^2 + v^2)/3, (u^2 - v^2)/3}

(* Example of both methods, here using scaled instead of bscaled *)

uvImplicit = ImplicitRegion[surface[u, v] \[Element] scaled, {{u, 0, 2}, {v, 0, 2}}]

uvInverseTransformed = InverseTransformedRegion[scaled, surface, 2]

(* Discretization test: *)

method["Implicit"] = ImplicitRegion[surface[u, v] \[Element] #1, {{u, 0, 2}, {v, 0, 2}}] &;
method["InverseTransformed"] = InverseTransformedRegion[#1, surface, 2] &;

inputregiontype["Mesh"] = scaled;
inputregiontype["BoundaryMesh"] = bscaled;

Table[
 Echo@AbsoluteTiming[{ (* Label: *) {disc, m, reg}, 
                       (* Output: *) disc[method[m][inputregiontype[reg]]]}
                     ],
 {disc, {DiscretizeRegion, BoundaryDiscretizeRegion}},
 {reg,  {"Mesh", "BoundaryMesh"}},
 {m,    {"Implicit", "InverseTransformed"}}
]

I'll keep you updated as these return...but in the meantime, maybe someone can suggest good option values for DiscretizeRegion that would speed things up, or a tweak to these approaches. I do think that an approach like this at least has the advantage of giving you "the right kind of thing" as a result—namely, a region. But I'm sure there are improvements that can be made.

So far, DiscretizeRegion[method["Implicit"][inputregiontype["BoundaryMesh"]]], i.e.

DiscretizeRegion[ImplicitRegion[surface[u, v] \[Element] bscaled, {{u, 0, 2}, {v, 0, 2}}]

is in the lead at just under 44 seconds.

Here's a fully generic way with ImplicitRegion.

First, we might want to simply use scaled, or we might instead want to convert to a BoundaryMeshRegion via

(As we'll see, this seems to improve performance.)

Second, with surface being your function of $u$ and $v$, we might want to use ImplicitRegion[surface[u,v] \[Element] reg, {{u, 0, 2}, {v, 0, 2}}].

(There is also InverseTransformedRegion[$reg,f,n$] instead of ImplicitRegion, which gives the region $\{p \in \mathbb{R}^n \mid f(p) \in reg\}$—in other words, exactly what you want, minus the bounds restrictions! However, as we'll see, there are computational problems with this option.)

Finally, if we discretize the result, we might want to discretize it via DiscretizeRegion or via BoundaryDiscretizeRegion.

Though, that part might depend on your application: if you don't need to discretize the resulting region, don't. Discretization of these regions takes some time! You might be able to use the region numerically even without discretizing it. some functions will work directly (and numerically) with the returned InverseTransformedRegion or ImplicitRegion anyway, so you might not need to discretize at all. However, if you're going to be doing a lot of different things with it, and if those computations are a bit slow, you might get better performance overall by discretizing from the outset.

I wrote a test to time all the different combinations:

• scaled vs. bscaled
• ImplicitRegion vs. InverseTransformedRegion
• DiscretizeRegion vs. BoundaryDiscretizeRegion.

The best was BoundaryDiscretizeRegion on an ImplicitRegion that referred to bscaled, which clocked in at around 44 seconds.

(* Best performance: *)

surface[u_, v_] := {(u^2 + v^2)/3, (u^2 - v^2)/3}

bscaled = BoundaryDiscretizeRegion[scaled];

BoundaryDiscretizeRegion[
  ImplicitRegion[surface[u, v] \[Element] bscaled,
                 {{u, 0, 2}, {v, 0, 2}}]
                        ]

Contrast with using scaled instead, which took almost 3 minutes. It performed about the same as using DiscretizeRegion on the outer level, but since BoundaryDiscretizeRegion seems to work better as an input, you might want to use BoundaryDiscretizeRegion anyway.

When I tried to discretize InverseTransformedRegion, it made my CPU temperature shoot up to 100ºC and wasn't any quicker. So if you're going to discretize it, I would recommend using ImplicitRegion instead of InverseTransformedRegion. However, there's always the chance that it will perform better with other operations; it all depends on how you'll be using this.

Maybe someone can suggest good option values for BoundaryDiscretizeRegion or DiscretizeRegion that would speed things up, or a tweak to this approach. I do think that an approach like this at least has the advantage of giving you "the right kind of thing" as a result—namely, a region. But I'm sure there are improvements that can be made.

For reference, here's the test I used. I wound up evaluating some individually when I wanted to skip some. I don't recommend running the test.

ClearAll[scaled, bscaled, method, reg, surface, uvImplicit, uvInverseTransformed, inputregiontype]

scaled = 
  TransformedRegion[
   RegionResize[
    DiscretizeRegion[
     GeoGridPosition[Entity["Country", "Portugal"]["Polygon"], 
      "Mercator"]], 1/2], TranslationTransform[{9.5, -43.4}]];

bscaled = BoundaryDiscretizeRegion[scaled];

surface[u_, v_] := {(u^2 + v^2)/3, (u^2 - v^2)/3}

(* Example of both methods, here using bscaled *)

uvImplicit = ImplicitRegion[surface[u, v] \[Element] bscaled, {{u, 0, 2}, {v, 0, 2}}]

uvInverseTransformed = InverseTransformedRegion[bscaled, surface, 2]

(* Discretization test: *)

method["Implicit"] = ImplicitRegion[surface[u, v] \[Element] #1, {{u, 0, 2}, {v, 0, 2}}] &;
method["InverseTransformed"] = InverseTransformedRegion[#1, surface, 2] &;

inputregiontype["Mesh"] = scaled;
inputregiontype["BoundaryMesh"] = bscaled;

Table[
 Echo@AbsoluteTiming[{ (* Label: *) {disc, m, reg}, 
                       (* Output: *) disc[method[m][inputregiontype[reg]]]}
                     ],
 {m,    {"Implicit", "InverseTransformed"}}
 {disc, {DiscretizeRegion, BoundaryDiscretizeRegion}}, 
 {reg,  {"Mesh", "BoundaryMesh"}} 
]

Here are a couple generic ways with region functions. I'm not sure which is best yet.

So far, DiscretizeRegion[method["Implicit"][inputregiontype["BoundaryMesh"]]], i.e. DiscretizeRegion[ImplicitRegion[surface[u, v] \[Element] bscaled, {{u, 0, 2}, {v, 0, 2}}] is in the lead at just under 44 seconds.

Here are a couple fully generic ways with region functions. I'm not sure which is best yet.

So far, DiscretizeRegion[method["Implicit"][inputregiontype["BoundaryMesh"]]], i.e.

DiscretizeRegion[ImplicitRegion[surface[u, v] \[Element] bscaled, {{u, 0, 2}, {v, 0, 2}}]

is in the lead at just under 44 seconds.