Borehole data are more commonly used than other geotechnical methods for determining the 3D geospatial information of a given site. However, the inadequate acquisition of continuous geotechnical information over an entire site limits the usefulness of this technique. In contrast, 2D continuous geophysical tomography over a large area reveals geophysical characteristics that can be transformed into geotechnical information. It is therefore possible to obtain more reliable stratal information by combining borehole datasets and geophysical measurements using geostatistical methods. In this study, a 3D geostatistical integration method incorporating site-specific geotechnical variability was proposed to construct a 3D geospatial database for developing geological unit boundaries for engineering projects. First, borehole datasets, reviewed by a geologist, were modified to fit a geostatistical trend of geological strata using a cross-validation-based outlier detection method. Then, site-specific geological layer criteria were determined based on proper statistical comparisons to convert geophysical values into boundary criteria for each geological layer. Finally, 3D geospatial information incorporating site-specific geomaterial characteristics was determined using indicator kriging. The proposed framework was established as functional modules in a geographic information system (GIS) platform and applied to three testbeds in Korea, incorporating a variety of site conditions such as spatial variability of coastal and mountainous areas, to validate its applicability depending on various geophysical characteristics. Application of the proposed framework indicates increases in estimation accuracy of 3D subsurfaces using geophysical and borehole information. However, the geostatistical integration of the two datasets has limited application in areas of geological complexity, such as areas with faults or anisotropic geological strata.