Accurate calculation of Wannier centers, position matrix, and composite operators using translationally equivariant and higher-order finite differences