Numerous seismic activities occur in North Korea. However, it is difficult to perform seismic hazard assessment and obtain zonal data in the Korean Peninsula, including North Korea, when applying parametric or nonparametric methods. Remote sensing can be implemented for soil characterization or spatial zonation studies on irregular, surficial, and subsurface systems of inaccessible areas. Herein, a data-driven workflow for extracting the principal features using a digital terrain model (DTM) is proposed. In addition, geospatial grid information containing terrain features and the average shear wave velocity in the top 30 m of the subsurface (VS30) are employed using geostatistical interpolation methods; machine learning (ML)-based regression models were optimized and VS30-based seismic zonation in the test areas in North Korea were forecasted. The interrelationships between VS30 and terrain proxy (elevation, slope, and landform class) in the training area in South Korea were verified to define the input layer in regression models. The landform class represents a new proxy of VS30 and was subgrouped according to the correlation with grid-based VS30. The geospatial grid information was generated via the optimum geostatistical interpolation method (i.e., sequential Gaussian simulation (SGS)). The best-fitting model among four ML methods was determined by evaluating cost function-based prediction performance, performing uncertainty analysis for the empirical correlations of VS30, and studying spatial correspondence with the borehole-based VS30 map. Subsequently, the best-fitting regression models were designed by training the geospatial grid in South Korea. Then, DTM and its terrain features were constructed along with VS30 maps for three major cities (Pyongyang, Kaesong, and Nampo) in North Korea. A similar distribution of the VS30 grid obtained using SGS was shown in the multilayer perceptron-based VS30 map.