Two three-dimensional numerical integration schemes are presented for molecular integrands such as matrix elements of one-electron operators occurring in the Fock operator and expectation values of oneelectron operators describing molecular properties. The schemes are based on a judicious partitioning of space so that product-Gauss integration rules can be used in each region. Convergence with the number of integration points is such that very high accuracy (8-10 digits) may be obtained with a modest number of points. The use of point group symmetry to reduce the required number of points is discussed. Examples are given for overlap, nuclear potential, and electric field gradient integrals.