NDVI Values Analysis
- Vecihi Asım Kartoğlu
- Jan 11, 2022
- 4 min read
Developing technology benefits the agricultural sector by using satellite images. Everything from irrigation, soil analysis, harvest analysis to pest detection can be improved by processing satellite images.
One of the most preferred methods for observing, reporting and revealing statistical data on plant class is to use the Normalized Difference Vegetation Index (NDVI) Normalized Vegetation Difference Index.
NDVI is an abbreviation of English Normalized Difference Vegetation Index. In short, we can say that the life indicator and analysis of plants.
What is NDVI? Technically speaking, green leafy plants absorb sunlight, which we call solar radiation, which they use as an energy source while photosynthesis (the chemical event that allows plants to make organic food with light energy is called photosynthesis). Therefore, living plants appear brighter in near-infrared (NIR) cameras. Chlorophyll in plant leaves strongly absorbs visible light (VIS) (0.4-0.7 µm) for use in photosynthesis. The special cell structure of the leaves strongly reflects near infrared light (0.7 to 1.1 µm).
As a result, the more leaves a plant has, the more the wavelengths of light are affected. If we say how to calculate NDVI in the light of this information;
NDVI = (NIR-Red)/(NIR+Red) (For Landsat OLI 8-9)
Thanks to the NDVI technology, you can learn the status of the crops in the piece of land you own (field, land, etc.), by analyzing the NDVI images, you can learn whether the vegetables you grow or the vitality of the plants in this area or the plant disease.
NDVI = -1 to 0 represent Water bodies
NDVI = -0.1 to 0.1 represent Barren rocks, sand, or snow
NDVI = 0.2 to 0.5 represent Shrubs and grasslands or senescing crops
NDVI = 0.6 to 1.0 represent Dense vegetation or tropical rainforest
A project has been prepared in order to create NDVI images and then show the minimum, maximum, mean, median, standard deviation values in a graph.
First of all, satellite images were downloaded via the USGS EarthExplorer website in order to download satellite images.
Satellite images of 10 different times were selected.
Satellite images belong to a single region and cover the borders of Istanbul/Turkey.
Landsat 8-9 OLI-TIRS C2 L2 images are preferred.
Land Cloud Cover is below 5%.
To create NDVI, Band 4-Red, Band 5-NIR are sufficient.
The information for all images is as follows.
2018/04/23, Cloud=0.01
2019/03/25, Cloud=0.05
2019/05/28, Cloud=0.66
2019/09/17, Cloud=0.50
2019/10/03, Cloud=0.11
2020/04/12, Cloud=0.73
2020/07/01, Cloud=0.05
2021/02/10, Cloud=2.53
2021/03/14, Cloud=1.48
2021/08/05, Cloud=0.01
After the satellite images are downloaded, the desired image can be created by using the NDVI formula in the QGIS software, using the Python code or using the auxiliary plugins.
I would like to show you the process of creating a new image with the NDVI formula with the Raster Calculator Plugin in QGIS software.
After performing the same process on all images, clipping was performed due to the fact that a very large area of the satellite images is the sea.
Clip Raster by Extend function in QGIS software is used for this process.
After all images were brought to the same format, a shapefile was created in the polygon feature type surrounding the region to be calculated in order to calculate the zonal statistic values.
It is sufficient to create a single shapefile. Because all satellite images belong to the same region.
Then, the python code was run to generate the minimum, maximum, mean, median and standard deviation values of the images. And a total of 5 different values for 10 images are saved in the .txt file.
Line graph NDVI Time Series data is drawn from txt file.
The operations described in this step are carried out with python code.
The code is shared below.
Operations were done via the Python Console in QGIS software.
import gdal
import ogr
import os
import numpy as np
import csv
from datetime import datetime
import matplotlib.pyplot as plt
def boundingBoxToOffsets(bbox, geot):
col1 = int((bbox[0] - geot[0]) / geot[1])
col2 = int((bbox[1] - geot[0]) / geot[1]) + 1
row1 = int((bbox[3] - geot[3]) / geot[5])
row2 = int((bbox[2] - geot[3]) / geot[5]) + 1
return [row1, row2, col1, col2]
def geotFromOffsets(row_offset, col_offset, geot):
new_geot = [
geot[0] + (col_offset * geot[1]),
geot[1],
0.0,
geot[3] + (row_offset * geot[5]),
0.0,
geot[5]
]
return new_geot
def setFeatureStats(min, max, mean, median, sd, sum, count, names=["min", "max", "mean", "median", "sd", "sum", "count"]):
featstats = {
names[0]: min,
names[1]: max,
names[2]: mean,
names[3]: median,
names[4]: sd,
names[5]: sum,
names[6]: count,
}
return featstats
mem_driver = ogr.GetDriverByName("Memory")
mem_driver_gdal = gdal.GetDriverByName("MEM")
shp_name = "temp"
inputPathTif = "F:/geo-insight/NDVI/datas/inputs/"
for images in os.listdir(inputPathTif):
fn_raster = os.path.join(inputPathTif, images)
r_ds = gdal.Open(fn_raster)
p_ds = ogr.Open("F:/geo-insight/NDVI/datas/polygon.shp")
lyr = p_ds.GetLayer()
geot = r_ds.GetGeoTransform()
nodata = r_ds.GetRasterBand(1).GetNoDataValue()
zstats = []
p_feat = lyr.GetNextFeature()
while p_feat:
if p_feat.GetGeometryRef() is not None:
if os.path.exists(shp_name):
mem_driver.DeleteDataSource(shp_name)
tp_ds = mem_driver.CreateDataSource(shp_name)
tp_lyr = tp_ds.CreateLayer('polygons', None, ogr.wkbPolygon)
tp_lyr.CreateFeature(p_feat.Clone())
offsets = boundingBoxToOffsets(p_feat.GetGeometryRef().GetEnvelope(),\
geot)
new_geot = geotFromOffsets(offsets[0], offsets[2], geot)
tr_ds = mem_driver_gdal.Create(\
"", \
offsets[3] - offsets[2], \
offsets[1] - offsets[0], \
1, \
gdal.GDT_Byte)
tr_ds.SetGeoTransform(new_geot)
gdal.RasterizeLayer(tr_ds, [1], tp_lyr, burn_values=[1])
tr_array = tr_ds.ReadAsArray()
r_array = r_ds.GetRasterBand(1).ReadAsArray(\
offsets[2],\
offsets[0],\
offsets[3] - offsets[2],\
offsets[1] - offsets[0])
if r_array is not None:
maskarray = np.ma.MaskedArray(\
r_array,\
maskarray=np.logical_or(r_array==nodata, np.logical_not(tr_array)))
if maskarray is not None:
zstats.append(setFeatureStats(\
maskarray.min(),\
maskarray.max(),\
maskarray.mean(),\
np.ma.median(maskarray),\
maskarray.std(),\
maskarray.sum(),\
maskarray.count()))
else:
zstats.append(setFeatureStats(\
nodata,\
nodata,\
nodata,\
nodata,\
nodata,\
nodata,\
nodata))
else:
zstats.append(setFeatureStats(\
nodata,\
nodata,\
nodata,\
nodata,\
nodata,\
nodata,\
nodata))
tp_ds = None
tp_lyr = None
tr_ds = None
dict = zstats[0]
dict_zstats = dict.get('min'),dict.get('max'),dict.get('mean'),dict.get('median'),dict.get('sd')
p_feat = lyr.GetNextFeature()
with open('F:/geo-insight/NDVI/datas/output/zonal_stats.txt', 'a+') as txtfile:
txtfile.write("{} ".format(fn_raster[33:-4]) + " ".join(map(str, dict_zstats)) + "\n")
txtfile.close()
print("Txt file written: {}".format(dict_zstats))
x = []
a = []
b = []
c = []
d = []
e = []
for line in open('F:/geo-insight/NDVI/datas/output/zonal_stats.txt', 'r'):
lines = line.split()
x.append(lines[0])
a.append(float(lines[1]))
b.append((float(lines[2])))
c.append((float(lines[3])))
d.append((float(lines[4])))
e.append((float(lines[5])))
plt.plot(x, a, color='r', marker='o', linestyle='dashed', linewidth = 2, markersize=4, label='min')
plt.plot(x, b, color='g', marker='o', linestyle='dashed', linewidth = 2, markersize=4, label='max')
plt.plot(x, c, color='b', marker='o', linestyle='dashed', linewidth = 2, markersize=4, label='mean')
plt.plot(x, d, color='c', marker='o', linestyle='dashed', linewidth = 2, markersize=4, label='median')
plt.plot(x, e, color='m', marker='o', linestyle='dashed', linewidth = 2, markersize=4, label='std')
plt.xticks(rotation = 20, size=10)
plt.xlabel("Date")
plt.ylabel("Values")
plt.title("NDVI Time Series")
plt.legend(['min', 'max', 'mean','median','std'], loc='upper right')
plt.savefig("F:/geo-insight/NDVI/datas/output/NdviGraph.png")The values printed to the txt file are given below.
The created graphic has been shared.
Different NDVI images of the region
2021/08/05
2020/07/01
This project was created to evaluate and demonstrate the benefits of analyzing vegetations through satellite images.
It would be useful to use statistical values at different times.
Thanks
Geomatics Insight
by Vecihi Asım Kartoğlu
Comments