We use adjoint tomography to invert for three-dimensional structure of the North Island, New Zealand and the adjacent Hikurangi subduction zone. Due to a shallow depth to the plate interface below the North Island, this study area offers a rare opportunity for imaging material properties at an active subduction zone using land-based measurements. Starting from a ray tomography initial model, we perform iterative model updates using spectral element and adjoint simulations to fit waveforms with periods ranging from 4–30s. In total we perform 28 L-BFGS updates, improving data fit and introducing Vp and Vs changes of up to ±30%. Resolution analysis using point spread functions show that our measurements are most sensitive to heterogeneities in the upper 30km. The most striking velocity changes coincide with areas related to the active Hikurangi subduction zone. Lateral velocity structures in the upper 5km correlate well with New Zealand geology. The inversion recovers increased along-strike heterogeneity on the Hikurangi subduction margin with respect to the initial model. In Cook Strait we observe a low-velocity zone interpreted as deep sedimentary basins. In the central North Island, low-velocity anomalies are linked to surface geology, and we relate velocity structures at depth to crustal magmatic activity below the Taupo Volcanic Zone. Our velocity model provides more accurate synthetic seismograms, constrains complex velocity structures, and has implications for seismic hazard, slow slip modeling, and understanding of volcanic and tectonic structures related to the active Hikurangi subduction zone.