Abstract: Regional characterization of the continental crust has classically been performed through either geologic mapping, geochemical sampling, or geophysical surveys. Rarely are these techniques fully integrated, due to limits of data coverage, quality, and/or incompatible datasets. We combine geologic observations, geochemical sampling, and geophysical surveys to create a coherent 3-D geologic model of a 50 x 50 km upper crustal region surrounding the SNOLAB underground physics laboratory in Canada, which includes the Southern Province, the Superior Province, the Sudbury Structure and the Grenville Front Tectonic Zone. Nine representative aggregate units of exposed lithologies are geologically characterized, geophysically constrained, and probed with 109 rock samples supported by compiled geochemical databases. A detailed study of the lognormal distributions of U and Th abundances and of their correlation permits a bivariate analysis for a robust treatment of the uncertainties. A downloadable 3D numerical model of U and Th distribution defines an average heat production of 1.5$^{+1.4}_{-0.7}$$\mu$W/m$^{3}$, and predicts a contribution of 7.7$^{+7.7}_{-3.0}$TNU (a Terrestrial Neutrino Unit is one geoneutrino event per 10$^{32}$ target protons per year) out of a crustal geoneutrino signal of 31.1$^{+8.0}_{-4.5}$TNU. The relatively high local crust geoneutrino signal together with its large variability strongly restrict the SNO+ capability of experimentally discriminating among BSE compositional models of the mantle. Future work to constrain the crustal heat production and the geoneutrino signal at SNO+ will be inefficient without more detailed geophysical characterization of the 3D structure of the heterogeneous Huronian Supergroup, which contributes the largest uncertainty to the calculation.