In the hydrodynamic model description of heavy ion collisions, the elliptic flow $v_2$ and triangular flow $v_3$ are sensitive to the quadrupole deformation $beta_2$ and octupole deformation $beta_3$ of the colliding nuclei. The relations between $v_n$ and $beta_n$ have recently been clarified and were found to follow a simple parametric form. The STAR Collaboration have just published precision $v_n$ data from isobaric $^{96}$Ru+$^{96}$Ru and $^{96}$Zr+$^{96}$Zr collisions, where they observe large differences in central collisions $v_{2,mathrm{Ru}}>v_{2,mathrm{Zr}}$ and $v_{3,mathrm{Ru}}<v_{3,mathrm{Zr}}$. Using a transport model simulation, we show that these orderings are a natural consequence of $beta_{2,mathrm{Ru}}ggbeta_{2,mathrm{Zr}}$ and $beta_{3,mathrm{Ru}}llbeta_{3,mathrm{Zr}}$. We are able to reproduce the centrality dependence of the $v_2$ ratio qualitatively and $v_3$ ratio quantitatively, and extract values of $beta_2$ and $beta_3$ that are consistent with those measured at low energy nuclear structure experiments. STAR data provide the first direct evidence of strong octupole correlations in the ground state of $^{96}$Zr in heavy ion collisions. Our analysis demonstrates that flow measurements in high-energy heavy ion collisions, especially using isobaric systems, are a new precision tool to study nuclear structure physics.