High-resolution subbottom profiler (SBP) data are commonly recorded during academic or industrial surveys and analyzed as 2D profiles. However, there is growing interest in 3D imaging of the shallow seafloor subsurface, e.g., site surveys prior to the installation of offshore wind farms. Despite ongoing advancements in 3D systems, SBPs continue to be used mainly for 2D profiling. In other seismic imaging applications, “pseudo-3D” seismic cubes have been successfully generated from dense 2D seismic surveys. We developed a novel open-source Python-based workflow to interpolate densely spaced 2D SBP profiles acquired with a hull-mounted parametric SBP into pseudo-3D cubes. This workflow is applied to a study area on the Chatham Rise east of New Zealand’s South Island, comprising numerous seafloor pockmarks and subsurface paleo-pockmarks. Our workflow consists of two stages: (1) processing of individual 2D profiles to compensate for existing vertical offsets and (2) sparse data interpolation by applying the iterative projection onto convex sets algorithm. We correct the residual static effect caused by the sea state during acquisition, compensate for the varying tidal elevations, and apply vertical time shifts to individual profiles to correct misties at profile intersections. The interpolation of the binned sparse 3D cube is performed in the frequency domain using the fast Fourier transform and executed in parallel processes, which significantly increases the computational efficiency. We validate our workflow by testing different bin sizes, sparse transforms, thresholding approaches, and varying total iterations, which noticeably impact the resulting interpolation quality. Our developed workflow generates pseudo-3D subsurface images from densely spaced SBP profiles, with numerous potential applications for academic and industrial surveys.