In simulation of fluid injection in fractured geothermal reservoirs, the characteristics of the physical processes are severely affected by the local occurence of connected fractures. To resolve these structurally dominated processes, there is a need to develop discretization strategies that also limit computational effort. In this paper we present an upscaling methodology for geothermal heat transport with fractures represented explicitly in the computational grid. The heat transport is modeled by an advection-conduction equation for the temperature, and solved on a highly irregular coarse grid that preserves the fracture heterogeneity. The upscaling is based on different strategies for the advective term and the conductive term, respectively. The coarse scale advective term is constructed from sums of fine scale fluxes, whereas the coarse scale conductive term is constructed based on numerically computed basis functions. The method naturally incorporates a coupling between the matrix and the fractures via the discretization, so that explicit transfer terms that couple solution variables in the fractures and the matrix are avoided. Numerical results show that the upscaling methodology performs well, in particular for large upscaling ratios, and that it is applicable also to highly complex fracture networks.