One of the important parameters in urban climate is Land Surface Temperature (LST), which directly controls the Urban Heat (UH) effect  . Voogt  stipulated that the properties of a surface govern the surface energy balance and in turn the temperature. In urban areas, the paving and building materials used generally have a lower albedo than vegetated areas. Urban materials reflect less and absorb more sunlight, resulting in higher surface and air temperatures. The albedo of a city or a town depends on the surfaces’ arrangement, materials used for roofs, paving, coatings, and solar position  . The geometrical arrangement and the albedo of individual reflecting surfaces from building-air volumes influence the reflection of short-wave radiation  .
Data about the earth’s surface over a wide area can be extracted from satellite imagery at high temporal and spatial resolutions. Changes in land cover, specifically vegetation health, have been carried out using near infra-red as the vegetation portion reflects highly in this range  . In urban areas where land cover changes from vegetated to impervious surfaces, thermal temperature studies are important in determining thermal comfort. Land surface temperature (LST) measured from satellite imagery has been used in deriving climate models, urban heat studies, at a local and global scale. Land surface temperature is sensitive to vegetation and soil moisture  . The use of remote sensing in land cover or land surface temperature studies has increased as it offers a synoptic and broad view over an area instantaneously  . Many studies have been carried out using NDVI to determine the relationship between urban surfaces and thermal temperature from satellite imagery  . When using NDVI, higher values are an indication of high vegetation fraction. The correlation relationship between LST and NDVI is known to be negative due to the process of evapotranspiration  . According to research by Voogt and Oke  , surface thermal characteristics can be described using LST. They also noted that LST, which replaces specific air temperature, can be used as a surrogate for surface thermal characteristics.
The aim of the paper is to investigate what impact these changes in land cover have had on the land surface temperature in Upper-Hill using indices derived from satellite products. Another aim is to analyze which indices would be a better predictor for land surface temperature. The hottest months in Kenya are February and March while the coldest is in July until mid-August. The research was carried out during the months of January and February due to availability of cloud-free imagery.
2. Data and Methods
2.1. Study Area
Upper Hill, is approximately 4 Km from the city center of Nairobi and it has seen rapid developments over the years. Figure 1 shows the study area enclosed by the larger red polygon covering a total area of 15 Km2, which is a 1 Km buffer around Upper-Hill, which is enclosed by the inner polygon and has an area of 4 Km2. Upper-Hill, classified as Zone 1E, with a ground coverage ratio of between 35% - 60% and plot ratio of 150% - 300%, has transformed from residential to commercial/residential/office space. Businesses are choosing to relocate their businesses out of the city center due to traffic congestion, lack of adequate parking space, availability of amenities, just to name a few. Land prices have also increased with this increased demand. Zoning policies were also revised to accommodate the rising demand of development. The rise in development has resulted in a change of the urban landscape, with urban greenery replaced by hard landscape. Karanja & Matara  investigated the areas under concrete in Upper-Hill over a ten (10) year period between 2002-2012 and it showed a percentage increase of 3.84% with a subsequent decrease in green areas. Land cover changes have occurred, with ground coverage and plot ratios increased to accommodate the rising demand for housing.
2.2. Landsat Satellite Imagery
Landsat imagery acquired on February 13th 1987, February 10th 2002, January 5th 2015 and January 26th 2017 were used as the primary data source for deriving the land surface temperature and indices that would be used to determine the land cover changes. Data sets were acquired during the hot- season but due to cloud cover, imagery within the same month was not available. Data analysis involved using spatial metrics to determine the relationship with land surface temperatures derived from Landsat time-series imagery. Statistical modeling using ordinary linear regression to obtain this relationship was used. The study area is approximately 15 Km2. A buffer of one kilometer from Upper-Hill boundary was
Figure 1. Location of Upper-Hill area and buffer zone in Nairobi City.
defined and used to demonstrate the aforementioned relationship. Wu et al.  developed a 3-Dimentional Urban Index (3DUI) to compare urban densities with air temperatures over Taipei metropolis. Results showed that volumes of man-made constructions within 1000 m circular buffer from weather stations had strong direct impacts on ambient air temperature. This informed the choice of the one kilometer buffer to analyze the changing surface temperature of surrounding areas. The main methodological approach is presented.
2.3. Calculating Land Surface Temperature
2.3.1. Converting Landsat 5 TM to Landsat 7 ETM+
Level 1T Imagery acquired from Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 sensors were downloaded as L1T data from the USGS  website having a cloud cover of less than 10%. Atmospheric correction was carried out on the data. The projection used was WGS, UTM Zone 37 North. The imagery was first resampled where 0 values were assigned NoData values for all bands. DN values were converted to reflectance to obtain LST. The imagery available was for the months of January and February. Landsat 7 SCR off data was not used for gap filling in land surface temperature as using past imagery and resampling data would have resulted to inaccurate results on temperature. Landsat 5 TM data for the year 1987 was converted to an equivalent of Landsat 7 ETM+ data so as to calculate the Top of Atmosphere (TOA). Vogelmann et al.  described the process of converting from Landsat 7 ETM+ to Landsat 5 TM. To transform Landsat 5 TM to Landsat 7 ETM+ Equation (1) was used  .
where DN 7 is Landsat 7 ETM+ and DN 7 is Landsat 5 TM equivalent DN data. The slope and intercept are the inverse of those described by Vogelmann et al.  and they are band specific. The values were calculated by Firl & Carter  and are given in Table 1.
The values obtained were then used to calculate the Top-of-Atmosphere using Equation (2).
2.3.2. DN Values to Top of Atmosphere (TOA)
Spectral information is in digital number (DN), which has to be converted to reflectance values for analysis. Landsat 5 and 8 use the same formula for their bands which is different from Landsat 7. ArcGIS 10.4 was used to compute the Land Surface Temperature for both day and night time imagery.
Landsat 7 ETM+ consists of two thermal bands, 6a and 6b, obtained at 120 m and resampled to 30 m. Band 6a was employed since it has low radiometric variance and used in areas where vegetation cover is present. Equation (2) shows the conversion from DN to Top of Atmosphere (TOA) radiometric values for Landsat 7 and the results are shown in Table 2:
Table 1. Slope and intercept values  .
Table 2. Converting to TOA.
Lλ: spectral radiance,
Qcalmin: minimum quantized calibrated pixel value in DN,
Qcalmax: maximum quantized calibrated pixel value in DN,
Qcal: DN value of the pixel,
Lmin: minimum radiance detected by the sensor,
Lmax: maximum radiance detected by the sensor.
2.3.3. OLI and TIRS at Sensor Spectral Radiance
Landsat 8 consists of two thermal bands, band 10 and 11. USGS  recommends quantitative data analysis to be carried out on band 10 as band 11 is more contaminated by stray light. However there are other processing methods that one can use band 10 and 11 together such as split window. The 16-bit integer values in Landsat 8 can be converted to Top of Atmosphere (TOA) radiance using the equation developed by USGS  as shown in Equation (3). Table 3 shows the values from the satellite metadata used to calculate the TOA and at-satellite brightness.
Lλ is Top of Atmosphere (TOA) radiance in (Watts/m2*um),
ML is Band-specific multiplicative rescaling factor (RADIANCE_MULT_BAND_x where m is the band number),
Qcal is the digital number,
AL is the band specific additive rescaling factor (RADIANCE_ADD_BAND_x where x is the band number).
2.3.4. At-Satellite Brightness
To obtain the at-satellite brightness, Equation (4) was used for the analysis:
Table 3. Metadata for converting to TOA and at satellite brightness.
TB is the satellite brightness temperature in degrees Celsius,
K1 is the band specific thermal conversion constant (K1_CONSTANT_BAND_x, where x is band 10),
K2 is the band specific thermal conversion constant (K2_CONSTANT_BAND_x, where x is band 10).
Data outputs were exported as Geotiff imagery from ERDAS IMAGINE to ArcGIS 10.4 for further data analysis.
To determine the Land Surface Emissivity (LSE), NDVI was first computed using the reflectance values of red and near infra-red (NIR) bands of the Landsat image. The 16-bit integer values in Landsat 8 were converted to Top of Atmosphere reflectance as shown in Equation (5)  using values indicated in Table 4.
= Top of Atmosphere (TOA) planetary spectral reflectance without solar angle correction and is unitless,
MK = band-specific multiplicative rescaling factor (REFLECTANCE_MULT_BAND_x where x is the band number),
Qcal = the digital number of band_x,
AK = the band specific additive rescaling factor (REFLECTANCE_ADD_BAND_x where x is the band number).
does not have the solar elevation angle correction hence it is not the true TOA. Using the solar elevation angle from the metadata, conversion to the true TOA is done using Equation (6)  .
Table 4. Converting DN to reflectance values.
= top of atmosphere (TOA) planetary reflectance and is unitless,
= solar elevation angle obtained from the metadata.
NDVI was then calculated using reflectance values of the red and infra-red bands using Equation (7)
Equation (8) calculates the vegetation portion to obtain the LSE as shown in Equation (9).
= vegetation portion,
= normalized difference vegetation index,
= minimum NDVI,
= maximum NDVI,
where the minimum NDVI is the value for pure soil normally given as 0.2 and maximum NDVI is the value of pure vegetation given as 0.5.
LSE is then computed using:
2.3.6. Land Surface Temperature
Using the at-satellite brightness temperature and the Land Surface Emissivity, LST was computed in degrees Celsius as shown in Equation (10).
LST = land surface temperature,
TB = at-satellite brightness temperature,
λ = wavelength of emitted radiance (λ = 11.5 µm),
ρ = h*c/σ (1.438*10−2 m∙K),
σ = Bolzmann’s constant (1.38*10−23 J∙K−1),
h = Planck’s constant (6.26*10−34 J∙s),
c = velocity of light (2.998*10−8 m∙s−1).
2.4. Calculating Urban Indices
2.4.1. Normalized Density Vegetation Index
This was calculated as indicated in Equation (7). This was reclassified into two classes, vegetated and non-vegetated, where values of more than 0.2 were classified as vegetated. This was used as a mask to visually analyze land surface temperature within the study area as well as validate results obtained by calculating BDI as shown in Equation (12).
2.4.2. Built-Up Density Index
The Normalized Density Building Index (NDBI) density index is analyzed using the difference between reflectivity within the mid-infra-red (MIR) range and near infra-red (NIR) band range as in Equation (11). This is a dimensionless value where bright values indicate high density of built-up areas. NDBI was developed by Zha et al.  as a sensitive indicator for built-up areas.
Built-up density index (BDI) is a novel method that was developed by Lee et al.  whereby urban areas are extracted with accuracy of 92.6% and it uses the difference between NDVI and NDBI as shown in Equation (12).
3. Results and Discussions
LST, NDVI and BDI analysis was initially undertaken at a resolution of 30 meters. Values were averaged to obtain 90 × 90 meter pixel cells, which were then converted to point data in each cell. Averaging was done to obtain the average value of each of the indices amongst nine pixels to increase processing speed and easier interpretation of results statistically and visually over the study area.
3.1. Land Surface Temperature
Land surface temperature was analyzed from different months due to availability of data. In 1987, cloud cover in some parts of the study area affected the land surface temperatures recorded after the analysis. Figure 2 shows the land surface temperature from 1987, 2002, 2015 and 2017. In 1987, areas covered by clouds were determined to have recorded temperatures of between 11.81˚C and 22˚C. These areas were excluded in the linear regression. Land surface temperature increased from 1987 to 2017, with 2002 having the highest range of temperatures. 2015 and 2017 recorded lower temperatures than 2002, with January 2017 recording lower temperatures than February 2015. It varied between the period 1987, 2002, 2015, 2017 which may have been due to changes in climatic conditions in each of the years. In Figure 2, higher temperatures were noted in areas where built-up densities are higher (Figure 4) within the central business district (CBD).
Results from Figures 2(a)-(d) showed that land surface temperatures increased from February 1987 to February 2002, reaching a maximum of 38.86˚C within the CBD in 2002. LST was also observed to have increased in 2002, thereby decreasing in 2015 and 2017 respectively. Changes in temperature occurred in areas having concentrations of vegetation and impervious surfaces as shown in Figure 3. Temperatures in the eastern part of all images were generally higher than surrounding geographic areas due to the central business district and surrounding dense areas as shown in Figure 4.
Figure 2. Showing results for land surface temperature in degrees Celsius (a) February 1987; (b) February 2002; (c) January 2015; (d) January 2017.
An analysis was undertaken to determine areas where LST increased or decreased by subtracting two raster images as shown in Figure 3.
The three figures namely Figures 3(a)-(c) show land surface temperature differences calculated between the years 1987-2002, 2015-2002 and 2017-2015 respectively. The 1987 raster image was subtracted from the 2002 raster image. Areas that had negative values indicated a decrease in temperature i.e. temperatures in 1987 were higher than in 2002. The same was done where 2002 and 2015 images were subtracted from 2015 and 2017 imagery respectively. Areas that had positive values implied an increase in temperatures while areas with negative values implied a decrease in temperature. In 1987, areas that had cloud cover recorded temperature differences of between 10.69˚C - 19.99˚C, and these were excluded from the analysis. However, temperature differences between the highest and lowest recorded difference were highest between 2002 and 2015 with
Figure 3. Land surface temperature differences in degree Celsius between (a) 1987-2002; (b) 2002-2015; (c) 2015-2017.
17.79˚C, and the lowest occurring between 1987 and 2002 with 12.22˚C. Overall temperatures had increased during the years, with the greatest temperature increase in parts of Upper-Hill and surrounding areas towards the north and western areas. However in the Central Business District that is located east of the study area showed a decrease in temperature. These changes in land surface temperature could be as a result in an increase in impervious surfaces which absorb heat during the day. Night-time land surface temperature in this period would need to be analyzed to determine whether there is an expected increase in temperature in these urban areas as they radiate absorbed heat.
3.2. Normalized Difference Vegetation Index (NDVI) and Built-Up Density Index (BDI)
Normalized Difference Vegetation Index (NDVI) indicates the health of vegetation at any given time during observation while built-up density index enhances built-up areas more than normalized difference built-up index (NDBI). These two indices were used to examine temperature variations and also determine their relationship with LST.
NDVI values range between −1 to 1 with very low NDVI values of −0.1 to 0.1 indicating presence of rocks, sand or impervious surfaces. Higher NDVI values of 0.2 indicate presence of vegetation, with denser vegetation having values of 0.6 to 0.9, and deep water having a value of −1. In Figure 4 minimum and maximum NDVI values increased and decreased respectively from 1987 to 2017 as the spatial coverage and volume of vegetation decreased. By visually inspecting NDVI and BDI as well as altering the ranges, they were found to coincide and represent vegetation and built-up densities.
Figures 4(a)-(d) show results from analyzing NDVI values in the year 1987, 2002, 2015 and 2017. Results showed an increase in minimum values and a decrease
Figure 4. Results for NDVI (a) February 1987; (b) February 2002; (c) January 2015; (d) January 2017.
in maximum values from 1987 to 2002. This is an indicator that the amount of impervious surfaces had increased from 1987. By visually inspecting the results, it was observed that the central part of the study area where Upper-Hill is situated had transformed and with values being indicative of a decrease in vegetation due to urbanization.
BDI values range between −2 and +2 as NDVI and NDBI ranges are between −1 and +1. Since BDI is calculated by the difference between NDVI and NDBI, the values obtained that are indicative of clusters of built-up areas are dependent on the values obtained in each. In Figure 5, the minimum and maximum BDI values increased and decreased respectively from 1987 to 2017 due to changes in values calculated in NDBI and NDVI as shown in Figure 4.
Figures 5(a)-(d) show BDI results for the years 1987, 2002, 2015 and 2017 respectively. Concentrations of built-up densities are shown in the north-western part of the study area, and gradually increase within the central part of the study area though the years.
3.3. Minimum, Maximum and Mean Values
Table 5 shows the minimum, maximum and mean values for LST, NDVI and BDI. Values in areas with cloud cover in the 1987 Landsat image were not in
Figure 5. Results for BDI (a) February 1987; (b) February 2002; (c) January 2015; (d) January 2017.
Table 5. Minimum, maximum and mean LST, NDVI and BDI values.
cluded in the values indicated. The table shows that minimum values in LST, NDVI and BDI had increased which applied also to the maximum value, which had resulted from a change in land cover. LST minimum temperature had increased by 1.6˚C (6.7%) while maximum temperature increased by 3.65˚C (10.13%) between 1987 and 2017. This was affected by changes in land cover where vegetated areas were replaced by impervious surfaces. These areas are a great storage of heat hence absorbing heat during the day and releasing this heat at night.
3.4. Regression Analysis
LST, NDVI and BDI raster imagery were aggregated to 90 × 90 meter cells for statistical analysis in ArcGIS. In 1987 due to the presence of cloud cover, the data values were excluded from the analysis as these values were outliers and would have affected the model. Linear regression was carried out for all years, with LST as the dependent variable and NDVI, BDI as the independent variables. Figure 6 shows the results of the scatter plot and the regression line while R2 values are in Table 6.
The scatter plot showed a negative correlation between LST and NDVI and a positive correlation with BDI in each of the four years. A histogram of the residuals indicated they had a normal distribution. Table 6 gives the results for R2 which is an indicator of the goodness of fit of the model.
The R2 values for each of the years indicate that BDI is a better predictor of LST compared to NDVI. Adjusted R2 for NDVI was the lowest at 26.1% in 1987, increasing in 2002 to 49.3% which was the highest and reducing to 48.1% and 15.5% in 2015 and 2017 respectively. For BDI, 1987 had the lowest R2 value at 41.4%, increasing to 58.1% in 2002 and reducing to 53.1% and 28.9% in 2015 and 2017 respectively. Values from Table 6 hence indicate that built-up density is a better predictor of land surface temperature and has a stronger relationship than NDVI. The R2 value does not increase progressively from 1987 through to 2017 due to changes in climatic conditions and other factors such as land use changes. The variations in R2 value are commensurate with the variations in LST as analyzed in Figure 2.
Figure 6. Regression lines for LST with NDVI and BDI.
Table 6. R2 values of LST with NDVI and BDI.
Cloud cover and cloud shadows within the data set may have affected the output results as areas that were covered by cloud were omitted in the analysis. These areas could have had important information that would have changed the ranges in temperature, vegetation or built-up indices. Minimum and maximum land surface temperature has increased indicating a change in climate as the study was taken over a 30 year period. Linear regression indicates that BDI is a better predictor of LST than NDVI, however the model values also reduced with the decrease in temperature, indicating that there were other factors that were affecting land surface temperature. Results showed that understanding effects of changing land cover over a specific area enables mitigative measures to be determined, rather than analyzing a large area that has different geographical conditions. Future work will focus at comparing night-time and day-time land surface temperature to determine the aspect of surface urban heat islands.
 Feizizadeh, B. and Blaschke, T. (2013) Examining Urban Heat Island Relations to Land Use and Air Pollution: Multiple Endmember Spectral Mixture Analysis for Thermal Remote Sensing. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 6, 1749-1756.
 Bouyer, J., Musy, M., Huang, Y. and Athamena, K. (2009) Cities and Climate, Chapter 5: Mitigating Urban Heat Island Effect by Urban Design: Forms and Materials. World Bank, Washington DC, 164-181.
 Yue, W., Xu, J., Tan, W. and Xu, L. (2007) The Relationship between Land Surface Temperature and NDVI with Remote Sensing: Application to Shanghai Landsat 7 ETM+ Data. International Journal of Remote Sensing, 28, 3205-3226.
 Dash, P., Gottsche, F.-M., Olesen, F.-S. and Fischer, H. (2002) Land Surface Temperature and Emissivity Estimation from Passive Sensor Data: Theory and Practice; Current Trends. International Journal of Remote Sensing, 23, 2563-2594.
 Yuan, F. and Bauer, M.E. (2007) Comparison of Impervious Surface Area and Normalized Difference Vegetation Index as Indicators of Surface Urban Heat Island Effects in Landsat Imagery. Remote Sensing of Environment, 106, 375-386.
 Karanja, F.N. and Matara, S. (2013) The Transformation from Green to Concrete Cities; A Remote Sensing Perspective. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences (ISPRS), XL-1/W1, 163-166.
 Wu, C.-D., Lung, S.-C.C. and Jan, J.-F. (2013) Development of a 3-D Urbanization Index Using Digital Terrain Models for Surface Urban Heat Island Effects. ISPRS Journal of Photogrammetry and Remote Sensing, 81, 1-11.
 Vogelmann, J.E., Howard, S.M., Yang, L., Larson, C.R., Wylie, B.K. and Van Driel, J.N. (2001) Completion of the 1990’s National Land Cover Data Set for the Conterminous United States. Photogrammetric Engineering and Remote Sensing, 67, 650-662.
 Zha, Y., Gao, J. and Ni, S. (2003) Use of Normalized Difference Built-Up Index in Automatically Mapping Urban Areas from TM Imagery. Internal Journal of Remote Sensing, 24, 583-594.