Accurate estimation of the time-averaged shear-wave velocity of the upper 30 m (VS30) is essential for seismic site characterization and hazard assessment. Direct in situ measurements are preferred for site-specific studies but are impractical at regional or national scales due to their high cost, logistical complexity, and time requirements, necessitating predictive modeling. This study develops a high-resolution VS30 map for Israel by integrating over 500 field measurements with optimized terrain- and geology-based proxies in a machine learning framework. Feature selection identified intermediate-resolution terrain variables (120–300 m) as the most effective predictors. Among the examined models, CatBoost (CB2) achieved the highest accuracy, outperforming traditional slope- and lithology-based proxies by better capturing complex, non-linear relationships between predictors and VS30. Residuals between measured and predicted VS30, interpolated via kriging, further improved local accuracy and quantified spatial uncertainty. The resulting VS30 map offers enhanced regional coverage with uncertainty estimates, and the methodology is readily transferable to other geologically complex regions.