RemoteSens_INdices is written in python and used to calculate different indices like NDVI, NDWI, NDMI, SAVI, EVI etc. for satellite imagery data.
Before running this code in your system you have to install numpy
and rasterio
library in your system. To install the numpy
library follow this link.
Geographic information systems use GeoTIFF and other formats to organize and store gridded raster datasets such as satellite imagery and terrain models. rasterio
reads and writes these formats and provides a Python API based on Numpy
N-dimensional arrays and GeoJSON
. To install the rasterio
, clik here.
Now we will breakdown and see how to execute the code:
import numpy as np
import rasterio
These are import statements. The numpy
library is imported with the alias np
, and the rasterio
library is imported. These libraries will be used to perform mathematical operations and read/write geospatial raster data, respectively.
input_file = r'E:/Python_Tutorial/Indices/Sentinel-2_Bardhaman.tif'
This line defines the path to the input raster file, which is a Sentinel-2 image. You replace the path as per your desired location.
ndvi_file = r'E:/Python_Tutorial/Indices/ndvi.tif'
ndwi_file = r'E:/Python_Tutorial/Indices/ndwi.tif'
ndmi_file = r'E:/Python_Tutorial/Indices/ndmi.tif'
evi_file = r'E:/Python_Tutorial/Indices/evi.tif'
savi_file =r'E:/Python_Tutorial/Indices/savi.tif'
These lines define the output file paths for the various indices that will be calculated. Five indices are being calculated: NDVI, NDWI, NDMI, EVI, and SAVI. For each index, a separate output file is created.
with rasterio.open(input_file) as src:
data = src.read()
profile = src.profile
profile.update(count=1)
Opening the input file using rasterio
, reading the data, and storing it in a numpy
array data. The metadata of the input file is stored in a dictionary called profile
, and the count
parameter is updated to 1 to indicate that the output files will have a single band.
red = data[2].astype(np.float32)
nir = data[3].astype(np.float32)
ndvi = (nir - red) / (nir + red)
with rasterio.open(ndvi_file, 'w', **profile) as dst:
dst.write(ndvi, 1)
This block of code calculates the NDVI
(Normalized Difference Vegetation Index) using the red
and nir
bands from the input raster data. The red band is at index 2 of the data array, and the nir band is at index 3. The astype
method is used to convert the data type of the red
and nir
arrays to np.float32
, which is necessary for performing mathematical operations on them. The NDVI
formula is then applied to the red
and nir
arrays, and the resulting NDVI
values are written to a new raster file using rasterio
.
green = data[1].astype(np.float32)
ndwi = (green - nir) / (green + nir)
with rasterio.open(ndwi_file, 'w', **profile) as dst:
dst.write(ndwi, 1)
The first line loads the green
band data from the Sentinel-2
image into a numpy
array and converts the data type to 32-bit float by using np.float32
.
The second line calculates the NDWI
by subtracting the near-infrared (NIR)
band data from the green
band data and then dividing the result by the sum of the green band data and the NIR band data. The third line creates a new GeoTIFF file using rasterio.open()
and writes the NDWI data to it with dst.write()
.
swir = data[4].astype(np.float32)
ndmi = (nir - swir) / (nir + swir)
with rasterio.open(ndmi_file, 'w', **profile) as dst:
dst.write(ndmi, 1)
Same as the previous lines of code. Using SWIR
band for calculate NDMI
.
blue = data[0].astype(np.float32)
evi = 2.5 * ((nir - red) / (nir + 6 * red - 7.5 * blue + 1))
with rasterio.open(evi_file, 'w', **profile) as dst:
dst.write(evi, 1)
the first step is to extract the Blue
band from the input data and convert it to float32 data type. Next, the EVI
formula is applied using the extracted bands (Red, Near Infrared, and Blue) to compute a vegetation index that is sensitive to changes in vegetation density and structure. This formula multiplies the difference between the NIR
and Red
bands by a constant (2.5)
and then divides the result by a combination of the NIR
, Red
, and Blue
bands. The combination of bands in the denominator is designed to minimize the effects of atmospheric interference and improve the sensitivity of the EVI
to vegetation density. Finally, the resulting EVI
raster is written to a file using the Rasterio
library. The output file is specified by the variable evi_file
, and the file format and metadata are specified by the profile variable.
L = 0.5 # the soil adjustment factor
savi = ((nir - red) / (nir + red + L)) * (1 + L)
with rasterio.open(savi_file, 'w', **profile) as dst:
dst.write(savi, 1)
The first step is to set the soil adjustment factor L
to a value of 0.5
. Next, the SAVI
formula is applied using the extracted bands (Red and Near Infrared) to compute a vegetation index that is sensitive to changes in vegetation density while correcting for the effects of soil background. This formula multiplies the difference between the NIR
and Red
bands by a factor that corrects for soil background (1 + L)
and then divides the result by the sum of the NIR
and Red
bands plus the soil adjustment factor (L)
. The soil adjustment factor helps to account for the effect of soil reflectance on the vegetation index. Finally, the resulting SAVI
raster is written to a file using the Rasterio
library. The output file is specified by the variable savi_file
, and the file format and metadata are specified by the profile
variable.
Here is the full code:
import numpy as np
import rasterio
# Define the input file path
input_file = r'E:/Python_Tutorial/Indices/Sentinel-2_Bardhaman.tif'
# Define the output file paths
ndvi_file = r'E:/Python_Tutorial/Indices/ndvi.tif'
ndwi_file = r'E:/Python_Tutorial/Indices/ndwi.tif'
ndmi_file = r'E:/Python_Tutorial/Indices/ndmi.tif'
evi_file = r'E:/Python_Tutorial/Indices/evi.tif'
savi_file =r'E:/Python_Tutorial/Indices/savi.tif'
# Load the data into a numpy array
with rasterio.open(input_file) as src:
data = src.read()
profile = src.profile
profile.update(count=1)
# Calculate NDVI
red = data[2].astype(np.float32)
nir = data[3].astype(np.float32)
ndvi = (nir - red) / (nir + red)
with rasterio.open(ndvi_file, 'w', **profile) as dst:
dst.write(ndvi, 1)
# Calculate NDWI
green = data[1].astype(np.float32)
ndwi = (green - nir) / (green + nir)
with rasterio.open(ndwi_file, 'w', **profile) as dst:
dst.write(ndwi, 1)
# Calculate NDMI
swir = data[4].astype(np.float32)
ndmi = (nir - swir) / (nir + swir)
with rasterio.open(ndmi_file, 'w', **profile) as dst:
dst.write(ndmi, 1)
# Calculate EVI
blue = data[0].astype(np.float32)
evi = 2.5 * ((nir - red) / (nir + 6 * red - 7.5 * blue + 1))
with rasterio.open(evi_file, 'w', **profile) as dst:
dst.write(evi, 1)
# Calculate SAVI
L = 0.5 # the soil adjustment factor
savi = ((nir - red) / (nir + red + L)) * (1 + L)
with rasterio.open(savi_file, 'w', **profile) as dst:
dst.write(savi, 1)
Calculating multiple vegetation indices (such as NDVI, NDWI, NDMI, EVI, SAVI) in a single time using Python can provide several benefits for remote sensing analysis:
-
Efficiency: By processing multiple vegetation indices in a single script, you can save time and processing power by avoiding redundant calculations and loading data into memory multiple times.
-
Comparative analysis: Having multiple vegetation indices calculated in one script allows for a more comprehensive analysis of vegetation conditions, which can help in comparative studies and change detection analysis.
-
Diverse data sources: Different vegetation indices respond to different vegetation characteristics, allowing for a more diverse set of data to be used in remote sensing analysis. With Python, you can calculate all of these indices from a single dataset, enabling a more holistic analysis.
-
Visualization: Python has many tools for data visualization, and having multiple vegetation indices calculated in a single script allows for easy visualization and comparison of the different indices.
-
Flexibility: Python is a flexible programming language that can be used to automate and customize workflows for remote sensing analysis. By processing multiple vegetation indices in one script, you can easily modify the workflow to suit your research question and the properties of the input data.
In summary, calculating multiple vegetation indices in a single Python script can help to save time, provide a more comprehensive analysis, enable the use of diverse data sources, facilitate visualization, and offer flexibility in workflow customization
-
Email: koushikghosh1a2s@gmail.com.
-
LinkedIn: Koushik Ghosh
-
ResearchGate: Koushik Ghosh