1. Libraries: The script uses rasterio for reading and writing raster data, and numpy for numerical operations. Install these.

2. Input Directories: Organize your NIR (Near-Infrared) and Red channel images in separate directories. For example:

   • NIR folder: Contains NIR channel images (e.g., LC08_L1TP_190045_20200801_20200807_01_T1_B5.TIF).

   • Red folder: Contains Red channel images (e.g., LC08_L1TP_190045_20200801_20200807_01_T1_B4.TIF).

3. Output Directory: Designate or create a directory where the NDVI output images will be stored.

4. Naming Convention: Ensure that the NIR and Red images for the same area and date have similar naming conventions to correctly pair them in the script.

   SCRIPT:

import rasterio

import numpy as np

import os

nir_folder = 'N:/Your NIR holder'

red_folder = 'N:/Your RED holder'

output_folder = ''N:/Your NDVI output holder''

nir_files = {f.split("_SR_B5")[0]: f for f in os.listdir(nir_folder) if f.upper().endswith('B5.TIF')}

red_files = {f.split("_SR_B4")[0]: f for f in os.listdir(red_folder) if f.upper().endswith('B4.TIF')}

print("Reading...")

print("Find NIR-files:", len(nir_files))

print("Find Red-files:", len(red_files))

for common_part in nir_files:

if common_part in red_files:

print(f"Couple processing: {common_part}")

path_to_nir = os.path.join(nir_folder, nir_files[common_part])

path_to_red = os.path.join(red_folder, red_files[common_part])

print(f"Open NIR: {path_to_nir}")

print(f"Open Red: {path_to_red}")

with rasterio.open(path_to_nir) as nir_src:

nir = nir_src.read(1)

with rasterio.open(path_to_red) as red_src:

red = red_src.read(1)

# Calculating NDVI by adding a small value to the denominator to avoid division by zero

ndvi = (nir.astype(float) - red.astype(float)) / (nir + red + 0.0001)

# Set NDVI to NaN where NIR + RED = 0

ndvi[(nir + red) == 0] = np.nan

ndvi_image_path = os.path.join(output_folder, f'{common_part}_NDVI.TIF')

print(f"Save NDVI: {ndvi_image_path}")

with rasterio.open(

ndvi_image_path,

'w',

driver='GTiff',

height=nir_src.height,

width=nir_src.width,

count=1,

dtype=rasterio.float32,

crs=nir_src.crs,

transform=nir_src.transform,

) as ndvi_dst:

ndvi_dst.write(ndvi, 1)

print("Done!")You can also use my .txt file to symbology your results (note that they are in Russian):# QGIS Generated Color Map Export File

INTERPOLATION:INTERPOLATED

0,0,0,0,0,255,Neutral/Bare Ground

0.20000000000000001,217,239,139,255, Medium vegetation

0.40000000000000002,166,217,106,255,Healthy vegetation

0.59999999999999998,102,189,99,255,High vegetation

0.80000000000000004,26,152,80,255,Very high vegetation

1.0,0,104,55,255, Maximum vegetation