Abstract A shape‐preserving multimoment scheme is designed to solve the advection equation on a hexagonal grid. Two types of local degrees of freedom are defined, including a volume‐integrated average and six point values for each cell. A quadratic reconstruction polynomial is constructed using a single‐cell stencil to implement a third‐order scheme. A flux‐correction algorithm is proposed for shape‐preserving simulations. By adjusting the mass leaving a cell through modifying the point values, the solution of volume‐integrated average can preserve the local extreme values in non‐divergent flows. A smoothness indicator based on the weighted essentially non‐oscillatory concept is adopted to avoid losing accuracy in smooth regions. The proposed scheme is checked by widely used benchmark tests, and numerical results demonstrate that it has third‐order accuracy on the planar and spherical hexagonal grids and is free of non‐physical oscillations around discontinuities. The proposed scheme has the practical potential to build accurate and scalable advection equation solvers on various spherical polygonal grids for atmospheric general circulation models.