# Numerical Integral with Boolean as part of argument

I am trying to numerically calculate a partial surface area integral of an oblate spheroid in spherical coordinates. The area I want to calculate is determined by a boolean function in the argument of my integral. However, this integral fails to converge using certain input values. Here are the relevant bits of my code:

Req = 1.46055*10^9;
obl = 0.183;
a = 1.49598*10^11;
e = 0.0;
i = 45.*Pi/180;
ω = 0*Pi/180;
asc = 0*Pi/180;

R[θ_] := Req / Sqrt[Sin[θ]^2 + Cos[θ]^2 / (1 - obl)^2];
rx[f_] := (a (1 - e^2)) / (1 + e*Cos[f])*(Cos[asc]*Cos[ω + f] - Sin[asc]*Sin[ω + f]*Cos[i]);
ry[f_] := (a(1 - e^2))/(1 + e*Cos[f])*(Sin[asc]*Cos[ω + f] + Cos[asc]*Sin[ω + f]*Cos[i]);
rz[f_] := (a(1 - e^2))/(1 + e*Cos[f])*(Sin[ω + f]*Sin[i]);

test = Table[{f,
NIntegrate[
R[θ]^2*Sin[θ]*Boole[rx[f]*Cos[ϕ]*Sin[θ]+ry[f]*Sin[ϕ]*Sin[θ] + rz[f]*Cos[θ]>=R[θ]],
{θ, 0, Pi}, {ϕ, 0, 2 Pi},
Method -> "AdaptiveMonteCarlo",
MaxPoints -> 100000000,
PrecisionGoal -> 4,
MinRecursion -> 15,
MaxRecursion -> 30
]},
{f, 0, 2 Pi, Pi/20}
];


The NIntegrate function seems to be struggling to recognize the area set to zero by the Boole inequality. Should I be using a different method or different settings or something? I feel like I have tried them all with very little luck.

By the way, this integral is evaluated very precisely when the constant a is only about 3% of its current size. At its current value, however, the test list of integral values are highly imprecise:

The above table of integrals is supposed to result in a nice, smooth sine curve. How can I make this integral evaluate more precisely?