The joint probability method (JPM) is the traditional way to determine the base flood elevation due to storm surge, and it usually requires simulation of storm surge response from tens of thousands of synthetic storms. The simulated storm surge is combined with probabilistic storm rates to create flood maps with various return periods. However, the map production requires enormous computational cost if state-of-the-art hydrodynamic models with high-resolution numerical grids are used; hence, optimal sampling (JPM-OS) with a small number of (~ 100-200) optimal (representative) storms is preferred. This paper presents a significantly improved JPM-OS, where a small number of optimal storms are objectively selected, and simulated storm surge responses of tens of thousands of storms are accurately interpolated from those for the optimal storms using a highly efficient kriging surrogate model. This study focuses on Southwest Florida and considers ~ 150 optimal storms that are selected based on simulations using either the low fidelity (with low resolution and simple physics) SLOSH model or the high fidelity (with high resolution and comprehensive physics) CH3D model. Surge responses to the optimal storms are simulated using both SLOSH and CH3D, and the flood elevations are calculated using JPM-OS with highly efficient kriging interpolations. For verification, the probabilistic inundation maps are compared to those obtained by the traditional JPM and variations of JPM-OS that employ different interpolation schemes, and computed probabilistic water levels are compared to those calculated by historical storm methods. The inundation maps obtained with the JPM-OS differ less than 10% from those obtained with JPM for 20,625 storms, with only 4% of the computational time.