Must be an easier way. When I use the change of coordinates you specified:
$\displaystyle u=x^2+y^2$
$\displaystyle v=x^2-y^2$
and transform the region to the u-v plane, I get the indicated region in the second plot below (the diagonal line emanating from the origin is an artifact of the plotting code) :
The domain in the u-v plane now looks a little easier to integrate over, but seems to me you'd still have to split it up into two integrals not to mention you'd have to figure out what each line is and adjust the integration limits accordingly. Also I can't say I'm 100% confident I've mapped it correctly as I used some pretty hack Mathematica code to make the mapping. Here it is if you're interested. We could of course easily change the code to mappings we do know are correct to verify the code if that was necessary.
Code:
rplot = RegionPlot[x >= 0 && y >= -1 &&
4 <= x^2 + y^2 <= 16 &&
1 <= x^2 - y^2 <= 9 && y <= 2,
{x, 0, 5}, {y, -2, 2}, PlotPoints -> 75,
AxesLabel -> {Style["x", 20],
Style["y", 20]}, Frame -> None,
Axes -> True]
transform = rplot /. GraphicsComplex[pnts_,
data__] :> GraphicsComplex[
({#1[[1]]^2 + #1[[2]]^2, #1[[1]]^2 -
#1[[2]]^2} & ) /@ pnts, data];
newplot = Show[transform, PlotRange -> All,
AxesLabel -> {Style["u", 20],
Style["v", 20]}, Frame -> None,
Axes -> True]
GraphicsArray[{{rplot, newplot}}]