Application of Remote Sensing Techniques and Geographic Information Systems to Analyze Land Surface Temperature in Response to Land Use/Land Cover Change in Greater Cairo Region, Egypt ()
1. Introduction
Cities are dynamic due to unavoidable changes that can be assigned to many factors. One of the main factors behind these changes is urban growth and population expansion [1] . As the population of a given area increases, the interest for new settlements continues increasing at the expense of other land cover classes, for instance, vegetation and barren lands. Land surface impacts that occur during the process of urbanization include, but are not limited to, soil compaction, vegetation reduction and change from permeable to impervious surfaces as buildings, parking lots and roads are constructed. A seeming lack of planning of land use/land cover (LULC) has been a problem at a local and regional scale, making it a major issue in the study of worldwide ecological change [2] [3] [4] . Such changes have many implications for human society, environmental resilience, and water issues, such as the alteration of runoff, infiltration and groundwater discharge [5] . In addition, poor water quality occurs due to a lack of planning with comprehensive arrangements or any consideration regarding their effects on nature. Nevertheless, an increase in land surface temperatures (LST) is one of the key effects of LULC changes [6] - [11] . Increases in LST over the past several decades are considered a major issue in urban regions, due to the conversion of vegetation cover to impervious cover [12] , which in turn has a negative impact on people [13] , affects many environmental processes and modifies the degree of solar radiation’s absorption, evaporation rates, desertification, air pollution, albedo, heat storage, wind turbulence and many aspects of the water balance [14] . Therefore, the impact on environmental processes cannot be well-understood and mitigated without understanding the impacts of climate change, the interaction between the earth and the atmosphere and knowledge of land use/land cover change at various scales that drive them [2] .
Using remote sensing data in conjunction with Geographic Information Systems (GIS) proved effective for mapping urban areas, modeling urban growth and monitoring the dynamic changes of LULC [15] [16] . Remote sensing (RS) provides medium and high spatial, spectral and temporal resolution data with consistent and repetitive coverage of the earth’s surface [17] , and a high capability to extract change information from satellite data [12] . However, LULC change and LST can be monitored by traditional surveys and land based observation stations, as well as satellite data, because it is a time and cost-effective technique that can provide more information with respect to land use’s geographical distribution [7] . Satellite RS techniques, therefore, have become prevalent in monitoring change detection in both rural and urban regions [18] [19] [20] [21] . As a result, they have been widely used to evaluate LULC change with useful outputs and different scales [22] [23] .
Landsat imagery, in particular, is among the most widely used satellite systems. These datasets are available since 1972 from seven satellites in the Landsat series and they have been a major component of NASA’s Earth observation program with four primary sensors: Multi-spectral Scanner (MMS), Thematic Mapper (TM), Enhanced Thematic Mapper (ETM+) and Operational Land Manger (OLM). Landsat datasets have provided high resolution visible and infrared data, with thermal data and a panchromatic image. In addition, it supplies an extraordinary level of information on the classification of several earth components at large scale [24] [25] using a variety of automated change detection techniques and commonly applied classification algorithms (i.e. principle component analysis (PCA), unsupervised clustering, Hybrid, Fuzzy, Bayes and supervised classification). These change detection and classification techniques require personal experience and additional ancillary data with respect to study areas, i.e. a very high resolution aerial images and ground data, which can be used to construct a trustworthy dataset for different classification algorithms that can be used further in training samples and accuracy assessments [7] . Although ground data are considered to be the most reliable reference data, such data are often either not accessible or very costly. Therefore, a pre-defined statistical characterizations file for the image is created to store a per-pixel signature of a certain land cover class. This uses the stored information and the raw digital number (DN) of each individual pixel in the scene and converts them to radiance values. Several researchers and scholars applied similar techniques to achieve satisfactory results. For example, Landsat satellite images themselves were used to evaluate the performance of classification algorithms used to map forest clear cuts in the Pacific Northwest [26] . In addition, supervised classification maximum likelihood algorithms were applied to detect land cover change in a watershed in Pakistan and India with 95% and 92% overall accuracy, respectively [27] [28] . Although a high accuracy is obtained from these results, the presence of other reference data are essential to evaluate the overall accuracy and performance of the created geospatial maps [7] .
Accurate mapping of LST is becoming more significant in providing information about surface physical properties [14] , and the use of satellite images has become the predominant way to monitor LST on local and regional scales [9] [16] . The use of thermal remote sensing data were first demonstrated by Rao in 1972 for monitoring urban areas in the mid-Atlantic coast of the USA [29] . The contributions of RS and GIS have since been used to evaluate and model LST in many regions with several climatic conditions by various scholars using a diversity of thermal infrared (TIR) sensors. For instance, LST and Normalized Different Vegetation Index (NDVI) were evaluated to compare the spatial occurrence of droughts over the geo-botanical zone of Mongolia using the NOAA- Advanced Very High Resolution Radiometer (AVHRR) [30] . Other scholars cal- culated the mean target brightness temperature and cloud cover fraction (CCF) derived from the Geostationary Operational Environmental Satellite (GOES-8) to find the correlation between LST and NDVI over North America [31] . In addition, different algorithms were applied to retrieve LST from different satellite data, for instance, Landsat TM, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) and Moderate Resolution Imaging Spectroradiometer (MODIS), in different regional scales [2] [32] [33] . Also, the patterns of LULC change were identified in some studies, followed by an investigation for the impact of these changes on LST [6] [34] .
Cities located in semi-arid and arid regions require more attention to be better evaluated and understood [35] . In many developing countries of the world, including Egypt, there are limited regional figures on land expense for monitoring urban expansion. Studies revealed that settlements in developing countries grow five times as fast as those in developed ones [15] . The present study focuses mainly on change detection evaluation of LULC and LST in the Greater Cairo Region (GCR) of Egypt for the past several decades, from 1990 to 2016. The GCR contains the largest portion of facilities and services, enabling the foundation of dwelling places for qualified work forces that are generally found close to and within the city [36] . Moreover, due to the suitable topographic and geologic setting, the areas surrounding GCR showed the highest proportion of urban expansion. This has contributed to the high rate of population growth, city expansion and extravagant development.
The main thrust of this paper is to analyze LULC change through classifications and post-classification change detection techniques by utilizing multi- spectral Landsat TM and OLI data of GCR for 1990, 2003 and 2016 through the integration of remote sensing and GIS. In addition, the use of these Landsat data to estimate LST in GCR will be evaluated through different models and algorithms, as described in detail in the methodology section. The current study aims to: 1) quantitatively delineate different LULC classes and evaluate the pattern of LULC change from 1990 to 2016 in GCR; 2) provide tools to reliably investigate the variation of LST values in relation to LULC change through time; 3) further evaluate the effect of vegetation on LST as derived from different algorithms for satellite imagery through an examination of the NDVI-LST and (Nor- malized Different Built-up Index) NDBI-LST correlation based on statistical analysis methods and the texture of LULC changes, to determine the main causes of these environmental changes; and 4) examine the potential and the accuracy of RS and GIS utilization in monitoring the spatial distribution of LULC changes. The information gleaned from the validated change detection outputs can help in understanding the dynamics of LULC change in order to help policy- makers predict and plan for future developments in GCR, achieve long-term sustainability of soil and water resources, address impacts of climate change, and therefore characterize the evolution of new hot spots for urban construction lands and infrastructure development.
Two of the most important underlying premises of the objectives tested in the investigation is the opportunity of obtaining a LULC and LST maps from synoptic view of Landsat images over large spatial areas; and improve thermal studies in Egypt that will be used to justify subsides for law seeking to reduce impacts on thermal comfort in exciting urban areas.
2. Study Area
The selected study area for this research is the Greater Cairo Region (GCR), Egypt, which includes three sectors. The main sector is the metropolitan Cairo city on the eastern bank of the Nile River, parts of Giza City on the eastern bank of the Nile, and Qalyoubia, north of Cairo. The study area is located at 30˚00'N and 31˚20'E, in the middle and southern part, i.e. apex, of the Nile Delta Region, covering an area of 845,137 hectares (Figure 1). The Nile forms the admistrative division between Cairo and Giza sectors, running through the study area in a floodplain 9 to 35 km wide. This is constricted by hills on both the eastern and western sides, with desert areas extending in the eastern and western direction [37] , characterizing it as a subtropical climatic region with high temperatures and solar radiation, dry and rainless summers, and cold, moist and rainy winters [38] .
The GCR was selected because of its unique location and climatic conditions, with a diversity of historical heritages, making it one of the most dynamic urban regions in Egypt. It now represents about 23% of the total population of Egypt and 43% of the urban population. While the population expansion has grown tenfold in the whole country, the GCR has grown more than thirty fold in the last century and a half [39] . Half of this expansion has taken place on vegetated and rich agriculture land, while the other half occurs on the agglomeration fringes located at the borders of GCR. Little sporadic growth in the form of new communities has been created on what was desert land on the eastern district. Based on Central Agency for Public Mobilization and Statistics (CAPMAS) estimation, it appears that GCR has a population of about 20 million as of 2016 [40] (Figure 2).
3. Data Used
3.1. Satellite Data
Landsat 5 TM for 1990 and 2003 and Landsat 8 OLI 2016 were selected due to
(a) (b)
Figure 1. Greater Cairo Region (GCR), Egypt: (a) Location map, (b) The study area (Source: (a) ESRI online, (b) Landsat 8 Pan-Sharpened with DEM from Shuttle Radar Topography Mission (SRTM) [41] ).
Figure 2. Growth of GCR urban population during 1990-2016 according to Central Agency of Public Mobilization and Statistics (CAPMAS) [40] .
Table 1. Specification of Landsat satellite images. (*TIR = 120 × 30, data is acquired at 100 m and resampled to 30 m.)
their high spatial resolution for both multispectral and thermal bands, which benefits an accurate location of different land uses and monitoring LST. Due to their availability, three cloud-free Landsat images were selected to detect changes in the study area: August 4, 1990; August 8, 2003; and August 11, 2016 with scenes along the same path. The details of the Landsat data used in the current study are furnished in Table 1. All the satellite images were acquired during the summer season, intermediate to the agricultural growth season, in which most agricultural fields are green and active, which maximizes the spectral difference between these agricultural fields, urban areas and barren lands [7] . High-quality Landsat data acquisition is available from private and public sources.
3.2. Auxiliary Data
Reference data were compiled for each of the three years and then randomly divided for use in either classifier training or for accuracy assessment. These data were used to collect sufficient information for image preprocessing, evaluate the ground truth of a certain type of land use with its imagining characteristics [2] , and to determine the major types of land cover in the study area. The generated reference data include Egyptian topographic maps with a large scale (1:50,000) prepared by the Egyptian Military Survey, geologic maps, Digital Elevation Model (DEM) [41] and road networks [42] . Grouping of different spectral classes was done on the basis of land-cover types obtained from FAO-Land Cover Classification System (LCCS) of 2004, knowledge-based approaches and incorporated information from organizations and institutions of the Egyptian Governments. The ground truth data were in the form of reference data points used for assessing accuracy of the classification that were selected using high resolution GeoEye and QuickBird imagery [43] in addition to points collected during a field survey using Global Positioning System (GPS) receivers.
Image processing, such as image extraction, rectification, atmospheric correction for Landsat data, restoration and classification, and GIS analysis and interpretation were performed using a set of software to assure higher accuracy: Earth Resources Data Analysis Systems (ERDAS) Imagine 2014, the Environment for Visualizing Images (ENVI 5.3), the Integrated Land and Water Information System (ILWIS), ArcGIS 10.4 (ESRI) software, Python and Statistical Analysis System (SAS) software.
4. Methodology
The analysis includes image preprocessing, image classification, land cover indices (NDVI and NDBI) derivation and the evaluation of LST using thermal bands in the Landsat dataset. A flowchart of the research process is described and summarized in Figure 3.
4.1. Image Preprocessing
Both Landsat TM and OLI data are composed of independent single-band images. Therefore, it is necessary to combine these single-band images to a multi- band image of TM and OLI using a layer stacking tool. Landsat images were processed to a level-one terrain (L1T) corrected product, which provided radiometrically calibrated and orthorectified images using GCPs and DEM to attain absolute geodetic accuracy [34] . Therefore, the end result is a geometrically rectified image, free from any distortion related to the sensor, satellite and Earth’s surface [44] . The input Landsat data were georeferenced using the World Geodetic System 1984 datum (WGS-84) and the Universal Transverse Mercator (UTM) projection (within zone 36 North), as the study area lay in this region.
Although the data acquisition dates had clear atmospheric conditions for the study area, the three images were captured in different periods of time resulting in different atmospheric conditions. Hence, atmospheric corrections were conducted using the FLAASH module [45] , which was implemented in ENVI software. These atmospherically corrected images were clipped and subset to occupy the study area using an image subset tool. Other radiometric and spatial enhancement techniques, like histogram equalization, principle component analysis (PCA), edge enhancement and spatial filtering, were carried out on each image to improve the visual interpretability of the images.
Figure 3. Data processing flow chart depicting procedures applied for preparation of LULC maps and LST evaluation from Landsat datasets.
4.2. Optimum Index Factor (OIF) Calculation
This study used the analytical method of the Optimum Index Factor (OIF) to determine the best RGB band combination emerging from all bands of the Landsat images [46] , without the thermal band. The OIF is a statistical approach to rank all possible RGB color combinations of multispectral remote sensing data according to total variance within bands and inter-band correlation coefficients. Its role is to provide spectral information of the object, i.e. the highest OIF has the highest variance and lowest duplication for the scene, and therefore, contains the highest amount of spectral information about the scene. The algorithm used to compute the OIF was [47] :
(1)
where
is the standard deviation of band i, and
is the absolute value of correlation coefficient of any two arbitrary bands. For the Landsat 5 TM data (1990 and 2003), the top ranked RGB band combinations were band1/band5/ band7 (157) with OIF values of 60.784 and 56.431 for 1990 and 2003, respectively. The OIF calculation indicated that the band2/band3/band4 (234) RGB band combination had the highest spectral information with OIF value of 8155.43 for Landsat 8 OLI 2016 (Figure 4). Overall classification accuracy was high when these bands were utilized in the classification process instead of using all bands [2] .
4.3. Land Use/Land Cover Classification
A classification scheme had to be established before image classification. By computing average spectra of each class, a spectral characteristic of each land use class in each of the acquired data had been recognized, resulting in a classification schema comprised of 4 LULC level classes described in Table 2. As highlighted below, a number of classification approaches were evaluated for their effectiveness in large area classification [18] .
Figure 4. RGB band combination according to the highest OIF values of (a) TM 1990, (b) TM 2003, and (c) OLI 2016.
Table 2. Classification schema of LULC in the study area.
4.3.1. Unsupervised Classification
A combination of unsupervised classification methods were used to classify the study area. Images were first classified using the Unsupervised Interactive Self- organizing Data Analysis (ISODATA) algorithm to identify spectral cluster information from image data and convert image data to thematic data. This information contains average spectra for each of the identified LULC stored in a signature file, which in turn makes use of analyst with the help of the ground truth points and first-hand knowledge of the study area to recognize and assign spectrally uniform training data for a subsequent application of different supervised classification algorithms [7] . This clustering process was repeated several times through much iteration until a threshold was reached and there was no significant change in the cluster statistics or the maximum number of iterations was reached [48] . Clustering processes are highly automated with no direction from the users, so are ideal for large study area application.
4.3.2. Supervised Classification
The data were processed further using different supervised classification algorithms after they were classified using an unsupervised ISODATA algorithm. Training samples were first digitized from different representative classes to identify pixels of a single class. Grouping different spectral and spatial classes was done on the basis of LULC classes by utilizing reference data obtained from GCPs; auxiliary information and knowledge-based approaches collected from various resources, as mentioned before, were used to evaluate the statistical signature files of each LULC class [49] and ensure that there was minimal confusion of the land use to be mapped [50] . Different supervised classification algorithms were then carried out on the Landsat images; the executed algorithms included: Parallelepiped, Minimum distance, Mahalanobis distance, Maximum Likelihood and Support Vector Machine. Several algorithms were applied to identify the best for the study area location.
In order to increase classification accuracy and reduce classification error caused by confusion in spectral response of specific classes, the generalized images were spatially reclassified and refined for classification validation. Spectral confusion occurred due to the fact that several LULC have similar spectral response with respect to sensor characteristics especially in urban areas [20] . Therefore, data reclassification has to be applied to consolidate different LULC types using the image spatial and contextual properties. Reclassification was carried out based on auxiliary data and several GIS functions, for instance: digitizing, overlaying and region of interest (ROI) functionality to produce the last version of LULC maps for different years.
4.4. Post Classification Smoothing
4.4.1. Accuracy Assessment and Validation
Quantitative statements about accuracy assessment for the classification processes were an essential approach to validate how well the classification represented the real world and ensure the reliability of the information derived from LULC maps. Confusion matrices were computed to evaluate the relationship between the reference data used and the resulting classified LULC maps. Confusion matrices are one of the most popular ways to evaluate the overall classification accuracy providing information about producer’s accuracy or errors of omission (percentage of a specific LULC class on the ground which is correctly classified) and user’s accuracy or error of commission (percentage of a certain pixel class on the produced map corresponding to the actual class on the ground) [13] [51] . Percentage of the overall accuracy was computed using the following formula [12] :
(2)
Congalton (1991) was the first to point out that 250 reference pixels (±5%) are needed to construct the confusion matrix and to estimate the actual mean of accuracy assessment [52] . Therefore, 300 randomly selected reference pixels, placed on the classified images, for each time period were generated, representing a specific coordinate of the image. These points, distributed using the stratified random method, were then listed in two classes, one representing the class or reference values, while the other represented the actual LULC type. The percentage of the actual agreement of the automated classifier over a purely random assignment to classes was determined using a non-parametric Kappa coefficient [49] to remove the contribution of correct classification due to chance [18] . The Kappa coefficient for the different classification algorithms was evaluated by the following simplified equation [53] [54] :
(3)
where P(A) is the observed accuracy and P(E) is the chance agreement.
4.4.2. Land Use/Land Cover Change Detection
A multi-date post-classification comparison change detection method was employed to quantify the temporal change in LULC in the area of interest [55] . Three change detection statistics were obtained over time from the independent classified images for this research by conducting cross-tabulation analysis on a pixel-by-pixel basis, i.e. thematic overlay of the classified images [56] . The possibilities were (1990-2003), (2003-2016) and (1990-2016) to evaluate the matrix table of “from-to” change information that revealed the main gains and losses in each category of the study site.
4.5. Derivation of Land Surface Indices (NDVI and NDBI)
NDVI and NDBI were utilized to characterize the LULC classes and evaluate the relationship between these classes and LST. NDVI is the most commonly used index to express information about the density of vegetation, predict crop production, monitor drought, map desert encroachment [57] , and measure surface radiant temperature [2] . NDBI was first developed [58] to investigate the extent of imperviousness and built-up areas and map these areas, as it can highlight the urban distribution with a typically higher reflectance in the short-wave infrared region band than that of the near-infrared one [59] . NDVI and NDBI were computed using the following formulas from different wavelength regions of the Landsat data:
(4)
(5)
where NIR, Red and SWIR are the reflectance in the Near-Infrared band (0.76 - 0.9 µm), Red band (0.63 - 0.69 µm) and Short-wave Infrared band (1.55 - 1.75 µm), respectively, for Landsat 5 TM. However, for Landsat 8 OLI these differed slightly: Near-Infrared band (0.85 - 0.88 µm), Red band (0.64 - 0.67 µm) and Short-wave Infrared band (1.57 - 1.65 µm).
4.6. Land Surface Temperature Retrieval from Landsat 5 TM Data
Atmospheric correction was first required to eliminate the atmospheric effect from thermal bands, as the satellite imagery measures the radiance of surface features modified by the atmosphere [9] . Therefore, the Top of Atmospheric (TOA) radiance correction model was applied on Landsat 5 TM imageries for both 1990 and 2003. TOA radiance is a simple model based on the scene calibration data available from the imagery header file. Based on [60] , brightness temperature from Landsat 5 can be obtained first by the conversion of the digital number of band 6 to Top of Atmospheric (TOA) radiance using the following equation:
(6)
where
is TOA radiance,
is highest radiance corresponding to
(DN = 255),
is lowest radiance corresponding to
(DN = 0), and
is the quantized calibrated pixel value of band 6 in DNs.
The thermal band can then be converted from TOA radiance to effective at-sensor brightness temperature under the assumption that the Earth’s surface is a blackbody with a uniform emissivity and includes atmospheric effects using the following expression:
(7)
where
is at-satellite temperature in Kelvin, K1 is a calibration constant 1 (W/m2 sr μm), and K2 is a calibration constant 2 in Kelvin (Table 3).
Thereafter, the TOA radiance (
) was converted to reflectance measures, as
does not consider atmospheric effects. Assuming that urban areas behave as
Table 3. Landsat thermal band calibration constants.
a Planck surface, the expression to convert the TOA radiance to surface reflectance, correcting for solar irradiance, solar zenith and atmospheric effects is [9] :
(8)
The correct evaluation of LST was constrained to an accurate estimation of surface emissivity. In this work, we considered NDVI to calculate emissivity using the following formula [61] :
(9)
where a and b are obtained by regression analysis based on a large dataset [62] , a = 1.0094 and b = 0.047.
Finally, the LST corrected, in Celsius, for spectral emissivity is computed using the following expression [62] :
(10)
where
is the wavelength of emitted radiance (the average wavelengths = 11.45 µm) [63] ,
(1.438 × 10−2 m×K) with:
is Boltzman constant (1.38 × 10−23 J/K), h is Planck’s constant (6.626 × 10−34 J×s), and c is the velocity of light (2.998 × 108 m/s).
4.7. Land Surface Temperature Retrieval from Landsat 8 OLI
In the case of Landsat 8, TOA spectral radiance was computed using the radiance rescaling factors corresponding to each band provided in the metadata file using the following equation [44] :
(11)
where
is the radiance multiplier,
is the pixel value in DN and
is the radiance additive scaling factor for the bands obtained from the metadata.
A web-based atmospheric correction model was used to evaluate surface temperature by first converting the previously calculated TOA radiance to surface-leaving radiance, taking into account the atmospheric correction of thermal regions of Landsat 8 OLI [64] :
(12)
where
is atmospherically corrected radiance,
and
are upwelling and downwelling radiance, respectively (W/m2×sr×μm), and
and
are transmissivity and emissivity, respectively. These parameters can be assessed using the Atmospheric Correction Parameter Calculator online tool (https://atmcorr.gsfc.nasa.gov/). This uses the MODTRAN radiative transfer code that integrates algorithms to estimate atmospheric global profiles and parameters for a certain date, time, and location as the input [65] . Land surface emissivity was computed according to Equation (9). Even though the emissivity was calculated via NDVI in both Landsat 5 TM and Landsat 8 OLI, it was also preferable to use the same emissivity model for both Landsat datasets, hence avoiding uncertainty in the change in LST. Additional emissivity models introduced by [66] were also applied; however, results obtained corresponding to Equation (9) were considered the most reliable and the closest to the real life after validation, with only small differences found between the models.
Thermal Infrared bands of Landsat 8 OLI are converted from spectral radiance to effective at-sensor brightness temperature by converting the radiance using the inverse Landsat Plank’s law [60] :
(13)
where K1 is the band specific-thermal conversation constant 1 (W/m2 sr μm), and K2 is a calibration constant 2 in Kelvin (Table 3). Lastly, the emissivity- corrected LST, in Celsius, was retrieved using Equation (10) with the replacement of BT instead of TSensor.
5. Results and Discussion
5.1. Spatial Distribution and Accuracy Assessment of LULC
The Support Vector Machine (SVM) and maximum likelihood algorithms provided higher overall accuracy and kappa coefficients than other supervised classification algorithms (Table 4). Utilizing this observation, image processing and spectral characteristics, the final product combining the unsupervised and supervised classifications in which the spatial distribution and patterns of the LULC changes and persistence for the years 1990, 2003 and 2016, are shown in Figure 5. Spatial distribution patterns reveal that the area was dominated by deserts and barren lands, vegetation in the northern region and urban cover in the middle. Due to the heterogenic and dense vegetation cover in the north central part of the study region, we could not obtain higher overall accuracies than the ones presented, even after repeated classification with different algorithms.
The overall classification accuracies achieved for the images were found to be 90.3%, 96.5% and 94.9% with kappa coefficients of 0.85, 0.94 and 0.86 for 1990, 2003 and 2016, respectively. Note that in all classification algorithms, the vegetation class was responsible for producer errors; however, the urban class was the main reason for user errors (Table 5). On the other hand, classes of barren land and water were classified relatively accurately, approximately 98% or higher. The overall classifications in 2003 and 2016 are higher because of the availability of more detailed and higher resolution aerial reference images. Meanwhile, the use of OIF and enhancement techniques prior to classification increased the overall accuracy by 15% to 20%.
Table 4. Accuracy assessment of different supervised classification algorithms used for LULC maps in GCR.
Table 5. Accuracy assessment of the LULC classification results for GCR.
Figure 5. LULC map produced by classification processes for the years 1990, 2003 and 2016 showing the change for types of classes within the study area.
5.2. Land Use/Land Cover Change Detection
Change detection statistics were computed from each consecutive pair of LULC maps (1990-2003 and 2003-2016), and the results of these changes are furnished in Table 6 and Table 7. This shows the nature of change with respect to each class obtained from a matrix algorithm. Change detection analysis results show a sharp growth of 128% in the urban class during the 26-year period (1990-2016) (Figure 6). Significant differences appear to be related to barren land and vegetation land cover classes. Vegetation cover was reduced by 17,665 ha (14.3%)
Table 6. “From-to” LULC change detection statistics for 1990-2003 for GCR in hectare.
Table 7. “From-to” LULC change detection statistics for 2003-2016 for GCR in hectare.
Figure 6. LULC change in GCR during 1990-2016.
during 2003-2016 as compared to 14,432 ha (10.5%) during 1990-2003. Meanwhile, the barren land had a major decline of 30,669 ha (4.8%) and 24,822 ha (4.1%) during the two periods of 1990-2003 and 2003-2016, respectively, with a total amount of 55,491 ha during the entire period. These massive changes are related to desert-urbanization activities and construction of new housing developments, initiated by the Egyptian government in the early 1980s and that have since been accelerated [7] . This increasing trend in urbanization enhances the effect of human interference and reinforces that socio-economic forces are the main stimulus on these anthropogenic land changes, specifically around streams coming out from the Nile River. However, the reduction of vegetation cover and agricultural area, especially for urbanization purposes, illustrates the poor planning of farmland protection laws and the ignoring of environmental and ecological legislation implemented in the urban master plan. Water bodies, on the other hand, increased in area during 1990-2003, and then decreased again, due to the use of surrounding land for agriculture. Results obtained from this study were similar to that evaluated by [67] . They used three satellite images (1984, 2003 and 2014) to produce three LULC maps in GCR using the SVM algorithm. Results indicated that 13% of the vegetation was lost to urban areas between the period of 1984 to 2003, and 12% was lost between 2003 and 2014. However, only 3% of desert became urban areas in the first period, jumping to 5% between 2003 and 2014.
In order to better understand these “From-to” relationships, further GIS and statistical analyses were conducted. A post-classification comparison was conducted through cross-tabulation GIS modules to overlay the two LULC maps (1990 and 2016) to produce a LULC change detection map pointing out the spatial pattern of change for the 26-year timespan (Figure 7). Figure 8 and Table 8 show the percentage of different land covers in the GCR for the three time periods considered in this study. Results highlighted from these analyses showed two clearly recognizable trends; a) barren land and vegetation cover declined gradually and b) urban area increased drastically and rapidly (at the rate of 128% as mentioned before). The conversion patterns between different land cover classes to urban land cover are illustrated in Figure 9. This reveals that barren land was the main contributor in shaping urban area by an amount of 8.70% followed by vegetation land cover by a rate of 7.73% within the 26-year timespan (1990-2016). This emphasizes the importance of RS in conjunction with GIS in the study of LULC change detection providing essential information about the dynamic nature and patterns of spatial change of land cover.
Areas classified as urban in the northern part of the study area, particularly in Qalyoubia Governorate, are mainly disseminated as urban encroachments within areas classified as vegetation land cover (Figure 5). However, the change in the central part of the study area, in the political capital of the Cairo Governorate, indicates a rough spherical sprawling tendency with large regions containing concentration of many areas of localized change with a dense and granulated
Figure 7. Land cover conversion in GCR from 1990 to 2016.
Table 8. Change areas of LULC classes in the three time periods 1990, 2003 and 2016 in the study area.
Figure 8. Percentage of land cover types in GCR for the three time periods.
Figure 9. Contribution to the net change in the urban land cover in GCR (Area percentage %).
texture, especially during the period of 2003-2016. These are related to the establishment and implementation of new settlements, industries and communities at the expense of desert land that relies on surface water from the Nile River, for instance, El-Obour city and 10th of Ramadan cities to the east of Cairo, and 6 of October City in the western part. In general, these intense expansions occurred to accommodate the increasing population which caused the need for creating new jobs and maintaining food security, and is confirmed by Census data discussed in section 2 (Figure 2).
5.3. Land Surface Temperature (LST) Change and the Relationship with Land Use/Land Cover (LULC) Change
Figure 10 shows the spatial distribution and the pattern of change in LST throughout the different time periods of the study (1990, 2003, and 2016). The enormous increase of LST for all LULC types is highly evident, in addition to the wide range of LST values over the period from 1990 to 2016. Due to the dominance of desert and barren land in the study area, LST ranged from 28.78˚C to 47.11˚C, with a mean of 38.4˚C in 1990, 27.02˚C to 53.84˚C, with a mean of 40.6˚C in 2003, and 29.35˚C to 52.71˚C, with a mean of 42.1˚C in 2016. The estimated LST from different Landsat data were cross validated with near-surface surface air temperatures obtained from two meteorological stations on the same day of the obtained satellite data in the study area (Table 9). The higher estimates of LST from satellite data over the three time periods were due to the effect of surface roughness on the surface temperature [2] . Landsat data can be
Figure 10. Spatial distribution of GCR LST for the years (a) 1990, (b) 2003 and (c) 2016.
Table 9. Cross validation of the estimated LST from Landsat data with meteorological data for GCR.
Table 10. Average LST distribution (˚C) over different LULC classes in GCR for the years 1990, 2003 and 2016.
used to calibrate the distribution of LST in such dense places as GCR. Different algorithms for LST evaluation from Landsat were applied in this study to obtain accurate results. The calibration of LST should be refined with more data and in situ measurements of LST in future studies, as suggested by [32] .
GIS analysis coupled with image interpretation can help us to evaluate the relationship between the thermal signature and LULC types and highlight the impact of land cover changes on LST. We expressed the results of average LST in degree Celsius (˚C) by LULC classes over three time frames in Table 10. LST distribution of 2016 demonstrated that these new high LST were related to non- evaporating impervious surfaces, like the industrial and residential areas that had been converted from other land cover types [9] . Results revealed that the highest maximum and mean LST were associated with barren land (mean value of 39.81˚C to 43.69˚C) and high density urban areas (mean value of 36.83˚C to 41.47˚C), followed by vegetation (mean value of 33.05˚C to 36˚C) and, lastly, water cover (mean value of 30.59˚C to 32.89˚C), in all three periods. Due to the urban-warming effect and rapid urban growth in GCR, urban areas, such as industrial districts and commercial centers in the eastern and western side of Cairo, experienced an increase in LST by 4.91˚C over the entire period. This implies that these noticeable increases were due to high emissions of pollutants and multiple artificial heat sources, and the replacement of native vegetation areas that reduce the amount of heat stored through transpiration with other non- transpiring and non-evaporating surfaces such as concrete, metals and stones. These high density building surfaces experience high radiance temperatures, confirming the phenomena of urban-warming effect in which man-made materials in dense urban areas alter the superficial temperature and strongly link the urban category to higher LST in GCR.
It is also noted that vegetation cover showed the highest standard deviation values, reflecting the heterogeneous and complex nature of the vegetation cover with a wide range of surface radiant temperature. On the other hand, barren lands exhibited the lowest standard deviation due to the dry nature of these surfaces and a lack of wide variation in surface radiant temperature.
Table 11 shows how the newly formed lands reacted with regard to LST after transformation, excluding water cover that accounts for less than 1% of the study area. It can be clearly noted that the newly developed barren lands and urban areas have measured higher temperatures when transformed from vegetation cover with rates of 2.60˚C and 2.06˚C, respectively. On the other hand, vegetation cover tends to decrease the radiant LST in the conversion from either urban areas (1.90˚C) or barren lands (0.43˚C). Note that the transformation between barren lands and urban areas (and vice versa) had minimal effect on LST. In general, different land covers had different influences on the thermal distribution with a different magnitude, and LST acted as an important function of the change in LULC. Therefore, it is vital to increase the green area to strengthen study area protection.
Table 11. Average LST change in different LULC change types from 1990 to 2016 in GCR.
5.4. Analysis of Land Indices and Relationship with LST
Two indices, NDVI and NDBI, as mentioned in section 4.5, were derived to quantify the relationship between LST and land indices. The visual depiction of the spatial pattern of both NDVI and NDBI are shown in Figure 11 and Figure 12, respectively. It is observed that higher NDVI values correspond to dense vegetation areas in the central north of GCR, while lower values were observed in urban areas and barren land (Table 12). NDVI values are in the range of −0.525 to 0.79 in 1990, have a mean value of 0.04 and standard deviation of 0.15. These values dropped to −0.444 and 0.681 with a mean of 0.026 and a 0.14 standard deviation in 2003, and gradually decrease to be between −0.528 and 0.681, with a mean of −0.02 and a 0.08 standard deviation in 2016. As documented in the literature [31] , higher levels of NDVI were associated with lower values of LST. On the other hand, NDBI values were found to increase over the study period. For 1990, 2003 and 2016, the average NDBI was −0.043, −0.039 and 0.021, respectively. In general, decreasing surface transpiration through the reduction green
Figure 11. Spatial distribution of NDVI for the years (a) 1990, (b) 2003 and (c) 2016.
Figure 12. Spatial distribution of NDBI for the years (a) 1990, (b) 2003 and (c) 2016.
Table 12. Mean values of NDVI and NDBI for the years 1990, 2003 and 2016 for GCR.
Table 13. Correlation coefficient matrix from the indices and LST for the years 1990, 2003 and 2016 for GCR.
Table 14. Pearson’s correlation between LST and two indices at 0.05 significance level.
canopy cover and increasing impervious surfaces modified thermal behavior and were essential to the reduced value of NDVI and increased NDBI. This pattern can be clearly seen in Table 13 and Table 14, showing the statistical analysis and the Pearson’s correlation coefficient between the indices and LST at a 0.05 significance level. The results revealed that NDVI was negatively correlated with LST, indicating the impact of green cover on LST is negative, in which the more green areas, the weaker LST will be (Figure 13). In comparison, NDBI presents a high positive correlation coefficient with LST over the three time periods of the study. Therefore, urban areas can strengthen urban heat effects and increase LST [32] (Figure 14). It was interesting to note that the high negative correlation between NDVI and NDBI in the three years could be explained by the action of establishing urban settlements in favor of green cover.
In order to disclose the variance among LST and land surface indices, we took the 2016 image as an example to quantify the relationship among LST, NDVI and NDBI. Figure 15 shows the derived pixel values of LST and the two other indices based on West/East profile from the 2016 image, showing that lower LST are usually found in areas of lower NDVI; however the peak values of LST are consistent with higher values of NDBI along the profile. A multivariate linear regression analysis between LST and the indices was performed, leading to a
Figure 13. Correlation between NDVI and LST for the years (a) 1990, (b) 2003 and (c) 2016 for GCR.
Figure 14. Correlation between NDBI and LST for the years (a) 1990, (b) 2003 and (c) 2016 for GCR.
Figure 15. Correlation among LST, NDVI and NDBI from a West/East profile in 2016 imagery.
relationship among the variables as shown in Equation (14), with a correlation coefficient of R2 = 0.80, p < 0.001, and Root MSE = 1.82 at a 0.05 significance level.
(14)
The finding of Equation (14) showed a higher correlation coefficient in the multivariate linear regression than those of simple linear ones for the same year (R2 = 0.78 for LST-NDVI and R2 = 0.79 for LST-NDBI).
Pearson’s Correlation Coefficient between LST, NDVI and NDBI by Different LULC Types
In order to analyze the influence of specific LULC on LST accurately, Pearson’s correlation was evaluated between mean LST/NDVI (Table 15) and mean
Table 15. Pearson’s correlation between LST and NDVI by LULC type at 0.05 significance level.
Table 16. Pearson’s correlation between LST and NDBI by LULC type at 0.05 significance level.
LST/NDBI (Table 16) for LULC, pixel-by-pixel, in 1990 and 2016.
Results showed that the correlation between LST and NDVI were all negative in both 1990 and 2016 except for water coverage in 1990, likely due to lower amounts of pollutants in the 1990 water. However, the lowest negative correlation was found on barren lands in 1990 due to the high area coverage in the study region. On the other hand, the highest negative coefficient of the regression function was found to be in vegetation covers (−0.620) and dropped slightly for urban cover (−0.596). In 2016, barren area still showed the lowest negative correlation, however urban cover experienced the highest correlation coefficient (0.652), slightly higher than that of the vegetation one. Results of this analysis are consistence with other studies discussing the relationship between LST and NDVI [68] [69] . This indicates that by increasing NDVI, the LST of both vegetation and urban areas decreases more quickly than that of barren land cover.
However, the NDBI index showed a positive correlation with all LULC types except for barren lands in 1990 which had negative or almost no correlation with LST. Similar to NDVI, the highest coefficient was related to vegetation in 1990 and urban cover in 2016. In general, the correlation coefficient between LST and land surface indices for the whole study area showed a higher correlation than in the case of using indices as indicators for LST according to each LULC type. However, the moderate relationship obtained for each LULC by surface indices and LST can provide important information for preliminary studies; for example, they can be simply used as proxies for temperature and for better planning by policy makers in large areas.
6. Conclusions
In this study, multi-temporal Landsat satellite data were used to accurately monitor the spatial and temporal change of LULC and to study the impact of rapid urbanization on land surface temperature in the GCR in Egypt. Three Landsat dates, TM 1990, TM 2003 and OLI 2016, were acquired at the same time of the season (summer) due to the availability of reference data and to keep the weather factor as constant as possible. The study showed the effectiveness of the remote sensing techniques in conjunction with GIS to enable us to delineate the urban expansion due to the establishment of new settlements and to produce an accurate landscape change map in the study area. Different image enhancements, atmospheric correction, information extraction techniques, and unsuper- vised and supervised classification algorithms were performed on each Landsat image to ensure accurate image classification and LULC mapping. This revealed that SVM and maximum likelihood gave higher accuracies with rates of 90.3%, 96.5% and 94.9% for the years 1990, 2003 and 2016, respectively. The post-clas- sification comparison change detection method was employed to quantify the spatial change of land cover units. In addition, statistical “from-to” information was applied to quantify the magnitude of change through the entire 26-year time- span. Results demonstrated the most distinct change was related to vegetation cover that drastically decreased by an amount of 32,097 ha (23.3%) from 1990 to 2016. In the same time period, significant reduction in barren land by 55,491 ha (8.70%) occurred. On the other hand, urban areas, due to the construction of new industrial and commercial settlements, showed a considerable increase by 87,689 ha (128.3%), particularly in the central and northern parts of the study area around water resources. These two land covers, barren lands and vegetation, were the main contributors to form new urban areas.
LST evaluation through Landsat satellite images is easily recognizable. Different algorithms were applied to the thermal infrared data of different Landsat images to accurately calculate LST of different land covers. A very minor shift was found from different emissivity models, however, results obtained in this study were considered the most reliable based on cross validation. Results showed that mean LST values were higher in barren lands and urban areas than in the surroundings over the entire period. These anomalies were associated with settlements and industrial and commercial areas that experienced dense populations. Moreover, the most typical impact of rapid urbanization on LST was investigated within a single LULC. The change in LULC modified the radiant temperature of the surface. It was believed that the change in LST and climatic response was strongly related to the removal of vegetation cover and their replacement with non-evaporative surface. This can be concluded from the increase in the magnitude of LST by a rate of 2.06˚C in the areas that transformed from vegetation cover to urban and 2.60˚C due to the transformation of barren lands from green areas in the entire period of study from 1990 to 2016. In general, rapid urbanization was considered the major contributor to urban climatic warming which created the major spatiotemporal variation in LST, particularly due to the vegetation reduction and pollution expansion that could be attributed to settlement expansion resulting in a large amount of waste heat in turn affecting the surface energy budget.
Results from remote sensing studies show that LST and land surface indices, NDVI and NDBI together can identify the pattern of temporal variation and spatial distribution in urban thermal environments. The highest NDVI was found in vegetated areas while the highest NDBI was found in barren lands and urban areas. Statistical analysis showed a strong inverse relationship between LST and NDVI in contrast to a high positive one between LST and NDBI along different profiles in the study areas. These relationships dropped in the case of quantitative analysis among LST, LULC pattern and land surface indices. By way of conclusion, the study area reveals comparatively higher LST and NDBI, and lower NDVI over the period of study. These findings recommend the establishment of measures that can mitigate the strong effect of increasing LST on sustainable developments, population density control that is not limited to horizontal growth only, green coverage improvements like parks and gardens, and roof top area cultivation with horticultural plants that can alleviate the effect of LST. More multi-date images from the same season are also recommended to be investigated and evaluated, in a manner of providing more evidence of the thermal behavior on urban areas for better understanding of the impact of urbanization on LST. Moreover, RS satellite images are likely to be affected by cloud cover and other atmospheric effects, in addition to surface roughness, that in turn affect the DN values and therefore, it is highly recommended for future work that the integration of RS imageries from different sources with more land surface meteorological data be explored, and more attention on surface roughness be considered for more accurate results [9] .
In general, results proved the potential of multi-temporal Landsat images that can accurately quantify the change pattern in LULC and LST in GCR in Egypt. In addition, the integration of RS and GIS can provide a valuable opportunity for surveying, environmental monitoring and the nature of land cover change. Hence, the information gleaned from the change detection outputs can help in understanding the dynamics of LULC change in order to help policy-makers to predict and plan for future developments in GCR, achieve long-term sustainability of soil and water resources and its impacts on climate change, and therefore characterize the evolution of urban construction lands.
Acknowledgements
The authors would like to thank Dr. El-Sayed Ewis for his technical support and the Egyptian Government General Scholarship Programme administrated by the Egyptian Cultural and Education Bureau, Washington, DC, USA, for their financial support of this research.