We present a nonlinear model order reduction (MOR) method for parametrized partial differential equations (PDEs) in aerodynamics that exhibit parameter-dependent shocks. Linear MOR is ineffective for shocked problems due to the Kolmogorov N-width barrier, which necessitates a nonlinear approximation. Furthermore, the presence of strong shocks can introduce instability in both the full-order model (FOM) and the reduced-order model (ROM) due to numerical oscillations. To overcome these challenges, we present a robust and reliable nonlinear model order reduction method for shocked flows which is composed of the following: i) a high-order discontinuous Galerkin scheme with shock-capturing based on PDE-smoothed artificial viscosity fields; ii) a nonlinear ROM based on registration-based nonlinear reduced approximation spaces and hyperreduction via empirical quadrature procedure [1]; iii) a dual-weighted residual (DWR) output error estimate for quantities of interest for the FOM and ROM; and iv) spatio-parametric adaptive mesh refinement and greedy sampling informed by the DWR error estimate. A particular emphasis is placed on understanding and effectively treating the interaction between the modified PDE stabilized by the artificial viscosity and the DWR error estimate, which plays a crucial role in certifying the FOM and ROM solutions as well as driving efficient adaptation and sampling. We demonstrate the proposed method using various two- and three-dimensional transonic and hypersonic aerodyanmic problems governed by the compressible Euler, laminar Navier--Stokes, and Reynolds-averaged Navier--Stokes equations.