Unfortunately the 1D reasoning will not directly carry into 2D and beyond. Depending on the domain for your integral, you will have considerable latitude in the node choice - but of course the quadrature weights can be generated from polynomials as usual. The closest you could achieve to a classic Gaussian-type rule is to coordinate-transform your domain to a product of 1D intervals and construct the rule by tensor product and coordinate transform back to get a rule for the domain of interest. This will generally not be an exact for polynomials, but should still retain good convergence rates provided the Jacobian of the coordinate transformation is well behaved.
Put briefly: in 1D you have a nice alignment of many goals for a quadrature rule in that they can be simultaneously achieved with a single rule.
- semi-analytic expression for the rule (e.g. zeros of orthogonal polynomials)
- Stability of rule (positive weights)
- Symmetry properties
- Fast convergence (e.g. exact on polynomials below a certain order)
It turns out that in 1D, all four of these conditions may be satisfied simultaneously with Gaussian quadrature rules. In 2D and beyond one must generally trade between the above features. For example often the fastest converging rules (those satisfying criterion (4)) can be incredibly expensive to compute, and often fail in (2) and (3). Rules that satisfy (1) (such as my example using coordinate transformations) often fail (2),(3) and (4).