Geological maps play a crucial role as fundamental sources of information in the Earth sciences, providing valuable insights into mineral and groundwater exploration, vulnerability to natural hazards, and various other applications. These maps are constructed using numerical or conceptual models, which extrapolate data from limited geological observations. Geostatistical techniques have traditionally been employed to generate reliable predictions that consider the spatial patterns inherent in the data. One such method is indicator cokriging, that estimates the probability of lithologic units at unsampled locations by analyzing its spatial correlation with continuous auxiliary variables, such as remote sensing images. However, these approaches are labor-intensive and have limitations compared to deep neural networks, particularly in terms of feature learning, which refers to the ability of models to learn relevant data representations for specific tasks. In this study, we propose an architecture based on U-shaped networks (U-Nets), a type of convolutional neural network designed for image segmentation, effectively exploiting the information contained in the secondary variables while producing spatially constrained predictions. Our architecture involves two U-Nets: the first focuses on learning the underlying patterns in the secondary data, and the second incorporates the ground-truth data along with the output of the first network. Given the sparsity of geological data, we must map the samples onto a regular grid. To achieve this, we interpolated field samples to align with the satellite images' grid, where each pixel represents the probability of a specific lithologic unit based on its frequency. The proposed methodology employs a multiscale loss function, combining categorical cross-entropy and a dilation operator. Additionally, to assess the aleatoric and epistemic uncertainty, we used Monte Carlo dropout as a Bayesian approximation in the model. We applied this methodology in eastern Quebec, Canada, specifically at the geological contact between the provinces of Churchill and Superior, identifying 20 distinct lithologic units. The secondary variables included multispectral, radar, and magnetic total resistivity data, interpolated to a pixel size of 300 m. To address spatial autocorrelation during validation, we divided the study area into spatial blocks. Through the implementation of the Spatially Constrained Bayesian Network (SCBN), we demonstrate its potential in generating accurate and field-data-constrained lithological maps, with an accuracy exceeding 0.7 in 11 out of 20 classes in the validation set. This study highlights the promising advancements of deep neural networks in geostatistics, particularly in handling complex spatial feature learning tasks, leading to improved predictive mapping outcomes.