For example, you can just use HCubature.jl with a change of variables to map the domain to a box. For example, see: 2D integration over non-rectangular domain using cubature - #7 by stevengj
(You could also do nested calls to a 1d quadrature routine like quadgk
, though this requires some care for the tolerance of the inner integrals.)