Learning the Hamiltonian underlying a quantum many-body system in thermal equilibrium is a fundamental task in quantum learning theory and experimental sciences. To learn the Gibbs state of local Hamiltonians at any inverse temperature (\beta), the state-of-the-art provable algorithms fall short of the optimal sample and computational complexity, in sharp contrast with the locality and simplicity in the classical cases. In this work, we present a learning algorithm that learns each local term of a (n)-qubit (D)-dimensional Hamiltonian to an additive error (\epsilon) with sample complexity (\tilde{O}\left(\frac{e^{\mathrm{poly}(\beta)}}{\beta^2\epsilon^2}\right)log(n)). The protocol uses parallelizable local quantum measurements that act within bounded regions of the lattice and near-linear-time classical post-processing. Thus, our complexity is near optimal with respect to (n,\epsilon) and is polynomially tight with respect to (\beta). We also give a learning algorithm for Hamiltonians with bounded interaction degree with sample and time complexities of similar scaling on (n) but worse on (\beta, \epsilon). At the heart of our algorithm is the interplay between locality, the Kubo-Martin-Schwinger condition, and the operator Fourier transform at arbitrary temperatures.