How To Analyse NDVI Trend

Analyzing the annual maximum NDVI trend using time-series Landsat images

Photo by Håkon Grimstad on Unsplash

The normalized difference vegetation index (NDVI) is a simple graphical indication that may be used to analyse remote sensing data, often from a space platform, to determine whether or not the objective under observation has living green vegetation.

Satellite data of the Normalized Difference Vegetation Index can be used to quantify changing patterns in ecosystem productivity (NDVI). The estimate of trends using NDVI time series, on the other hand, varies significantly depending on the studied satellite dataset, the related spatiotemporal resolution, and the statistical approach used. We analyze the efficacy of a variety of trend estimate approaches and show that performance diminishes as inter-annual variability in the NDVI time series increases.

ArcGIS

The Environmental Systems Research Institute maintains ArcGIS, a geographic information system for dealing with maps and geographic information. It is used to create and utilize maps, collect geographic data, analyze mapped information, distribute and find geographic information, use maps and geographic information in a variety of applications, and manage geographic data in a database.

products/services from ArcGIS (Author Created)

I utilize Arcpy.ia and the ArcGis API for Python in this scenario. arcpy.ia is a Python module for managing and processing images and raster data. The module also provides all of the geoprocessing methods given by the ArcGIS Image Analyst extension, as well as enhanced functions and classes that allow you to automate your raster processing workflows.

Without any further ado, let’s get started. I’ll describe the procedure as I explain the code. First, I import the required modules into my workflow.

import os
import matplotlib.pyplot as plt
import arcpy
arcpy.CheckOutExtension("ImageAnalyst")

Create Raster Collection from time series Landsat images

Get a list of all the Landsat scenes in the folder.

I set up my data folder. I collected this data from U.S. Geological Survey, and you can get similar data by registering with U.S. Geological Survey. I also created four arrays for future use.

U.S. Geological Survey website

folder = r"\\landsat\\LE70820552011359EDC00"
raster_lists = []
acquisition_time = []
names = []
quality_bands = []

In the code below, I retrieve all tif files from the source folder and add them to the raster dataset. I also include the acquisition time and filename for each array. However, before adding the filename, I change it with .xml.

for directory in sorted(os.listdir(folder)):
for file in os.listdir(folder+"\\"+directory):
if file.endswith("xml") and len(file.split(".")) == 2:
raster_lists.append(folder+"\\"+directory+"\\"+file)
time = file.split("_")[3]
acquisition_time.append(time[:4]+"-"+time[4:6]+"-"+time[6:])
names.append(file.replace(".xml", ""))
quality_bands.append(folder+"\\"+directory+"\\"+file.replace(".xml", "") + "_pixel_qa.tif")

Create Raster Collection

To generate a raster collection, I utilize the RasterCollection method. The RasterCollection object easily sorts and filters a bunch of rasters and prepares a collection for subsequent processing and analysis. The RasterCollection object has six methods for generating statistics for each pixel across matching bands inside the collection’s rasters (max, min, median, mean, majority, and sum).

rc = arcpy.ia.RasterCollection(raster_lists, {"Varieble": "Landsat7", names, "STDTime": acquisition_time, "Quantity": quality_bands})

Filter

The dataset we’re utilizing has data covering 18 years. So I’m just going to use only 10 years of data. We may accomplish this by using the filterByCalanderRange function.

filtered_rc = rc.filterByCalenderRange(calender_field = "YEAR", start = 1997, end = 2007)

Visualize Raster Collection

No, we can view the Raster Collection.

filtered_rc

The result is as follows. As we can see, it’s displaying variables, names, and so on.

result (Author Created)

Visual Raster in item

The code below retrieves the raster from the collection’s final entry. I use the arcpy.ia.Render function to render the result and the exportImage method to display it as an image.

To improve symbology, use the Render function to adjust the appearance of a raster object. When working in a Jupyter notebook, where data presentation is a crucial advantage, this function comes in useful. The function generates a raster object using the specified rendering rule or color map. At least one rendering rule or color map must be given.

first_item = filtered_rc[0]['Raster']
rendered_raster = arcpy.ia.Render(first_item_raster, rendering_rule = {"bands": [3, 2, 1], "min": 500, "max": 1000, colormap='NDVI'})
rendered_raster.exportImage(width = 500)

This is the resulting picture. As we can see, this image has a lot of clouds. As a result, we must clean up this image.

Result image (Author Created)

Pre-processing

In the Pre-processing phase, I’m going to process the previously obtained image and remove clouds, shadows, water, and snow.

In the code below, I first retrieve the raster and the quality of each raster. Then I save the pixel values that I need to remove from the raster. Following that, I use arcpy.ia.Apply with bandRaster to apply to each and every raster. bandRaster will mask the pixels in the specified RGB bands. Then I assign a cleaned rater and assign it to return.

def clean(item):
raster - item["Raster"]
qa_raster = arcpy.Raster(item["Quantity"])
values_tb = (qa_raster != 68) & (qa_raster != 72) & (qa_raster != 80) & (qa_raster != 96) & (qa_raster != 132) & (qa_raster != 136) & (qa_raster != 112) & (qa_raster != 144) & (qa_raster != 160) & (qa_raster != 176)
masked_clean_band_list = []
for bandRaster in raster.getRasterBands([1, 2, 3, 4, 5, 6]):
masked_clean_band = arcpy.ia.Apply((values_tb, bandRaster, 0), "Local", args = {"Operation": "StdTime"})
masked_clean_band_list.append(masked_clean_band)
masked_clean_raster = arcpy.ia.CompositeBand(masked_clean_band_list)
return {"raster": masked_clean_raster, "Name" : item["Name"], "StdTime" : item["StdTime"]}
cleaned_rc = filtered_rc.map(clean)

visualize

Using this code, I retrieve the raster within the collection’s second item.

first_item_raster = cleaned_rc[0]['Raster']
rendered_raster = arcpy.ia.Render(first_item_raster, rendering_rule = {"bands": [3, 2, 1], "min": 500, "max": 1000, colormap='NDVI'})
renderes_raster.exportImage(width = 500)

As shown in the image below, the whole water area has been masked out.

Result image (Author Created)

Calculate NDVI

I define a method called clacNDVI to calculate NDVI, and I use arcpy.ia.NDVI to create an NDVI raster from an item raster and two band values. Finally, I return parameters like item, name, and standard time.

def clacNDVI(item):
ndvi = arcpy.ia.NDVI(item['Raster'], nir_band_id = 4, red_band = 3)
return {"raster": ndvi, "Name": item["Name"], "StdTime": item["StdTime"]}
ndvi_rc = cleaned_rc.map(calcNDVI)

When we call the function, we get the NDVI raster.

Visualize

Let’s visualize the raster in the same way that I did earlier.

first_item_raster = ndvi_rc[0]['Raster']
rendered_raster = arcpy.ia.Render(first_item_raster, rendering_rule = {"bands": [3, 2, 1], "min": 0, "max": 1.0, colormap='NDVI'})
rendred_raster.exportImage(width = 500)

It’s the same raster as the one in previous. However, it is now displaying NDVI.

Result Image (Author Created)

Get Max NDVI per year

I converted the multidimensional raster collection using the code below. This method returns a multidimensional raster dataset, with each item in the raster collection representing a slice in the multidimensional raster.

anual_max_ndvi_mdraster = arcpy.ia.Aggregate(ndvi_mdraster, demension_name = "StdTime", raster_function = "MaxIgnoreNoData", aggregate_definition = {"interval": "yearly"})

I use the arcpy.ia.Aggregate function to determine the maximum NDVI each year. This method generates a raster object based on the provided multidimensional raster. Aggregation is computed by default for all variables linked with the chosen dimension. In this case, I’m using STDTime as the dimension time and MaxIgnoreNoData as the raster function with an annual interval time.

Generate Trend

Now I apply linear trend for the annual maximum multidimensional raster. The linear trend line is a best-fit straight line for estimating basic linear relationships. A linear trend denotes a pace of change that increases or decreases at a consistent rate.

max_ndvi_linear_trend = arcpy.ia.Apply(Annual_max_ndvi_mdraster, "TrendAnalysis", raster_function_arguments = {"DimensionName": "StdTime", "RegressionType": 0})

Visualize

We can retrieve the slope of the trend raster using the code below. getRasterBands returns a Raster object for each band supplied in a multiband raster dataset. In this situation, I use 1 as the raster band.

slope - max_ndvi_linear_trend.getRaterBands(1)

Remap

To modify negative pixel values, I utilize the remap function. The remap function allows you to alter or reclassify raster data pixel values. This may be accomplished by either supplying a range of pixel values to map to an output pixel value or by mapping the pixel values to output pixel values using a table. As a result, the output will be a boolean map.

increase_decrease_trend_map = arcpy.ia.Remap(slope, input_ranges=[-10000, 0, 0, 10000], output_values=[-1, 1])

Render

I use the render function with different colors to render (green and gray).

rendered_raster = arcpy.ia.Render(increase_decrease_trend_map, colormap = {"values":[0, 1], "colors": ["#B0C4DE", "green", "gray"]})
rendered_raster.exportImage(width = 500)

When we look at this image, we can see that the NDVI is increasing in the green area. The grey spots on the image represent NDVI-decreasing areas.

Result Image (Author Created)

Let’s recap

NDVI is a pixel-by-pixel mathematical computation performed on an image with GIS tools. It is a plant health indicator determined by measuring the absorption and reflection of red and near-infrared light. NDVI may be used to study land all over the world, making it suitable for both focused field studies and continental or global-scale vegetation monitoring. ArcGIS, a geographic information system for dealing with maps and geographic information, is maintained by the Environmental Systems Research Institute. ArcGIS can be used to calculate NDVI.

Reference

What is NDVI (Normalized Difference Vegetation Index)?

ArcGis Documentation

ravindu kavishwara - Machine Learning Researcher - Self-employed | LinkedIn