Identiﬁcation of Shift in Sowing and Harvesting Dates of Rice Crop ( L. Oryza sativa ) through Remote Sensing Techniques: A Case Study of Larkana District

: The present study aimed to determine the impact of climate variability on rice crops in terms of sowing and harvesting dates and crop period. The identiﬁcation of sowing and harvesting dates were spotted by mask identiﬁcation, variations in land surface temperature (LST) on a temporal scale in the respective months, and a ﬁeld-level social inquiry. The study was conducted during a time period (1994–2017), in which geo-referenced crop samples, farmer’s perception survey data, Landsat satellite images, and climate data of district Larkana were used. The analysis of satellite imageries revealed that on 20 June 1994, the rice was transplanted on 14.7% of the area of the region while it was only 7.1% of the area in 2017. Similarly, the area under rice crop in the ﬁrst week of July 1994 was 18.3% compared to 8.15% during the same period in 2017. However, in the ﬁrst week of October 2017, the rice crop was standing on 46.8% of the area while it was on 34.6% of the area during the year 1994 on the same date. This LST variation depicts a delay in the sowing and harvesting of the rice crop. This changing pattern is further conﬁrmed through mean LST. Mean LST ( ◦ C) has been increasing in the sowing period of rice crop from 31.9 ◦ C in June 1994 to 35.8 ◦ C in June 2017, and from 32.8 ◦ C in July 1994 to 36.8 ◦ C in July 2017. Furthermore, the LST decreased during the harvesting period of rice crop from 31 ◦ C in October 1994 to 28.6 ◦ C in October 2017. The present study quantiﬁes a delay of 15–30 days in sowing and harvesting dates of the rice crop in the district due to climate variability.


Introduction
Pakistan is located in a sub-tropical region and has arid and temperate climate. The country enjoys two seasons for agriculture, summer (April-September) and winter (October-December)-both of which are suitable for growing a variety of crops with optimal production. About 26% of the area of Pakistan is estimated to be cultivable [1]. Agriculture is the primary source of income for about 53% of the population of Pakistan [1]. Agriculture, specifically, rice production is one of the sectors, which plays a crucial role in Pakistan's economic development. During the 1950s, Pakistan's agricultural sector contributed about 53% of the gross domestic product (GDP). While between 2012-2014, it was reported that 21.4% dramatically dropped due to certain problems such as climate variability (changes in regional climate patterns, i.e., temperature and precipitation) impacts on agricultural production, deteriorating investment in the agricultural sector, urban migration, and a growing services sector in developing countries [2]. The fiscal year of 2012 to 2013 proved to be the third consecutive year in which the worst impacts of climate variability damaged crops [2]. Climate variability is severely affecting both agricultural and socioeconomic activities in many countries [3]. It has impact on crop yield in most developing countries, where climate is an essential element of agricultural production [4]. Past studies have also shown that climatic variability affected agricultural production and food security [5].
Climatic variability is one of the world's leading problems and Pakistan is also amongst the most vulnerable on the list of countries facing climate extremes [6]. Zhu [7] calculated that rising temperatures usually reduce yields by accelerating plant growth, causing early crop maturation and a decrease in high-yield production periods. Rice is more climate-sensitive than wheat, and other crops, especially those that are unable to adapt to these changes in climatic conditions, suffer the worst. Such climate variability has a direct influence on farmers incomes and thus on Pakistan's economy. Climate variability can influence the sowing and harvesting dates of crops as well as their growing period. In this study, we quantify the change in sowing and harvesting dates (crop calendar) of rice crops in the Larkana district region. The triangulation approach is used to investigate the change in the rice sowing and harvesting calendar from 1994 to 2017. The area of rice yields was estimated through the Landsat images and later this land use and land cover changes were verified through changes in land surface temperature (LST) from the area of study. Then, a social investigation was carried out using a farmer's perception survey to further validate the hypothesis.
LST was used to identify the changes in vegetation [8], urban heat island due to rapid urbanization [9,10], drought and soil moisture monitoring [11], and high-temperature risk assessments of crops at the regional scale [12]. Stoyanova et al. [11] studied the spatial-temporal variability of land surface dry anomalies in the climate aspect, and identified drought-prone areas through LST maps. The results derived from satellite LST data were used to identify the drought occurrence and its severity. They reported that the Soil Moisture Availability Index (SMAI) has a meaningful statistical relationship with LST [11]. The LST is a key factor controlling most physical, chemical, and biological processes on earth. Knowledge of the LST is necessary for many agricultural, environmental, and management activities related to surface resources. LST results from an energy exchange at the surface. Satellite-based monitoring is regarded as an essential prerequisite of regional and global observations of surface water, energy, and radiation budget. Different studies investigated that land use and land cover changes influenced the land surface energy balance significantly [13][14][15][16][17]. Irrigated agriculture (e.g., type of crop, planting pattern, water conditions, and crop canopy) has influenced this energy balance of the earth's surface environment [18] substantially. For example, irrigated agriculture impacted the range of diurnal air temperature [19]. Guoming et al. [20] studied the relationship between LST and paddy rice expansion. They found that LST decreases with the increased fraction of the rice paddy area more rapidly when the rice paddy is aggregated and accounted for by more than 80% of each study grid. On the other hand, land surface temperature was substituted for air temperature, and Kimura [21] examined the rice quality using LST satellite data from the Moderate Resolution Imaging Spectroradiometer (MODIS) during the ripening stage (August) in the Tottori prefecture.
LST datasets were applied in the lower Indus delta area to identify the land use and land cover changes. Different studies found that changes in the vegetative cover significantly impacted the LST in the Indus delta [22,23]. Rehman et al. [22] found that LST in the Keti Bandar area of the Indus delta is increasing due to the decrease in the vegetation of the area. They also reported that due to the mangrove's plantation drive of the Sindh government, LST was significantly decreasing in 2014 as compared to 2010. Siyal et al. [23] studied the Indus delta region and found a negative statistical correlation coefficient of determination R 2 0.65 between LST and the Normalized Difference Vegetation Index (NDVI).
The Crop Reporting Services (CRS) department collected the crop data annually through a randomized farmers survey. The discrepancies present in the field survey and data reporting creates a severe doubt to use such data for research purposes. Siyal et al. [24] studied Larkana district rice crop acreage and found that the Crop Reporting Sensing (CRS) overestimated (19-24%) the rice crop in Larkana district from 2006-2010. Thus, the annual crop reporting data is not reliable to use in finding the sowing and harvesting dates for rice crops. To prevail over these deficiencies in secondary data, we estimate the rice crop area using rice Crop Mask over Landsat imageries available in 1994, 2000, 2010, and 2017.

The Study Area
Larkana District was selected for this study because it is the leading rice-producing district in Pakistan's Sindh Province ( LST datasets were applied in the lower Indus delta area to identify the land use and land cover changes. Different studies found that changes in the vegetative cover significantly impacted the LST in the Indus delta [22,23]. Rehman et al. [22] found that LST in the Keti Bandar area of the Indus delta is increasing due to the decrease in the vegetation of the area. They also reported that due to the mangrove's plantation drive of the Sindh government, LST was significantly decreasing in 2014 as compared to 2010. Siyal et al. [23] studied the Indus delta region and found a negative statistical correlation coefficient of determination R 2 0.65 between LST and the Normalized Difference Vegetation Index (NDVI).
The Crop Reporting Services (CRS) department collected the crop data annually through a randomized farmers survey. The discrepancies present in the field survey and data reporting creates a severe doubt to use such data for research purposes. Siyal et al. [24] studied Larkana district rice crop acreage and found that the Crop Reporting Sensing (CRS) overestimated (19%-24%) the rice crop in Larkana district from 2006-2010. Thus, the annual crop reporting data is not reliable to use in finding the sowing and harvesting dates for rice crops. To prevail over these deficiencies in secondary data, we estimate the rice crop area using rice Crop Mask over Landsat imageries available in 1994, 2000, 2010, and 2017.

The Study Area
Larkana District was selected for this study because it is the leading rice-producing district in Pakistan's Sindh Province (  During the summer season (April to September), the average maximum temperature of the district is 42 • C, and the average minimum temperature is 31 • C, while during the winter season (October to March) the average maximum temperature is 21 • C and the average minimum temperature is 11 • C.  [26,27]. The reference crop evapotranspiration for the study is 1026 mm annually, and the potential crop water requirement for rice crops is 1303 mm. The cropping intensity reported was 88% annually and 50% during the summer season. The canal water allocation is 1860 L per second for a 1000 ha area [28].

Rice Crop Acreage Estimation
Landsat 5, 7, and 8 satellite imageries of path 152 and row 41 for June, July, September, and October corresponding to the years 1994, 2000, 2010, and 2017, were downloaded from the USGS GloVis (Global Visualization Viewer) website (http://glovis.usgs.gov/). The Landsat satellite imagery for 22 July 2000 was available, but the imagery has cloud coverage of more than 10%. Since cloud coverage interferes with image classification for estimating the vegetation area, the nearest Landsat imagery from 15 August 2000, was used to represent this year (see Table 1). Table 1. Landsat images used in the study and their radiometric coefficients, DOY (day of the year), d (earth-sun distance in astronomical units) and θ s (sun elevation) taken from the metadata files of the selected images.

S. No
Landsat Based on Landsat satellite imagery, rice crop masks for June, July, September, and October corresponding to the years 1994, 2000, 2010, and 2017 prepared to determine spatial and temporal variations in the rice crop area in Larkana District. The images were classified into rice and non-rice crop areas using visual interpretation of the multi-temporal Landsat imagery in combination with high-resolution imagery from Google Earth, and employing local expert knowledge, with the maximum likelihood classification algorithm in Arc GIS 10.3.1. The Landsat data must be preprocessed before analyzing for any change detection between rice and non-rice areas [29]. Therefore, the Landsat data were preprocessed in ERDAS Imagine 2016 and ENVI 5.1 to gap-fill the scan lines of missing data in the Landsat 7 imageries. An extract the area of interest (AOI), and the GIS extract-by-mask tool were used. After preprocessing the data, the Landsat imageries were classified using the ArcGIS 10.3.1. The supervised classification method was used to train the satellite data in this study due to the non-availability of field data. Initially, the Landsat data was classified into six different classes: water, rice, mango, guava, barren land, and vegetation. These divisions were determined using supervised classification with a maximum likelihood algorithm using field data samples; for every image, signature files were created. The six classes were then aggregated into two categories: rice and non-rice. The classified images were converted into a vector format from a raster to simplify area calculations. From the attribute table of the vector datasets, the rice crop area was determined based on the area covered by each class. The procedure was repeated for all selected Landsat imageries. After converting the classified Landsat images into vector format, the area under each class was calculated using a query builder tool. The accuracy level of image classification methods is reported in Table 2. The accuracy of the image classification depends on the classification technique and the number of training samples of each class. The classification method based on support vector machine accuracy is 6% better than the maximum likelihood method, but in this study maximum likelihood algorithm (MLA)-supervised classification method was adopted. The results obtained through this method may vary up to 6% and can be improved by adopting a better classification approach [30].

Land Surface Temperature
The mean LST of Larkana District was calculated from thermal bands of Landsat satellite imagery for the years 1994, 2000, 2010, and 2017. The mean LST was calculated by averaging minimum and maximum LST of the study area [23]. Landsat 5, 7, and 8 images were used in this study due to the long temporal time series. The significant difference between these products is the resolution of the products. The LST is retrieved based on the number of algorithm and emissivity model. There are different algorithms such as the Mono Window Algorithm (MWA) [31], the Radiative Transfer Equation (RTE) [32] method, the Single Channel Algorithm (SCA), and the Split Window Algorithm (SWA) [33,34]. The surface emissivity models used were categorized as: classification-based emissivity method, NDVI-based emissivity method, multichannel emissivity methods, and physically based emissivity methods. In this study, we used a Land Surface Emissivity (LSE) model proposed by Sobrino et al., [35] for LST retrieval and for a comparative purpose in Landsat 5, 7, and 8 images. Thus, accuracy difference based on different Landsat product is not to be of importance as for as LST is concerned. LSE proposed by Sobrino et al., [35] is considered the best model for all Landsat missions.

Farmers Perception Survey
A field survey of the rice crop was conducted in Larkana District during Kharif season 2017, aided by a map and Garmin 64 s GPS device. The field survey was conducted between September and October 2017. Before the field visit, a proforma for collecting the required geo-referenced data was prepared. The tabular proforma contained longitude, latitude, date, union council names, Taluka names, and district and crop name. The proforma was completed during the field visit using a Garmin 64 s GPS device for taking coordinates of the survey respondent locations. The perception survey of farmers was conducted to validate the sowing and harvesting shift in the rice crop plantation and the Sustainability 2020, 12, 3586 6 of 15 total rice crop grown in the district. The sample size (100) was calculated through a Cochran's formula because it is best when the population size is large and not known.
where: e is the desired level of precision (i.e., the margin of error), p is the (estimated) proportion of the population which has the attribute in question, q is (1 − p) [36]. We selected small, medium, and large farmers from all union councils of Larkana District for the farmer's perception survey. Thus, we can say it is a multistage random sampling. During the field survey, nearly all the union councils (UCs) of the district where rice crop production occurs were covered. The spatial distribution of the field data sampling points in Larkana District are shown in Figure 1.

Meteorological Data
Meteorological data from 1994-2017 in Larkana District was collected from the Pakistan Meteorological Department (PMD) in Karachi. The data consists of the mean monthly and annual minimum and maximum temperature in degrees centigrade ( • C), as well as the mean monthly and annual sum of precipitation in millimeters (mm). The crop season (June-October) maximum and minimum temperature and rainfall series were analyzed for trend analysis and changes in the trend slope. For trend detection, the Mann Kendall non-parametric test was used and the trend slope estimated through Sen's estimator. It is widely used to detect trends in time series data. It is a distribution-free test and does not require normally distributed data [37][38][39][40]. The N values time series are ranked from smallest to largest and computed by the Sen's estimator. Sen's estimator of the slope is the median of these N values.

Temperature and Rainfall
The annual temperature and annual rainfall data of Larkana District from 1994 to 2017 were analyzed on an annual basis as well as for the crop season months (June-October) to identify temporal changes in temperature and rainfall patterns. The crop season maximum and minimum temperature and the rainfall are graphically represented in the following Figure 2.  Trend analyses for monthly mean maximum and minimum temperature series reveal that minimum temperature series during the crop season shows an increasing trend (0.105 • C/Year). In contrast, the maximum temperature series do not have a significant trend, but the Sen's slope is 0.022 • C/Year. The mean precipitation trend for the crop months also does not have a significant trend in the series, while the Sen's slope shows a 0.087 mm/year increase annually during the crop season months. Salma et al. studied the different climatic zones of Pakistan and found that the southwest zone of Pakistan does not have a significant trend in temperature and rainfall patterns [41]. The tendency of air temperature to rise also influences the LST, as [42] suggests that air temperature alone explained 81-98% variance in LST under clear sky and cloudy conditions.

Quantification of the Change in Sowing and Harvesting Dates
Images for the months June, July, September, and October for the years 1994, 2000, 2010, and 2017 were classified to determine temporal variation in the area under rice crop in Larkana District and are graphically represented in Figure 3. For accuracy assessment of the crop mask for rice and non-rice areas, the results were compared with the previous study in the same district. Siyal et al. [24] Figure 3. The graph shows that the area under rice crop in June from 1994 to 2017 decreased due to the shift in sowing dates. The results reveal that there is a 106% decrease in the rice crop area in June from 1994 to 2017, this gap narrowed down to 14% in July. We analyzed the satellite imageries for these sowing months (June and July) with only one-day difference. On 20 June 1994, rice grown on a 42% area in the district as compared to 19 June 2017, image, only 16% completed the sowing. The scenario for July was the same when 52% completed the sowing process in 1994 as compared to only 35% in 2017. These stats clearly depict the shift of rice sowing pattern towards first and second quarter of July. Harvesting months September and October showed that harvesting started from the first quarter of October in 1994 as compared to 2017, in which harvesting started in the third quarter. This indicates that the area under rice crop starts reducing as compared to September in the

Sowing and Harvesting Dates Identified through Land Surface Temperature
To further validate the sowing and harvesting shift of the rice crop pattern, LST was analyzed for the selected satellite imageries. The results of the mean LST for the selected months for 1994, 2000, 2010, and 2017 was calculated from the satellite images shown in Table 3.
The LST has a strong correlation with vegetative growth on the field. The area under rice crop in different months and LST has a significant negative correlation, as shown in Table 3. Furthermore, the functional relationship between rice crop area and LST was established using a regression analysis method. The linear and nonlinear regression models between two variables were performed and the function of the relationship is shown in Figure 5. The logarithmic functional relationship between these two variables has the highest correlation coefficient (R 2 = 0.7318), as compared to the linear relationship (R 2 = 0.7119). Different studies also had similar findings [20,22,23]. It shows that the more vegetative growth on the land, the lower will the LST be, and vice versa. The impact of the shift in sowing and harvesting dates on the temporal scale influenced the LST of Larkana District for June and July (sowing period) and September and October (harvesting period) for the years 1994, 2000, 2010, and 2017. The results reveal that the mean LST of Larkana District has increased during June and July from the years 1994 to 2017. Furthermore, the mean LST has decreased during the harvesting period from 1994 to 2017. The LST of the study area for the targeted months and years are graphically shown in Figure 6. Mean LST ( • C) has been increasing in the sowing period of rice crop from 31.9 • C in June 1994 to 35.8 • C in June 2017, and from 32.8 • C in July 1994 to 36.8 • C in July 2017. Furthermore, the LST decreased during the harvesting period of rice crop from 31 • C in October 1994 to 28.6 • C in October 2017. This study found that LST is 11% more on 19 June 2017, as compared to 20 June 1994. This difference in temperature is due to a difference in the area under rice cover. The difference between the areas under rice crop also confirms from the difference in LST between 1994 to 2017, as shown in Figure 6, and the same trend is revealed in July. This temporal variation of LST also proves the shift of sowing of rice crop. The delay of rice crop harvesting also depicts from LST variation between October 1994 to 2017. As we have shown that there is more area under rice crop in October 2017 as compared to 1994, this will impact the LST to be 8% (28.6 • C) less in 2017 as compared to 1994. For the last 24 years, mean LST in June and July has increased 3.9 • C and 4 • C respectively, whereas in October the temperature decreased to 2.4 • C due to more area under crop as compared to 1994.

Crop/Socioeconomic Survey of Larkana District
The responses and results obtained from the crop/socioeconomic survey of Larkana District (i.e., Taluka Ratodero, Taluka Larkana, Taluka Bakrani, and Taluka Dokri) are discussed and graphically represented as follows: From all the responses in Table 4 during the rice crop sowing and harvesting survey of Larkana District, the study concluded that the sowing and harvesting dates of the rice crop in the district shifted 15-30 days during the last 25 years. Secondly, the responses regarding rainfall and temperature reveal that the amount of rainfall in the district has decreased, while the temperature has increased during the last three decades. 61% of the people responded that the temperature has increased in summer since the last 25 years, while 39% of the people had responded with mixed ideas. 6 Has rainfall increased/decreased in summer during the last 25 years?
75% of the people responded that the rainfall has increased in summer since the last 25 years, while 25% of the people had no idea about that.

Discussion
LST and its contributing factors in an irrigated environment have practical and theoretical implications for managing the agricultural production system better. Changing patterns of LST-due to multiple reasons, has a more significant influence on crop growth. Different studies find a significant negative correlation between NDVI and LST. The correlation coefficient we have found in the different months was stronger as compared to the reported results of Siyal et al. [23] and relatively weaker as compared to Guoming et al. [19]. Guoming et al. [20] reported that the correlation coefficient between paddy rice and LST was stronger (greater than 95%)-except for August, as compared to our study (less than 90%). The reason for this difference in statistical correlation between Guoming's findings and our study is probably due to study area features. Northern Sanjiang Plain area is a relatively cooler area as compared to Larkana. Different studies found that LST influencing factors are intricate and may depend upon the surface air temperature, wind speed, and land terrain, etc., [42][43][44]. It might be possible that the variation of this statistical correlation is due to the study area specific characteristics. The same is the case for study with Siyal et al. [23]; they found a lower value of correlation as compared to our study. The Indus delta has a different landscape, with a large area of tidal floodplain area and a relatively larger number of water bodies as compared to Larkana. These contributing factors of LST influence the relationship with NDVI.

Conclusions and Recommendations
The use of geospatial tools for the estimation of the shift in the sowing and harvesting dates of the crop is a new field in the Pakistani context. For this purpose, literature for the comparison purpose is a huge limitation we faced. This study is unique in the local context of Sindh and Pakistan. Based on the ground truth surveys, results obtained from Landsat satellite imageries, historical rice crop data, meteorological data, and crop/socioeconomic surveys, the following conclusions are drawn. On 20 June 1994, the rice was transplanted on 14.7% of the area of the district while it was only 7.1% of the area in 2017. Similarly, the area under rice crop in the first week of July 1994 was 18.3% compared to 8.15% during the same period in 2017. Thus, showing a delay/shift in the sowing of the rice crop. However, in the first week of October 2017, the rice crop was standing on 46.8% of the area while it was on 34.6% of the area in the year 1994 on the same date. This LST pattern depicts a delay in the harvest of the rice crop. Mean LST ( • C) has been increasing in the sowing period of rice crop from 31.9 ( • C) in June 1994 to 35.8 ( • C) in June 2017, and from 32.8 ( • C) in July 1994 to 36.8 ( • C) in July 2017. Furthermore, the LST decreased during the harvesting period of rice crop from 31 ( • C) in October 1994 to 28.6 ( • C) in October 2017. One of the possible reasons for the delay in sowing and harvesting dates is soil moisture anomalies and their relationship with LST [10]. This hypothesis needs to be further examined in future studies in this area. The study of LST is important to understand the impact on rice crop growth, crop tiller, and consequently, on yield. In future research, these aspects of being considered as control variables so that future cropland planning may account for these changes.
Based on the results, we recommend the remote sensing technique for the estimation of actual crop grown area, and it is more viable and gives a reliable estimation as compared to the crop reporting services. Canal commands of Pakistan's cultivated on the rotation basis; thus, late sowing of the one crop may influence the optimum sowing dates of the other crops. For the proper crop scheduling, there is a need to optimize the sowing dates of rice-wheat cropping pattern as an adaptation strategy of the future climate variability scenarios.