Back to the fields? Increased agricultural land greenness after a COVID-19 lockdown

The literature on the assessment of agricultural systems and productivity through the use of satellite data has proven the role that satellite imagery can play at a variety of scales (Jia et al 2019) and the depth of the insights that they can provide. This is especially true for developing countries (Lobell et al 2018), such as Indonesia (Yamamoto et al 2019). For instance, a large body of literature has shown that composite indexes based on multi-spectral satellite observations—such as the NDVI and EVI indexes—can be effective proxies to detect and forecast yield variation of a multitude of crops such as soybeans, rice, and wheat, even in smallholder agriculture, with results consistent with official yield statistics and correlation maps of historical yield data. Vegetation indexes are also well suited to study scattered agricultural topography and to inform management decisions on the optimal timing for planting and harvesting.

A recent review by Jia et al (2019) points out how advances in machine learning and the availability of data from earth observing satellites can dramatically improve the ability to map croplands over long periods and over large regions. This includes the capacity to link vegetation productivity to remotely sensed estimates of absorbed photosynthetically active radiation.

Lobell et al (2018) use data from a survey experiment in Uganda to evaluate the accuracy of Sentinel-2 remotely-sensed plot-level maize yields with respect to ground-based measures relying on farmer self-reporting. Similarly, Wolanin et al (2019) apply machine learning to a process-based model to estimate crop primary productivity from Sentinel-2 and Landsat 8 satellite data. They highlight the benefit of using the red edge bands from Sentinel-2 in increasing estimation accuracy. Cai et al (2019) integrate satellite and climate data to predict wheat yield in Australia using machine learning approaches. They highlight a generally very high estimation accuracy, with optimal prediction performance can be achieved two-month lead time before harvest. Burke and Lobell (2017) demonstrate the potential to track smallholder maize yield variation in western Kenya using a combination of 1-m Terra Bella imagery and intensive field sampling on thousands of fields over 2 years. They obtain R² coefficients of up to 0.4. Their method can also detect positive yield responses to fertilizer and hybrid seed inputs.

Kubitza et al (2020) review literature on the use of satellite data to assess adoption and impacts of agricultural management practices in developing countries. They find that in general satellite data are successful in detecting cropping intensity, tillage, crop residue cover, irrigation, and soil and water conservation. They highlight that only a limited number of studies has engaged with estimating yield impacts and the economic implications, irrespective of solid methodologies having been developed. Sjöström et al (2011) appraised the potential of MODIS EVI for modeling gross primary production across African Savannah ecosystems. In their study, EVI was found to correlate well with estimated GPP on a site-by-site basis.

Son et al (2014) carry out a comparative analysis of multitemporal MODIS EVI and NDVI indexes data for large-scale rice yield estimation in the Mekong River Delta. They find that the approach could be used to estimate rice production prior the harvesting period through the use of physical crop yield models, and that EVI-based models were more accurate than those from NDVI-based models. Their correlation coefficients (R²) range from 0.62 to 0.71 for spring-winter and 0.4 to 0.56 for summer-autumn rice crops, respectively. Figueiredo et al (2016) estimated correlation maps to assess historical yield during the soybean crop cycle from EVI data in Paraná State, Brazil. Their model was able to explain 96% to 98% of the variance in estimated yield from correlation maps.

Duveiller et al (2012) demonstrated how it is possible to monitor the growth and evolution of winter wheat over a decade using the green area index (GAI) even on a highly fragmented topography characterized by non-wheat areas within the spatial resolution of MODIS. Their results are congruent with ground trough data at both the regional and the filed level with a relative root mean square error (RRMSE) of 25.7% and 37.6%, respectively.

Jain et al (2017) developed a 30-meter resolution yield maps to explore the causes of yield gaps of smallholders in the Indo-Gangetic Plains, India. Using a crop model simulation, the authors were able to rely only on satellite vegetation indices without the need for any ground data. Finally using Random forest the authors proved that roughly 80% of the variance in yields is explained by differences in the weather conditions and management practices (such as the saw timing)Their results suggest that yield could be increased by 11 to 32% if more focus is given to the later.

The COVID-19 pandemic is exerting enormous pressure on the global economy in both higher and lower income countries. The observed labour market responses have however been remarkably different (Bartik et al 2020): in high income countries a significant share of the population is employed in the tertiary sector and could therefore move to home working without significant changes in the work productivity. Workers of more affected sectors could receive compensation and income support schemes from national governments, which overall mitigated the short-run household income plunge.

Conversely, in middle and lower income countries many people's incomes largely depend on agriculture, handcraft or tourism (Danquah et al 2020). The latter was one of the most affected sectors: according to recent estimates (WTTC 2020), up to 75 million tourism sector workers are at immediate job risk as a result of COVID-19, with a travel tourism GDP loss in 2020 of up to $2.1 trillion. Yet, lower income countries' governments often lack the means to support workers who lost their jobs as a consequence of COVID-19 limitations on social interactions. In these regions, populations were thus hit more harshly by the pandemic, and faced the necessity to adopt autonomous adaptation actions to try to compensate for the reduction or lack of their previous income (Ataguba 2020). Reports have shown that many tourist workers have decided to go back to their towns of origin to pursue agricultural activities and mitigate the labour demand shock in their previous employment sectors (UN World Food Programme 2020, Laula and Paddock 2020).

Badung is one of the 9 regencies in the island of Bali, Indonesia. It spans over 418.52 km² (figure 1) from north to south and consists of 6 districts and 46 villages (table 1. Like the rest of the island, Badung Regency relies on tourism as its main source of income but according to Indonesian Central Bureau of Statistics (BPS 2020), it is also the richest regency (the revenues in 2019 Badung were higher than the total revenues from other 8 regencies in Bali), the one with the highest minimum wage (IDR 2,9 million), the lowest percentage of people in poverty (1.78%) and the one that has experienced the highest annual population growth rate from 2010-2019 (2.36%). Based on recent statistics, around 140,000 people in Badung regency work in the tourism sector, and only 3,000 work in agricultural-related activities (BPS 2020).

Zoom In Zoom Out Reset image size

According to the Central Bureau of Statistics (BPS 2020), its two biggest districts are Petang and Kuta Selatan. The southern districts (Kuta, Kuta Selatan, Kuta Utara) are characterised by a strong tertiary sector that has been constantly growing together with the population living in the three districts. On the other hand, the northern districts (Mengwi, Abiansemal, Petang) are mostly agricultural areas where the population has remained overall steady in the last few years and tourism has grown slowly as compared with the south. The population of Kuta Selatan, the center of tourism activities of the island, grew by 4.4% in 2018-2019, while the population of Petang, the center of agricultural activities declined by −0.22% in the same period.

Badung, and in particular its northern districts, is also an important agriculture hotspot: 45% (190.49 km²) of its territory (418.52 km²) is devoted to agricultural activities (see figure appendix A.1). Most rice fields are located in the southern area of Badung (which has lower altitude). In contrast, non-paddy commodities are located in the northern area of the regency, where the most important horticulture commodities are coconut (201,840 ha), rice (12,943.13 ha), chayote (152 ha) and chilli (105 ha) (BPS 2020).

A tropical climate characterises Badung and more in general Bali, with monsoon from October to April (wet season). According to data from Meteorogical, Climatological, and Geophysical Office—Region III Denpasar (BPS 2020), January is the period with the most abundant precipitations (about 340 mm) and highest number of rainy days (18 days).

The analysis is based on remotely sensed environmental and atmospheric data detailed in table 2. The data includes the MODIS EVI and the variables used as predictors in the model. A time series is assembled based on data from the beginning of 2010 to 9th of July 2020. Raw data is processed in Google Earth Engine (Gorelick et al 2017), where each remotely-sensed variable is first masked over cropland pixels (as from the GFSAD1000 gridded cropland product (Teluguntla et al 2015); then, the mean value of each variable is extracted for each of the 6 districts of Badung for each day of the period considered.

The processed data are imported into the R scientific programming environment, where the data. Table (Dowle et al 2019), caret (Kuhn et al 2008), and Xgboost (Chen et al 2015) packages are used for statistical analysis.

Before the model training stage, missing data (19%) are imputed through a chained random forest Stekhoven and Buhlmann (2011) as implemented in the package MissRanger Mayer (2019) where each variable containing missing values are imputed by using all the other variables as predictors iteratively until the mean out-of-bag (OOB) error reach its minimum value. Finally, a Savitzky-Golay filter is applied to remove the noise that characterises satellite data (see Chen et al (2004)).

The technological advances and the need to account for the specific features of a given study area have brought researchers to develop a multitude of vegetation indexes (Viña et al 2011). For example, the EVI (Enhanced Vegetation Index, equation (1)) was developed to better account for atmospheric conditions and light reflection from the ground that causes several issues with the traditional NDVI (Normalized Difference Vegetation Index) measurement (equation (2)).

$$\begin{eqnarray}&&{EVI}=2.5\times {NIR}-{RED}\times ({NIR}+6\times {RED}-7.5\times {BLUE})+1\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{NDVI}=\displaystyle \frac{({NIR}-{RED})}{({NIR}+{RED})}\end{eqnarray} \tag{ 2 }$$

where NIR, RED, and BLUE represent the wavelength bands measured by a satellite.

In particular, the EVI includes additional wavelengths and adjustment parameters to correct for the inaccuracies of NDVI (Huete et al 2002). Variations in solar incidence angle, atmospheric conditions like distortions in the reflected light caused by the particles in the air, and signals from the ground cover below the vegetation are corrected for using EVI. Moreover, the NDVI is mostly chlorophyll sensitive, while the EVI is more responsive to canopy structural variations (Morcillo-Pallarés et al 2019).

3.3.1. Time-Series Prediction of EVI in a machine learning framework

3.3.2. Lockdown effect on EVI estimation

After model training and testing to generate the synthetic projection, we estimate the effect of the lockdown as the difference between the synthetically estimated time series and the actual data (Δ_{Observed,Counterfactual} ) before and after the lockdown. An analogous methodology is applied in Granella et al (2020) to estimate the effect the lockdown on pollution levels.

We subtract our prediction to the observed value of EVI and we estimate the difference-in-differences linear regression using ordinary least squares (OLS) as per equation (3):

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{Observed},{{Counterfactual}}_{i}t}={y}_{{it}}-{\rho }_{{it}}=\alpha +\beta {T}_{{it}}+{\epsilon }_{{it}}\end{eqnarray} \tag{ 3 }$$

where y_it is EVI measured on district i on day t, ρ_it is the predicted value, and T is a dummy equal to 1 during the lockdown (from 1^st April 2020 to the 9^th July 2020) and 0 prior to it (from 2^nd April 2019 to the 31^st March 2020). We weight each observation based on the cropland area (measured in square meters) for each district to account for differences in agricultural land surfaces.

Several consideration support our choice of assessing the lockdown impact starting from April 1^st 2020:

• It is the date of the official lockdown beginning in Bali.
• It roughly ( ± one week) coincided with the date when most EU countries, the United States, and Australia—the main origin countries of international tourists—put their lockdown policies into force (Hale et al 2021).
• Even for the few countries where the lockdown begun about two weeks earlier (e.g. Italy), many tourists were still in Bali when the lockdown had begun in their countries of origin.

We consider the formal conclusion of the lockdown of July 9^th, as after that date the Indonesian government allowed internal flights. Hence local tourists were able to come to Bali. Extending our analysis further would capture a mix of effects, namely the lockdown effect and the local tourism effect.

To ensure that our results are not driven by a different treatment than the lockdown, we run a placebo test. We estimate equation (3) from April 2nd 2019 to March 31st 2020, with T being equal to 1 from April 2nd 2019 to July 9th 2019 and 0 afterwards. To carry out the test we: (i) drop all the outcomes for treated observations after they receive treatment (the beginning of the lockdown); (ii) insert a phantom treatment event in the middle of the remaining data (we select one year before the lockdown); (iii) run the linear estimation model and evaluate the significance of the β coefficient from equation (3).

To further appraise the attribution of the measured effect to the lockdown, we test additional specifications of equation (3) where monthly, weekly, and daily fixed effects are included as additional regressors.

Table 3 presents the R-squared (R²), root mean square percentage error (RMSPE), and mean absolute percentage error (MAPE) for both the training dataset (from the 1st of January 2010 to 1st April 2019) and test dataset (from the 2nd of April 2019 to the 1st April 2020) for the whole study area as an average of the 6 districts. Figure 2 illustrates the average observed and predicted value of EVI for the whole regency of Badung before and after the lockdown.

Zoom In Zoom Out Reset image size

The results show that our model is able to closely learn the historical EVI variation over agricultural land in Badung based on our set of atmospheric and seasonal predictors, with a training accuracy of 95%. When testing the trained model on unseen data (the future, compared to the training data time span), the model performs relatively well, with a testing accuracy of 74%.

Table 4 reports the results of the effect estimation model of equation (3). We find that the lockdown has on average increased the value of EVI on the whole Badung regency by 0.003 52 which corresponds to an increase of 1.17%, from the average values that would have been observed had the epidemic not occurred. This result is highly statistical significant.

The results of additional specifications of equation (3) where monthly, weekly, and daily fixed effects are added as additional regressors are found in the table appendix A.1. The coefficient for the lockdown remains positive and highly significant, reinforcing our attribution of the observed increase in the EVI index to the lockdown.

To explore the heterogeneity effect of the lockdown on the north and southern districts we run separate regression models for the two areas. It must be remarked that in the study area, the South is dominated by tourism and the tertiary sector while the North is mostly agricultural. Figures 3–4 illustrate the average observed and predicted value of EVI for the two geographical areas before and after the lockdown.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the north the EVI has increased by 1.06% while no statistically significant effect is measured in the south. This is coherent with the hypothesis that workers formerly employed in the tourism sector in the south might have moved back to their region of origin in the north to work in the fields and partially compensate for the income lost from previous employment.

Table 5 and figure 5 report the result of the regression model based on equation (3) used to estimate the lockdown effect alongside the result of the placebo test in which the pseudo treatment is set one year before the real lockdown period.

Zoom In Zoom Out Reset image size

The results of the robustness check show that the placebo estimation is not statistically significant, while the true lockdown impact has a positive and highly significant impact on the measured EVI index over agricultural land in Badung. The finding thus allows as to further strengthen the causal interpretation of the analysis.