Satellite-based evaluation of temporal change in cultivated land in Southern Punjab (Multan region) through dynamics of vegetation and land surface temperature

2.1 Study area and cropping calendar

The study area selected was Multan division located in Southern Punjab (Pakistan). This area lies between latitude 29° 27′ 21″ N to 30° 45′ 30″ N and longitude 71° 00′ 54″ E to 72° 58′ 43″ E approximately (Figure 1). The Multan division has four districts namely Multan, Khanewal, Lodhran, and Vehari districts. The study area is bounded by districts Jhang and Toba Tek Singh on the north, the districts Bahawalnagar and Bahawalpur on the south, the district Pakpattan on the east, and the district Muzaffargarh on the west. Multan division lies at Nili Bar, which is situated among Ravi, Beas, and Sutlej rivers. Many canals are present here that irrigate Multan region. The area is an alluvial plain and flat, which is an ideal area for agricultural activities. This irrigation network of canals makes the land very fertile. In the monsoon season, land area close to the Chenab river is usually flooded. All area with very hot summers and mild winters features has an arid climate. The average rainfall is roughly 186 mm (7.3 inch). The area witnesses some of the most extreme weathers in the country with highest recorded temperature approximately 52°C (126°F), during 2010 and the lowest recorded temperature approximately −1°C (30°F) [40]. Dust storms are a common phenomenon in the area. Multan region lies in cotton-wheat cropping zone of Punjab. Wheat is the major crop of Rabi season, and cotton is the major crop of Kharif season in Multan region.

Figure 1

Location of the study area (Multan division).

The cropping calendar is a tool, which provides complete data about crop phenology from sowing up to harvest; it contributes to devising the management plans required for sound agronomic practices and yield production assessments, etc. Cropping pattern of Multan division has two growing periods, namely kharif (summer) and rabi (winter). Rice, cotton, and kharif fodder (maize and sorghum) are grown in kharif season, whereas wheat and rabi fodder (mostly barseem) are grown in rabi season [51]. Sugarcane is an annual crop, which is mostly cultivated in February and September as well. Figure 2 represents the crop calendar of main crops grown in Multan division.

Figure 2

Cropping calendar adopted in the Southern Punjab. Rabi expresses the winter crop and kharif expresses the summer crop.

2.2 Image classification

Landsat images of Landsat 4, 5 Thematic Mapper (TM), Landsat 7 Enhanced TM Plus (ETM+), and Landsat 8 OLI were used for the analysis of cropping pattern and LST changes (Table 1). For this research, Landsat images were downloaded for 4 years (1990, 2000, 2010, and 2020) from official national aeronautics and space administration (NASA) Earth Explorer United States Geological Survey website (earthexplorer.usgs.gov).

Landsat images were preprocessed for georeferencing, layer stacking, and mosaicking by using GIS software: Arc GIS 10.4 and ERDAS imagine 15 [42]. Process of sub setting was used by using Arc GIS 10.4 software based on the study area [43,44]. Landsat images were then processed through supervised classification with digital topographic map of the Multan division. For our research, LULC classification was achieved with 0.97 convergence threshold and 50% iterations [45]. To assist Iterative Self Organizing Data Analysis Technique (ISODATA) algorithm, more refining of results was done by seeking farmers’ opinion considering various cropping patterns in Multan division. Using ISODATA clustering process at 30 m × 30 m spatial resolution, LULC classes were treated separately for both kharif and rabi cropping seasons. The comparative benefit is that this classification is appropriate for forestry, water resources, and agricultural studies. Furthermore, spectral signature of cultivated area for separate fodder did not exist from other source to confirm greatest accuracy of postclassification. Bare soil and built-up area were combined to one class (built-up/barren); therefore, five LULC classes were defined for rabi seasons including sugarcane, wheat, water bodies, built-up/barren, and rabbi fodder (barseem and vegetables), and six classes were defined for kharif season comprising sugarcane, cotton, rice, water bodies, built-up/barren, and kharif fodder (sorghum, maize, and vegetables, etc.) in the study area.

A survey of 420 farms was undertaken to collect information about crop information and LULC changes. A total of 35 union councils were selected from each district; from every union council, three villages were selected by simple random method. GPS Essential software was used to attain the latitude and longitude points. It is smartphone-based application, which permits collecting the georeferenced and digital field information. Areal distribution for all LULC classes was calculated for both cropping seasons from 1990 to 2020 over the whole study period. The relationship of temperature with NDVI for each cropping season was drawn. The details of methodology are given in Figure 3.

Figure 3

Flow diagram showing methodological and analytical steps.

2.3 Estimation of NDVI

Different Landsat images were used for the analysis of NDVI changes (Table 1). The NDVI was computed as follows [47,48]:

where NIR is near-infrared band (TM and ETM + band 4, OLI band 5), and RED is the red band (TM and ETM + band 3, OLI band 4). NDVI value ranges between −1 and +1.

Vegetation indices such as NDVI (time-based) profiles were used to combine various classes as well as to identify crop growth phases like sowing and harvesting [46]. The NDVI were derived for the mentioned study period (30 years) for Multan region. Such NDVI images were reclassified by using software Arc GIS 10.4 to compare the changes found in the study duration.

2.4 Classification accuracy assessment

Error matrix is a quite general widespread means to show accuracy outcomes. Different statistical processes of thematic accuracy can be drawn from the error matrix; for example, percentage of user’s accuracy, producer’s accuracy, and overall accuracy address the error produced by chance [49,50].

In equation (2), N sc is the number of sampling classes classified correctly, and N rc is the number of reference sampling classes. The KHAT ( k ˆ ) values are an estimate of how well RS classification is correct and agree with reference data conceptually, k ˆ can be defined as follows [51]:

In equation (3), A o is the observed accuracy, C a is the chance assessment, and C ag is the chance agreement. This statistic serves as an indicator of the extent to which the percentage correct values of an error matrix are due to “true” vs “chance” agreement.

2.5 Estimation of LST

LST is a general term relating to combined temperature of intact objects existing on the land. Numerous researchers used defined measurement on Landsat images to estimate the LST [52]. The LST was estimated from the OLI Landsat 8 thermal infrared band (band 10) with spatial resolutions of 30 m, respectively [53]. First, Lλ values were estimated to spectral radiance by equation (4) [54].

Lλ denotes spectral radiance and QCAL is quantized calibrated pixel value in digital number.

In the second step, the spectral radiance value was converted to temperature by using equation (5) [55].

where K ₁ and K ₂ values are 607.76 and 1260.56 for Landsat 4/5 (TM) and Landsat 7 (ETM+), and 772.88 and 1321.07 for Landsat 8 (OLI), respectively, and ln shows natural logarithm. In the last step, temperature in Kelvin (K) was converted to Celsius (°C) through equation (6) [56].

3.1 Areal distribution of major crops

Assessment of areal distribution for different LULC classes forms vital information for hydrological modeling mainly in watershed agriculture. The maps of LULC were obtained to identify the various crops in kharif and rabi seasons in Multan region. Figure 4 indicated the dominance of arable vegetation area in Multan division. Wheat and cotton were dominated crops during rabi and kharif periods, respectively. Rice was dominated in the downstream region along with sugarcane. For rabi seasons, a characteristic pattern in crop-cultivated land was perceived at different intervals in the study area. The cultivation of wheat showed most relative variation in area, that is, 1165106.73 ha (65.27%) and 1070775.22 ha (59.98%) in 1990 and 2020, respectively. Wheat area was decreased 5.28% during 1990–2020 in Multan region (Table 2). Minimum and maximum areas of sugarcane were 273827.24 ha (15.34% during 1990) and 221766.66 ha (12.42% during 2020), respectively. The cultivated area under sugarcane decreased 2.92% from 1990 to 2020. Minimum and maximum areas for rabi fodder (like barseem and vegetables etc.) were observed 10.4% during 1990 and 15.72% during 2020, respectively, whereas built-up/barren area were 7.03% during 1990 and 10.53% during 2020 (Figure 4). Build-up area increased 3.5% during 1990–2020 in Multan region. The water bodies showed relative variation in area, that is, 35054.01 ha (1.96%) and 23919.65 ha (1.34%) in 1990 and 2020, respectively, which decreased by 0.62% during 1990–2020. The decrease in water resources has also been one of the main reasons in the reduction of rice and sugarcane area in Multan region. It is worth observed that during the last 30 years, almost the (8%) wheat and sugarcane area transformed to fodder, cash crop maize, and other valuable vegetable crops as well as built-up area.

Figure 4

Crop identification maps for rabi and kharif seasons from 1990 to 2020.

During all kharif seasons, cotton was the main class in Multan division with a total cropped land ranging between 625973.33 ha (35.01% in 1990) and 508358.75 ha (28.48% in 2020). Cotton area decreased 6.59% during 1990–2020 in the study area. Rest of the area was under cultivation of sugarcane and rice with values ranging 15.54–12.65% for sugarcane and 21.57–19.96% for rice during the study period (Table 3). Minimum and maximum areas for kharif fodder (like maize, vegetables, etc.) were observed in 1990 (17.85%) and in 2020 (26.34%), respectively. The built-up area increased from 7.31 to 10.95% during 1990–2020. The shortfall in production is attributed to the decline in cotton area sown due to water shortage, cotton leaf curl virus, and high prices of synthetic spray of insecticides and pesticides. The cultivation and production of the maize crop increased due to the increase in area sown and as good economic returns and comparatively timely payment encouraged the growers to bring more area under maize cultivation. According to an earlier study [34], cropping pattern of different crops has changed and LST has increased during the last few years in the Multan district.

3.2 Accuracy assessment

Accuracy assessment in Table 4 shows the year wise detail of producer and user accuracies for various LULC classes including the kharif and rabi crops, water bodies, and barren land/buildings. The average of producer accuracy for kharif crops of cotton, rice, sugarcane, and fodder were 84.7, 84.2, 83.6, and 82.0%, respectively. Similarly, the average user accuracy for the four crop classes were 80.9, 80.9, 84.1, and 80.4% in the respective order. Similarly, the producer and user accuracies for rabi crops wheat, sugarcane, and fodder were 89.1 and 87.7, 80.3 and 80.0, and 78.6 and 81.3%, respectively.

In Table 5, we see year wise trend of overall accuracy and K values in addition to producer and user accuracies. The maximum overall accuracy of kharif season (84.8%) and rabi season (87.2%) were obtained in 1990, whereas the highest K value of kharif (85%) and rabi (87%) were obtained in 2000 and 2010, respectively. In Table 6, a brief overview of overall accuracy and K values of previously conducted studies (based on the use of different satellites and their sensor) is given. This is a sort of comparison and it is clear that Landsat had better values (58, 59, 60, 50, and 34) than MODIS. The current study had also reasonable overall accuracy (82.5–84.1%) and K value (80.0–81.5%). These values are compatible with already conducted studies and it shows the reliability of our results. The lower accuracy values in this comparison stem due to minor plot and combination of pixels sizes due to mixed cropping pattern during kharif season [57].

3.3 Cropping pattern and NDVI

Two depressions and two peaks in one cropping year were detected in Multan division. The two depressions seemed at end of month of April (on harvesting of wheat) and earlier in month of November (before sowing of wheat). The first peak corresponded to the extreme growth of wheat in February to March. The second peak corresponding to cotton at its extreme growth was observed during August to September. Sowing of wheat started at the end of November in the study area; however, an NDVI value became greatest during the month of February (Figure 5). The class “rice” in kharif period had a unique tendency. The early part was lengthy and bit slack because of growth of rice from June to July; whereas, the late part attained extreme height due to fast rice growth during the last week of August (Figure 5). Sugarcane was mostly cultivated during February. Its graph remained low and static in rabi season and reached to maximum values during kharif season due to increased vegetative growth. The class “building” showed the minimum NDVI values during the year due to less reflectance, which is according to earlier results by Cheema and Bastiaanssen [57]. Outcome of this satellite-based study showed that during 1990–2020 cropping seasons, various cropping patterns were executed by the local farmers in Multan division, namely “Rice–Wheat–Rice,” “Cotton–Wheat–Cotton,” “Wheat–Sugarcane–Wheat,” and “Wheat–Maize–Wheat.” This pattern of cropping was similar with the actual dominant system as shown during survey of the study area.

Figure 5

NDVI time-based tendencies for main crops from 1990 to 2020.

As regards NDVI, its values ranged in 1990 from −0.05 to +0.6; however, in 2000, NDVI showed the minimum value −0.09 and maximum +0.58 (as shown in Table 7). The NDVI values ranged from −0.12 to 0.55 in 2010 and from −0.15 to 0.54 in 2020 in the study area. This NDVI trend showed that average NDVI values decreased from 0.27 to 0.19 during 1990–2020 in the study area (Table 7). The NDVI is one of the commonly used indices to distinguish between healthy vegetation and vegetated regions from stressed and nonvegetated lands. The NDVI values identify the amount of chlorophyll content found in vegetation, wherever maximum NDVI value show dense and healthy vegetation and least NDVI value show bare soil and sparse vegetation. Therefore, as water absorbed electromagnetic radiations, water areas gave negative NDVI values. Similarly, Figure 6 is based on categorization of different land uses on the basis of NDVI in the month of March. It shows that vegetation was dense in most parts of the Multan region and the buildings occupied the less area in 1990 (as shown by green color having NDVI values >0.50). This situation transformed during subsequent years. We see that the vegetation area was occupied by more buildings in 2000, 2010, and 2020. In 2020, the highest number of buildings and/or barren land could be seen especially in Multan district (as shown by red color having NDVI <0). More towns have been constructed in the suburbs of cities where vegetation previously existed, or the vegetation area has been occupied by public and private property developers leaving the land barren at a vast area [58]. This situation could be alarming from agricultural point of view if no proper planning is executed to accommodate useless or less productive land for town development instead of productive area [60].

Figure 6

NDVI map of the study area during March for 1990–2020.

3.4 LST

The previous data of maximum and minimum temperatures for the 30 years (1990–2020) of Southern Punjab were attained from NASA website (http://power.larc.nasa.gov). Minimum temperature trend showed a slight decrease from 1990 to 1994, a trend reversed in 1996; considerable changes in the mean minimum temperature from 19 to 21°C were observed from 2014 onward. Maximum temperature trend during 1990–1995 showed a slight decrease, a similar trend as maximum temperature during 2014, an absolute increase in maximum temperature onward 2020 (Figure 7). The average temperature was observed 27.5°C (1990) and 28.6°C (2020), so an increase of 1.1°C in the average temperature was observed in Multan division during 1990–2020.

Figure 7

Average temperature trend of the study area from 1990 to 2020.

LST in Multan region for the month of March is shown in Figure 8. The results showed that LST values increased during 1990–2020, which was due to more built-up area producing maximum thermal signal. For the year 1990, LST ranges dominated in Multan region were from 18.5 to 36.7°C, and for the year 2000, LST ranged from 19.6 to 36.2°C in the study area. However, in 2010, LST values ranged from 19 to 37.3°C, and in year 2020, LST values ranged from 19.2 to 37.8°C (Table 7). In Multan region, high temperature was perceived in the built-up area because this area is more developed in terms of infrastructure than surrounding area. It is observed that the average LST values increased from 27.6 to 28.5°C during 1990–2020 due to increasing built-up area in Multan region. These facts are also supported by the comparison of Figure 6 with Figure 8, and we see that in southern part of Multan region, the infrastructural development and fallowing of cropping land for town purpose has not only led to the decrease in the value of NDVI but also resulted in more LST. It is most useful to have knowledge and information of LST to estimate and assess the urban surrounding climate. According to an earlier study [34], high temperature was observed in built-up area and bare soil because this area is more developed in terms of infrastructure than the surrounding area. Similarly, low temperature was recorded in vegetation area such as forest, wheat, rice, cotton, and sugarcane. This is more of true that in areas proximal to city buildings and roads in Multan region, the agricultural production is appreciably decreased due to increasing temperature and less rainfall; hence, the vegetation pattern has been changed and people have shifted from agronomic crops (wheat, cotton, etc.) to horticultural crops (vegetables, etc.) and fodder [61]. According to Amin et al. [62], uneven increase in precipitation and asymmetric increase in temperature extremes in the future would also increase the risk associated with management of climatic uncertainties. According to Amin et al. [63], southern Punjab has much agricultural production and has great contribution in total economic condition of the country; therefore, any temperature variations may affect the agriculture productivity.

Figure 8

LST maps of the study area during 1990–2020.

3.5 Relationship of NDVI and LST

To identify and assess the relationship of NDVI and LST in the study area, regression analysis was used between NDVI and LST values from the survey points for the month of March; the regression R ² analysis showed that LST and NDVI were negatively connected. In cropping seasons, nature of control (R ² = 0.82, 0.77, 0.74, and 0.71) was observed in 1990, 2000, 2010, and 2020, respectively (Figure 9). Most of the cluster of NDVI values was found to be from −0.05 to +0.6 in 1990; however, it was shifted to −0.15 to +0.54 in 2020 (Table 7). The NDVI values for Multan region, was higher in 1990 to +0.6 and lowest in 2020 to −0.15. The LST values for Multan region, was higher in 2020 to 37.8°C and lowest in 1990 to 18.5°C. The LST depends on soil moisture, thermal characteristics and particle size of a soil, whereas NDVI trend depends on physio-chemical features of leaf and plant texture. Policy orientation, market, and other factors affect the change of agricultural planting structure and also affect the change of LULC. The change of LULC leads to the great change of underlying surface, and then affects the change of surface temperature. These results show that the vegetation-covered areas can also be decreased because of the surface temperature. The negative relationship between NDVI and LST identifies that the greater biomass of vegetation cover has lower LST [34]. It also implies that the territories showing the most reduced NDVI values have a less vegetal spread due to the urban extension, whereas the higher NDVI values have a thick vegetal spread, and in this way, LST increase with the abatement in the vegetal density [64]. Some studies show that intact land-covers slowed down the increase in the surface temperature for mid and high latitudes in the region of North America and Eurasia over a period of 20 years. But this finding may be effective for a narrow range of temperature increases. So, if the temperature increases too quickly, or if the temperature increase crosses some critical value, then the vegetation may decline. Statistical estimates are consistent with this potential and show that there is an inverted U-shaped relation between summer temperature and summer NDVI for most of the land covers in the region of Eurasia [65,66].

Figure 9

Regression analysis of NDVI and LST from 1990 to 2020.