Geothermal resources have become an increasingly important source of renewable energy for electrical power generation worldwide. Combined Three Dimension (3D) Subsurface Temperature (SST) and Land Surface Temperature (LST) measurements are essential for accurate assessment of geothermal resources. In this study, subsurface and surface temperature distributions were combined using a dataset comprised of well logs and Thermal Infrared Remote sensing (TIR) images from Hokkaido island, northern Japan. Using 28,476 temperature data points from 433 boreholes sites and a method of Kriging with External Drift or trend (KED), SST distribution model from depths of 100 to 1500 m was produced. Regional LST was estimated from 13 scenes of Landsat 8 images. Resultant SST ranged from around 50 degrees C to 300 degrees C at a depth of 1500 m. Most of western and part of the eastern Hokkaido are characterized by high temperature gradients, while low temperatures were found in the central region. Higher temperatures in shallower crust imply the western region and part of the eastern region have high geothermal potential. Moreover, several LST zones considered to have high geothermal potential were identified upon clarification of the underground heat distribution according to 3D SST. LST in these zones showed the anomalies, 3 to 9 degrees C higher than the surrounding areas. These results demonstrate that our combination of TIR and 3D temperature modeling using well logging and geostatistics is an efficient and promising approach to geothermal resource exploration.