A search is presented for the pair production of higgsinos χ˜ in gauge-mediated supersymmetry models, where the lightest neutralinos χ˜10 decay into a light gravitino G˜ in association with either a Higgs h or a Z boson. The search is performed with the ATLAS detector at the Large Hadron Collider using 139 fb⁻¹ of proton–proton collisions at a centre-of-mass energy of s=13TeV. It targets final states in which a Higgs boson decays into a photon pair, while the other Higgs or Z boson decays into a bb¯ pair, with missing transverse momentum associated with the two gravitinos. Search regions dependent on the amount of missing transverse momentum are defined by the requirements that the diphoton mass should be consistent with the mass of the Higgs boson, and the bb¯ mass with the mass of the Higgs or Z boson. The main backgrounds are estimated with data-driven methods using the sidebands of the diphoton mass distribution. No excesses beyond Standard Model expectations are observed and higgsinos with masses up to 320GeV are excluded, assuming a branching fraction of 100% for χ˜10→hG˜. This analysis excludes higgsinos with masses of 130GeV for branching fractions to hG˜ as low as 36%, thus providing complementarity to previous ATLAS searches in final states with multiple leptons or multiple b-jets, targeting different decays of the electroweak bosons.