In geosciences, kriging is leading spatial interpolation, and co-kriging is the most commonly used method for accomplishing spatial interpolation of a target variable by incorporating information from a secondary variable. Co-kriging relies on the assumption of spatial stationarity, which may not hold true in all geospatial contexts, leading to potential inaccuracies in interpolation. The effectiveness of co-kriging can be compromised in areas with sparse data, impacting the reliability of interpolated results. Moreover, it can be resource-intensive when used for interpolation with a substantial volume of data, especially in the case of 3D interpolation. In this paper, we introduce a new method for spatial interpolation that considers two variables using a generative deep neural network. This approach utilizes a convolutional neural network with an encoder–decoder architecture, featuring a single encoder and two decoders to handle the two variables. Additionally, we introduce a loss function that facilitates the control over the relationships between the two variables. Traditional Deep Learning methods require prior training and labeled data, whereas the proposed approach eliminates this requirement and simplifies the interpolation process. In order to assess the performance of our method, we use two real-world datasets. The first one is a 2D dataset of total soil organic carbon combined with the Normalized Difference Vegetation Index. The second one is a 3D dataset that combines concentrations of Hydrocarbon and Fluoride obtained from hyperspectral analysis of soil cores with very limited number of boreholes. The experimental results demonstrate that the proposed method outperforms ordinary kriging and co-kriging, showing a significant improvement when both variables are used. We also demonstrate how the inclusion of the auxiliary variable serves as a means to mitigate the overfitting of the model.