Analyze past land cover change
To understand how Ethiopia has changed due to population growth in the past few decades, you’ll use the Change Detection Wizard to compute land cover change from 1992 to 2018.
Explore land cover layers
First, you'll download the compressed .zip file that contains the data you'll use in this lesson.
- Go to the page for the lesson data, and click Download.
The data file might take some time to download, since it contains rather large raster files.
- Locate the downloaded ChangeInEthiopiaData.zip compressed folder on your computer, and move it to a location of your choice, such as the Documents folder. Right-click it to extract its content.
Depending on your web browser, you may be prompted to choose a file location before you begin the download. Most browsers download to your computer's Downloads folder by default.
Next, you’ll create an ArcGIS Pro project.
- Open ArcGIS Pro. If necessary, sign into your ArcGIS Online account.
If you don't have ArcGIS Pro or an ArcGIS account, you can sign up for an ArcGIS free trial.
This lesson was most recently tested for ArcGIS Pro 2.7. If you're using a different version of ArcGIS Pro, you may receive different results.
You'll create an project using the Map template.
- Under New, choose Map.
- In the Create a New Project window, for Name, type Change in Ethiopia. For Location, click Browse and choose a location of your choice. Click OK.
The new project opens. You'll now add the land cover data you downloaded.
- On the ribbon, on the Map tab, in the Layer group, click the Add Data button.
- In the Add Data window, browse to the location on your computer where your unzipped ChangeInEthiopiaData folder is located.
- Double-click ChangeInEthiopiaData to open it. While pressing the Ctrl key, click the Ethiopia_LandCover_1992.crf and Ethiopia_LandCover_2018.crf datasets to select them. Click OK.
The two raster layers are added to the map.
Each layer is a land cover dataset derived from the European Space Agency (ESA) Climate Change Initiative land cover maps of the world. For more information, see the Item Details page. The ESA created a global land cover map for every year from 1992 to 2018.
The top layer is the 2018 land cover map. Both layers show the following generalized land cover classes: cropland, forest, shrubland, grassland, surface water, urban areas, and bare soil. Notice that much of the bare soil is in northern Ethiopia (in the highlands), while southern Ethiopia is made up predominantly of shrubland. The agricultural region, symbolized in pink, is in the center of the country, surrounding the dense urban area of Addis Ababa in dark red.
You'll now compare the two layers.
- In the Contents pane, select the Ethiopia_LandCover_2018.crf layer.
The .crf extension on the layer tells you that the dataset is in Cloud Raster Format (CRF). This is Esri’s native raster format, which is optimized for writing and reading large files in a distributed processing and storage environment.
- On the ribbon, on the Appearance tab, in the Compare group, click Swipe.
- Starting at the top of the map, drag your pointer down to reveal the layer underneath. Go back and forth to compare the two layers.
The layer on top is the land cover from 2018, and the layer underneath is the land cover from 1992.
Notice that some parts of the country have visibly changed. For example, the capital city, Addis Ababa, is indicated by a group of pixels in dark red in the center of the country. You can see that it has expanded noticeably in the 26-year period.
Both layers are categorical raster datasets. Categorical raster data is raster data where each pixel has a value that is representative of a class or category. It is sometimes referred to as discrete data, thematic data, or discontinuous data, and it is most often used in GIS to represent land cover, land use, or other zonal information such as risk level. In this case, the categories represented are land cover types, such as cropland, forest, water, and urban.
- On the ribbon, on the Map tab, in the Navigate group, click Explore to exit the swipe mode.
- Save your project.
You created a project in ArcGIS Pro, added two raster datasets to it, and compared them visually. In the next section, you'll compute the change between the two layers.
Compute land cover change
You’ll now use the Change Detection Wizard to detect change throughout the country between 1992 and 2018, with a focus on change that is likely due to population growth.
- In the Contents pane, select the Ethiopia_LandCover_1992.crf layer.
- On the ribbon, on the Imagery tab, in the Analysis group, click the Change Detection button, and choose Change Detection Wizard.
The Change Detection Wizard pane appears.
- In the Change Detection Wizard, on the Configure pane, click the Change Detection Method drop-down menu to see the options available for change detection.
The Categorical Change option is used to identify change that has occurred between two thematic (or categorical) rasters, such as land cover or risk zone. The Pixel Value Change option is used to calculate the difference in pixel values between two continuous rasters, such as temperature rasters or multiband imagery. Finally, the Time Series Change option is used to identify the date of change in a time series of images.
The Categorical Change method is selected by default, because the raster layer that was selected in the Contents pane when you launched the wizard is categorical raster data.
- For From Raster, verify that the Ethiopia_LandCover_1992.crf layer is selected.
- For To Raster, select the Ethiopia_LandCover_2018.crf layer.
This means that the Ethiopia_LandCover_1992.crf layer will be compared to the Ethiopia_LandCover_2018.crf layer.
- Click Next.
On the Class Configuration pane, you can choose which type of filtering to perform, the classes to include in the analysis, and the rendering method for the results. You are interested in seeing only the areas that have changed, and only the changes that are likely due to population growth.
- For Filter Method, verify that Changed Only is selected.
- In the From Classes list, keep all classes selected.
- In the To Classes list, hover over the Urban category and click only.
Now, the Urban class is the only class selected in the list. However, urban growth is not the only class that may indicate change due to population growth. An expansion of cropland could also indicate population growth.
- Check the box next to Cropland to include this land cover type in the To Classes list as well.
In summary, you want to detect all the areas that have changed to the urban or cropland land cover types.
- For Transition class color method, keep the default option of Average.
This will determine how the output classes are rendered.
- Click the Preview button.
In the Contents pane, a Preview_ComputeChange layer is added. This layer is dynamically generated and is not saved to disk. You'll generate the permanent change layer later in the workflow.
- In the Contents pane, turn off Ethiopia_LandCover_2018.crf and Ethiopia_LandCover_1992.crf to better see the new layer.
- On the map, zoom in with the mouse wheel button to view the capital city of Addis Ababa.
Notice the clusters of pixels indicating areas of change.
- In the Contents pane, click Preview_ComputeChange layer to select it.
- On the map, click several pixels indicating change. In the Pop-up pane, confirm the type of change that has occurred.
It seems that most of the changes are from Cropland, Shrubland, Grassland, or Forest to Urban areas. The city has clearly expanded between 1992 and 2018.
- Zoom out and zoom back in to the town of Nejo, approximately 350 km west of Addis Ababa, between Mendi and Gimbi.
Notice the clusters of pixels that converted from Forest to Cropland.
- When you’re finished exploring, close the pop-up window.
- In the Change Detection Wizard, on the Class Configuration pane, click Next.
The Post-processing pane appears. You'll now save your output to disk.
- For Smoothing Neighborhood, keep None.
That parameter allows you to smooth your results for better visualization. In this case, you don’t want to smooth the results because you are interested in calculating the land cover area, and smoothing results would change pixel values.
- For Save Result As, verify that Raster Dataset is selected.
- For the Output Dataset, click Browse and double-click Folders and Change in Ethiopia. For Name, type Ethiopia_LandCoverChange_1992_2018.tif. Click Save.
- Click Run.
The change dataset is added to your map.
- On the Change Detection Wizard pane, click Finish.
- In the Contents pane, right-click the Ethiopia_LandCoverChange_1992_2018.tif layer and choose Zoom To Layer.
- Save your project.
Analyze the results
In the previous section, you generated a land cover change raster. You'll now explore the results in more depth and create a chart.
- In the Contents pane, right-click the Preview_ComputeChange layer and choose Remove, as you won't need it any longer.
If you didn’t click Finish to close the Change Detection Wizard pane, the Preview layer cannot be removed.
- In the Contents pane, right-click the Ethiopia_LandCoverChange_1992_2018.tif layer and choose Attribute Table.
The attribute table opens.
- Review the attribute table.
The Class_name field lists different transitions to the Cropland and Urban classes. The Count field indicates the total number of pixels in each category. The Area field indicates the total area this represents (in square meters). The area can be calculated because the dataset is in a projected coordinate system, with linear units of meters.
When computing areas, it is important to start from raster datasets that are in a projection that preserve areas, also known as equal area. In this case, the land cover layers use the Africa Albers Equal Area Conic projection.
- In the attribute table, right-click the Area field and choose Sort Descending.
The rows are now sorted according to area, where the transition with the largest area is listed first.
- Look at the first two rows in the table.
The first row, with the Class_name value Other, represents all the transitions that occurred from 1992 to 2018 that were not included in the analysis. In the second row, No Change represents the pixels that did not transition but stayed the same.
You don’t need these rows, so you will delete them from the table.
- In the attribute table, press the Shift key, and click the start of the both rows to select them. Press the Delete key.
- On the ribbon, on the Edit tab, in the Manage Edits group, click Save to save all the edits.
- In the attribute table, review the first few rows.
In the attribute table, the class transition with the largest area is now Shrubland to Cropland. According to your analysis, 5,537,592,079.12 m2 (or approximately 5,538 km2) of Shrubland was converted to Cropland from 1992 to 2018. Next, a sizeable amount of Forest was also converted to Cropland. This points to cropland expansion into natural vegetation to support population growth.
Two rows down, you can see that a large amount of Cropland was converted to Urban land cover. This is consistent with findings that rapid urban expansion is threatening the fertile agricultural land surrounding Addis Ababa (Deribew, 2020).
You'll now create a bar chart summarizing these results.
- In the Contents pane, right-click the Ethiopia_LandCoverChange_1992_2018.tif layer and choose Create Chart, and click Bar Chart.
The Chart Properties pane appears and a blank chart appears at the bottom of the project.
- In the Chart Properties pane, set the following options:
- For Category or Date, choose Class_From.
- For Aggregation, choose Sum.
- For Numeric field(s), click Select, and check Area. Click Apply.
- For Split by, choose Class_To.
The chart updates. The x axis shows the Class_From land cover types, and the y axis shows the area (in m2) of each category that transitioned to Cropland (light blue bars) or Urban (dark blue bar).
Next, you'll improve the appearance of the chart to match the symbology in the data.
- In the Chart Properties pane, open the Series tab. Select the Stacked series option.
Cropland and Urban bars are now stacked on top of each other.
- In the Series table, click the symbol for Cropland and select color properties.
- In the Color Editor window, set the Red, Green and Blue values to 247, 198 and 196, respectively. Set the Transparency to 0%.
This is the color of the Cropland category in the land cover data.
- Click OK to apply the color.
- Repeat the same steps to change the color of the Urban land cover bars, using Red, Green and Blue values of 175, 55, and 46, respectively.
This is the color of the Urban category in the land cover data.
- In the Chart Properties pane, click the General tab, and set the following options:
- For Chart title, type Cropland and Urban Growth in Ethiopia
- For X-axis title, type Original Class (1992).
- For Y-axis title, type Total Area (m2).
- For Legend title, type New Class (2018).
- Turn off Description.
The chart updates to its final appearance.
- Review the final chart.
You can see that much of the loss of Bare, Forest, Grassland, Shrubland, and Water categories went to the Cropland class. If you hover over the Cropland bar, you can see that around 485 million square meters (468 square kilometers) of cropland was converted to urban areas. Meanwhile, the largest contributor to cropland expansion was Shrubland, followed by Forest, and then Grassland. Based on the chart, it seems that population growth in Ethiopia between 1992 and 2018 has primarily contributed to a substantial increase in agricultural land use. Urban growth is present but secondary.
Analyze recent vegetation change
To understand how Ethiopia is being impacted by a massive locust invasion, you will use Landsat 8 satellite imagery to compare vegetation index values from before and after the start of the invasion, which began in Kenya in December 2019 and spread to surrounded countries in the months since.
Explore imagery layers
First, you’ll create a map within your project, and add the two Landsat 8 images to it.
- If they are still open, close the attribute table and chart.
- On the ribbon, on the Insert tab, in the Project group, click New Map.
A new map, Map1, is added to the project next to the first map.
- Make sure Map1 is selected.
- On the ribbon, on the Map tab, in the Layer group, click Add Data. Browse to the location where your unzipped ChangeInEthiopiaData folder is located.
- Press the Shift key and select the Landsat8_2019_10_15.tif and Landsat8_2020_11_18.tif datasets. Click OK.
The two Landsat 8 images are added to the map. The first layer is an image that was captured on October 15, 2019, before the locust invasion began. The second image was captured on November 18, 2020, after the invasion had swept through the region. The images cover the city of Addis Ababa and surrounding rural areas up to the border of the Aledeghi Wildlife Reserve.
You'll optimize the rendering of the two images and compare them.
- In the Contents pane, select the Landsat8_2020_11_18.tif layer.
- On the ribbon, on the Appearance tab, in the Rendering group, click the Symbology button.
The Symbology pane opens and the Primary symbology option is set to RGB. Landsat 8 multispectral imagery has 11 spectral bands originally, but 7 have been provided in the images you downloaded. The Red, Green, and Blue channels are currently set to bands 1 (coastal aerosol), 2 (blue), and 3 (green), respectively. This a default band combination, due to the original band order. You will change the symbology to view the image in natural color rendering, which is composed of bands 4 (red), 3 (green), and 2 (blue).
- Set the Red channel to sr_band4. Set the Green channel to sr_band3, and the Blue channel to sr_band2.
The layer updates in the map. Now, vegetation is displayed in green, bare soil in brown or brownish gray, water is blue or bluish gray, and urban areas are bright grey.
Clouds and cloud shadows on the images have been set as NoData values using the Quality Assessment (QA) band available with Landsat 8 surface reflectance. Those areas appear empty, and you might see the layer below at those locations.
- Repeat the step above for the Landsat8_2019_10_15.tif layer.
You'll use the swipe tool to compare the two images from before and after the locust invasion.
- In the Contents pane, click the Landsat8_2020_11_18.tif layer to select it.
- On the ribbon, on the Appearance tab, in the Compare group, click Swipe.
- Drag the pointer from top to bottom to peel off the top image and reveal the second image underneath.
You can see that there is significantly more vegetation in the 2019 image compared to the 2020 image.
- When you are done exploring, on the ribbon, on the Map tab, in the Navigate group, click Explore to exit the swipe mode.
Compute the pixel value difference
Now that you’ve visualized the difference between the two images, you’ll compute the difference in vegetation using the NDVI vegetation index in the Change Detection Wizard.
- On the ribbon, on the Imagery tab, in the Analysis group, click Change Detection, and choose Change Detection Wizard.
- For Change Detection Method, verify that Pixel Value Difference is selected.
This time, this option was selected by default because the imagery is continuous.
- For From Raster, select the 2019 imagery layer. For the To Raster, select the 2020 image.
- Click Next.
The Band Difference pane appears. It will allow you to choose several options specific to the pixel value change mode.
- For Difference Type, keep Absolute.
The Absolute difference is the mathematical difference between the pixel values from each image. The Relative difference, by comparison, accounts for the magnitude of the values being compared. In this case, the values of the vegetation index you will use are already normalized (they range from -1 to 1), so there is no need to use relative difference.
- For the Band Difference Method, select Band Index Difference.
This option allows you to first compute a band index on each image before performing the comparison. In this case, you’ll use the NDVI index, which is used to compare vegetation coverage.
- For the Band Index, keep the default value of NDVI.
Normalized Difference Vegetation Index (NDVI) is an index that is commonly used to assess the presence or absence of healthy green vegetation in imagery. It uses spectral reflectance information from the Red and Near Infrared (NIR) bands and computes a ratio with the following formula:
NDVI = (NIR – Red) / (NIR + Red)
You'll specify which bands of your images correspond to the NIR and Red spectral bands.
- Set the Near Infrared Band Index to 5 - sr_band5 for both images. Set the Red Band Index to 4 - sr_band4 for both images.
- Click Next.
The Change Detection Wizard calculates NDVI on both images and computes the difference between them, so it may take a minute to complete. The Preview_Mask layer is then added to the map. This layer shows the difference in NDVI values.
In the Change Detection Wizard, the Classify Difference pane has also appeared. It has a histogram showing the distribution of difference values between the two dates. Positive values indicate increased NDVI (that is, an increase in healthy vegetation), while negative values indicate a loss in NDVI (that is, a loss in healthy vegetation).
You'll change the symbology of that layer to better understand what it shows.
- In the Contents pane, click the Preview_Mask option.
The Symbology pane appears.
- In the Symbology pane, for Color scheme, expand the drop-down list and check the Show names box. Scroll down the color ramp options, and choose Yellow-Green-Blue (Continuous).
- Check the Invert box.
The Preview_Mask updates. Areas where the healthy vegetation has decreased are dark blue or medium blue. Areas where the healthy vegetation has increased are light yellow. A few areas of the layer do not have a color because the imagery has clouds and they are displayed as NoData.
- Close the Symbology pane.
- In the Classify Difference pane, in the Explore Differences histogram, drag the maximum handle arrow to approximately 0 so that only the negative values of the histogram are selected between the minimum and maximum handles.
The Preview_Mask layer updates in the map to show only the pixel values between the minimum (-1.4) and the maximum (0) values selected in the histogram. You can see that the vast majority of the values is below 0, this means that most areas have lost vegetation. A small loss in NDVI is expected between two dates, especially since the capture dates are more than one month (and one year) apart. You are interested in identifying areas of significant loss in NDVI.
- Drag the maximum handle to approximately -0.25.
The layer updates. Now, only the areas that experienced a loss in NDVI of 0.25 or more are displayed in the map. You'll consider this to be significant loss of NDVI.
- On the Classify Difference pane, click Add New Class.
The minimum and maximum values get added to the Classify Output table. This functionality allows you to extract and classify a specific range of values from the difference raster. Instead of computing the difference between two datasets, you can highlight the phenomenon you are interested in.
- In the Classify Output table, set the Output value to 1, Class Name to NDVI Loss, and the Color to red.
- Click the Preview button at the bottom of the pane.
The Preview_ClassifiedDifference layer is added to the map. The red pixels represent all the areas that have experienced a significant loss in NDVI.
- Click Next.
You'll now save your output.
- In the Post-processing pane, choose the following parameter values:
- For Smoothing Neighborhood, keep None.
- For Save Result As, keep Raster Dataset.
- For Output Dataset, click Browse, and double-click Folders and Change in Ethiopia. For Name, type NDVILoss_2019_2020.tif. Click Save.
Specifying the .tif extension determines the output format for the raster dataset will be a TIFF file. For the list of all supported raster formats, see the Raster file formats documentation.
- Click Run.
The new dataset is added to the map.
- In the Contents pane, turn off the Preview_ClassifiedDifference, Preview_Mask, Landsat8_2020_11_18.tif, and Landsat8_2019_10_15.tif layers, to only view the NDVI_Loss_2019_2020.tifand basemap layers.
If the Preview_ClassifiedDifference layer looks different from the final result, it’s because the preview layer is generated using raster functions, which compute the results dynamically using a resampled pixel size depending on the current display of the dataset.
To compare at the source resolution, in the Contents pane, right-click the Preview_ClassifiedDifference layer, and select Zoom to Source Resolution. At that resolution, the two layers are exactly the same.
- Save your project. Do not close the Change Detection Wizard.
Small changes in NDVI are expected from year to year, but a large loss of NDVI such as the one you identified can only be attributed to a disruptive event. It is likely that the locust invasion is the cause of the loss you see in the imagery. Locusts swarm vegetated areas, and cropland is particularly vulnerable. This has resulted in a devastating loss for millions of people.
Perform the same analysis again and again
Because this change detection workflow is completed using raster functions, you can save the output (and the output from the previous module) as a raster function template. You can then use the raster function template on other images, for a fast and repeatable analysis that can be used in multiple areas, or for different years.
Next, you'll create your NDVI comparison workflow as a new raster function template.
- In the Contents pane, turn off all layers except for the Landsat 8 images and the basemap layers.
- In the Change Detection Wizard, in the Post-processing pane, click the Save Result As drop-down menu and select Raster Function Template. Click Run.
The Raster Function Template editor window appears, prepopulated with the functions you used to run the analysis.
- Right-click the top Raster input and rename it From Raster. Rename the bottom Raster input To Raster.
- Double-click the top Band Arithmetic function to open its Properties pane. Click the Variables tab and check the box to set the From Raster parameter to IsPublic. Click OK.
This means the From Raster parameter will be visible in the final raster function. When running the function, it will be possible to choose what raster should be used as the From Raster parameter.
- Repeat the same step for the bottom Band Arithmetic function, checking the IsPublic box for the To Raster parameter.
- In the Raster Function Template1 window, click Save.
- In the Save window, enter the following values:
- For Name, type Landsat 8 NDVI Loss.
- For Category, choose Custom.
- For Description, type Compares two Landsat 8 images and extracts a loss in NDVI of 0.25 or more.
- Keep all other defaults, and click OK.
- Close the Raster Function Template1 window. If you get a message asking you to save the edited function chain, click No.
- In the Change Detection Wizard, click Finish.
You'll now test the function that you just created.
- On the ribbon, on the Imagery tab, in the Analysis group, click the Raster Functions button.
The Raster Functions pane appears.
- In the Raster Functions pane, click the Custom tab, and expand the Custom1 category.
- Click Landsat 8 NDVI Loss to open the function you just created.
- In the Landsat 8 NDVI LossProperties pane, for From Raster, select the Landsat8_2019_10_15.tif layer. For To Raster, select the Landsat8_2020_11_18.tif layer.
- Click Create new layer.
The resulting function raster layer is added to the map. It used the exact same processing steps that you used in the Change Detection Wizard. You can provide any two Landsat 8 images to generate a similar result.
If you wanted to choose a different sensor (that is another imagery type), you must change the band index values in the Band Arithmetic functions to ensure the correct bands are used for NDVI calculation.
- Save your project, and close ArcGIS Pro.
In this module, you used the Change Detection Wizard to compute the pixel value difference between Landsat imagery layers and identify areas with vegetation loss. You also created a reusable raster function template to apply the same analysis to other data.
Perform change analysis where you live (optional)
You can use the Change Detection Wizard to compute the difference between two rasters from your own data collections. You can also use the wizard to compare two layers from an global image service, with the option of choosing to analyze any chosen location, including the area where you live.
Analyze change with the Global Land Cover image service
The Global Land Cover 1992-2018 image service is hosted on ArcGIS Living Atlas and can be configured for specific years to see how land cover has changed in the area where you live. As a stretch exercise, you can use the steps below to see how your region of interest is changing.
- Launch ArcGIS Pro and create a project using the Map template.
- If prompted, sign in to your ArcGIS Online organization account.
You'll access a dynamic layer that contains information like the one you've used, but for the whole world.
- On the ribbon, on the View tab, in the Windows group, click Catalog Pane.
- In the Catalog pane, click the Portal tab, and click Living Atlas. Search for Global Land Cover.
- Right-click the Global Land Cover 1992-2018 image service, and choose Add To Current Map.
- On the ribbon, on the Map tab, in the Inquiry group, choose Locate. In the Locate pane, enter a location of interest. For example, type Denver, Colorado, USA.
- When you have found your area of interest, close the Locate pane.
You can zoom out a bit until you see the whole area.
- In the Contents pane, right-click the Global Land Cover 1992-2018 layer and choose Copy. Right-click the Map heading, and select Paste.
A duplicate layer appears in the Contents pane. You'll customize each layer to show the 1992 and 2018 data only.
- Rename the first layer Land Cover 1992. Rename the second layer Land Cover 2018.
- In the Contents pane, right-click the Land Cover 1992 layer and choose Properties.
- In the Layer Properties window, click the Definition Query tab.
- Click New definition query and enter the definition query: where Start Dateis equal to 1/1/1992.
- Click Apply, and click OK to close the Layer Properties window.
- Repeat the same steps for the 2018 layer, setting the definition query to: where Start Date is equal to 1/1/2018.
You now have two global land cover maps in your project: one showing land cover from 1992, one from 2018. You'll now limit the processing extent to your area of interest.
- On the ribbon, on the Analysis tab, in the Geoprocessing group, click Environments.
- In the Environments window, for Processing Extent, choose Current Display Extent to limit the processing area to your region of interest. Click OK.
This step is important because you are limited in the number of pixels you can export from an analysis on an image service.
- On the ribbon, on the Imagery tab, in the Analysis group, click Change Detection, and Change Detection Wizard.
- In the Change Detection Wizard, select the following parameter values:
- For Change Detection Method, keep Categorical Change.
- For From Raster, choose Land Cover 1992.
- For To Raster, choose Land Cover 2018.
You can now proceed to analyzing land cover change in your region of interest.
The global land cover service is projected in the Web Mercator Auxillary coordinate system, which is not an area-preserving projection. To perform area calculation, you can download the data for your study area, and choose to reproject it to an equal area projection, such as Albers Equal Area Conic. To do that, in the Contents pane, right-click the global land cover layer, and choose Data > Export Raster.
In this lesson, you compared land cover datasets from 1992 and 2018 to calculate change due to urban growth and agricultural expansion in Ethiopia. Your findings indicate that population growth in Ethiopia has mostly contributed to an increase in agricultural land use, and urban growth is secondary. Although agriculture has expanded, a massive locust invasion in 2020 has resulted in large losses in crops, indicated by a drop in NDVI. You compared two Landsat 8 images and extracted regions that lost NDVI.