Using the data on coherent $J/psi$ photoproduction in Pb-Pb ultraperipheral collisions (UPCs) obtained in Runs 1 and 2 at the Large Hadron Collider (LHC), we determined with a good accuracy the nuclear suppression factor of $S_{Pb}(x)$ in a wide range of the momentum fraction $x$, $10^{-5} leq x leq 0.04$. In the small-$x$ region $x < 10^{-3}$, our $chi^2$ fit favors a flat form of $S_{Pb}(x) approx 0.6$ with approximately a 5% accuracy for $x=6 times 10^{-4} - 10^{-3} $ and a 25% error at $x=10^{-4}$. At the same time, uncertainties of the fit do not exclude a slow decrease of $S_{Pb}(x)$ in the small-$x$ limit. At large $x$, $S_{Pb}(x)$ is constrained to better than 10% precision up to $x=0.04$ and is also consistent with the value of $S_{Pb}(x)$ at $langle x rangle =0.042$, which we extract from the Fermilab data on the $A$ dependence of the cross section of coherent $J/psi$ photoproduction on fixed nuclear targets. The resulting uncertainties on $S_{Pb}(x)$ are small, which indicates the potential of the LHC data on coherent charmonium photoproduction in Pb-Pb UPCs to provide additional constraints on small-$x$ nPDFs. We explicitly demonstrate this using as an example the EPPS16 and nCTEQ16 nuclear parton distribution functions, whose uncertainties decrease severalfold after the Bayesian reweighting of the discussed UPC data.