A search for gauge-mediated supersymmetry (SUSY) in final states with photons and large missing transverse momentum is presented. The data sample of pp collisions at $\sqrt{s} = $ 13 TeV was collected with the CMS detector at the CERN LHC and corresponds to an integrated luminosity of 35.9 fb$^{-1}$. Data are compared with models in which the lightest neutralino has bino- or wino-like components, resulting in decays to photons and gravitinos, where the gravitinos escape detection. The event selection is optimized for both electroweak (EWK) and strong production SUSY scenarios. The observed data are consistent with standard model predictions, and limits are set in the context of a general gauge mediation model in which gaugino masses up to 980 GeV are excluded at 95% confidence level. Gaugino masses below 780 and 950 GeV are excluded in two simplified models with EWK production of mass-degenerate charginos and neutralinos. Stringent limits are set on simplified models based on gluino and squark pair production, excluding gluino (squark) masses up to 2100 (1750) GeV depending on the assumptions made for the decay modes and intermediate particle masses. This analysis sets the highest mass limits to date in the studied EWK models, and in the considered strong production models when the mass difference between the gauginos and the squarks or gluinos is small.