In modeling spatial processes, a second-order stationarity assumption is often made. However, for spatial data observed on a vast domain, the covariance function often varies over space, leading to a heterogeneous spatial dependence structure, therefore requiring nonstationary modeling. Spatial deformation is one of the main methods for modeling nonstationary processes, assuming the nonstationary process has a stationary counterpart in the deformed space. The estimation of the deformation function poses severe challenges. Here, we introduce a novel approach for nonstationary geostatistical modeling, using space deformation, when a single realization of the spatial process is observed. Our method is based on aligning regional variograms, where warping variability of the distance from each subregion explains the spatial nonstationarity. We propose to use multi-dimensional scaling to map the warped distances to spatial locations. We assess the performance of our new method using multiple simulation studies. Additionally, we illustrate our methodology on precipitation data to estimate the heterogeneous spatial dependence and to perform spatial predictions.