Ultracold atoms in optical lattices are one of the most promising experimental setups to simulate strongly correlated systems. However, efficient numerical algorithms able to benchmark experiments at low-temperatures in interesting 3d lattices are lacking. To this aim, here we introduce an efficient tensor network algorithm to accurately simulate thermal states of local Hamiltonians in any infinite lattice, and in any dimension. We apply the method to simulate thermal bosons in optical lattices. In particular, we study the physics of the (soft-core and hard-core) Bose-Hubbard model on the infinite pyrochlore and cubic lattices with unprecedented accuracy. Our technique is therefore an ideal tool to benchmark realistic and interesting optical-lattice experiments.