Standard Environmental Drone Monitoring Techniques: NDVI, ENDVI and SAVI

With a standard filter we have demonstrated the use of a drone under single deployment to create Near-Infrared images.

The next step is to create indices which are important for the specific application in environmental monitoring. Of particular interest is the monitoring of regions which are vulnerable to natural and human induced disasters such as forest fires

We are also interested in comparing regions based on the amount of soil exposed due to factors such as grazing, human activity, etc.

To do this we must consider that there are many environmental indices to choose from in the field of NIR monitoring (and even more with other spectral bands).

Using experimental flights taken with an NIR camera using a bandpass filter with the following spectrum:


We have developed 3 important classification indices to consider using standard NIR filters on converted IR Cameras deployed on a drone used in environmental monitoring.

NDVI - Normalised Differential Vegetation Index

We have discussed NDVI before on this blog and will not go through it in detail again, however it is important to point out the differences between the satellite use of NDVI and the Drone/Aerial Monitoring implentation of the NDVI with regards to the different but essentially equivalent formula.

With satellite imagery we use an NDVI formula which has the best possible leveraging between a single band in the visible spectrum vs the NIR band. In an RGB image taken of the Earth's surface by a satellite, the blue portion of the visible spectrum (roughly 400nm-450nm) will diminished the most significantly by the atmosphere of the Earth, which scatters blue light more strongly than green or red. Hence, water aside*, features of the surface which reflect blue are relatively diminished in intensity. Hence, most satellite cameras used in Earth observation of vegetation ignore the blue region altogether and focus on Red, Green and several key IR bands.

Therefore, for satellite imagery, we use what has been come to be known as the standard NDVI formula:

With the drone imaging we perform, we can use a different but equivalent formula which gives a far better leveraging between NIR and the Blue channel for several reasons.

• There is less atmosphere between the surface feature being examine and the drone CCD sensor, so blue light scattering is not an option
• Some vegetation is in fact Red when Healthy, due to pigments. Blue however is not reflected by healthy plants, due to chlorophyll A and B both having prominent absorption in the blue region of the visible spectrum.
• New CCD sensors used in drones (such as the Hasselblad) are generally designed to have equal sensitivity between the Red and Blue bands 
• The Red band can be used as the NIR band in a converted NIR camera, making for easier data acquisition with existing cameras

So for Drone imaging we use the following NDVI formula, sometimes called the Infrablue NDVI formula.

Using our code in Matlab (see Appendix) from the following Near-Infrared DNG image

We process the following NDVI image

We can also use a similar code in Python that can do essentially the same processing, just with a different colormap scale:

At this point it all depends on personal preference which format to use, Matlab or Python. Python tends to be better for preserving bit-depth for use in further applications we will discuss in a future article on 3D reconstruction from NDVI processed images. I also prefer the colormaps that I have been able to create in Python myself as they provide a much sharper definition of the different NDVI scales in an image overall.

*(Water can reflect blue light in higher intensity because it is, in effect, reflecting the color of the sky (which is itself blue because of blue light being scattered in the first place) ).

ENDVI - Enhanced Normalised Differential Vegetation Index

Extended Normalized Difference Vegetation Index (ENDVI) data gathered for vegetative health
monitoring can be used to provide similar, but spectrally different information as compared
to traditional NDVI data.

Soil background, differing atmospheric conditions and various types of vegetation can all influence the reflection of visible light somewhat differently.

ENDVI uses Blue and Green visible light instead of the Red only method of the standard NDVI algorithm. This allows for better isolation of plant health indicators.

A normal healthy plant will reflect both visible Green and NIR light while absorbing Blue visible light, hence plants appear Green in the Visible Spectrum and will Glow Brightly in the Near-Infrared.

In the ENDVI measurement we are calculating the relationship between high absorption of this Blue light and high reflectance of Green and NIR radiation as an indicator of plant health.

Hence ENDVI analysis may, at times, be able to impart more accurate or reliable information regarding plant or crop health by additional reference of information in the green portion of the spectrum with respect to NIR.

It also has an added advantage in terms of distinguishing the differences between algae, lichens and plants on ground and, more importantly, in water which is critical for the monitoring of algal blooms and other contaminants which reflect in different intensity in the Near-Infrared than vascular plants,

As discussed in a past article, Vascular plants reflect NIR more strongly when healthy due to the presence of spongy mesophyll tissue which algae and lichens lack regardless of their individual health. Thus ENDVI gives more ability to classify plants and objects that appear bright in the NIR overall and with particular respect to their position in the ecosystem.

Again, using our code in Matlab on the same NIR DNG image we process the image from NIR to ENDVI

And we can do the same with generating ENDVI with our Python code:

SAVI - Soil-adjusted vegetation index

The SAVI is structured similar to the NDVI but with the addition of a soil reflectance correction factor called L that leverages the overal NIR brightness of soils so as to not wash out any potential features in the soil.

Just as in the case with NDVI, the satellite imaging version of the SAVI formula contains the Red band for use in leveraging against the NIR:

For our drone version, as before with NDVI, we use the blue channel as our leveraging against the NIR, which has been swapped with the camera's red channel using our InfraBlue Filter. The formula is hence amended to:

L is a constant and is related to the slope of the soil-line in a feature-space plot.

Hence the value of L varies by the amount or cover of green vegetation: in very high vegetation regions, i.e. thick forest, L=0; and in areas with no green vegetation, i.e. deserts, steep slopes, canyons, deserts, beaches etc then L=1. Generally, an L=0.5 works well in most situations (i.e. mixed vegetation cover, the red region in the NIR vs VIS(Blue) Plot).

So L=0.5 (half) is usually the default value used.
Note: When L=0, then SAVI = NDVI.


Again, using our code in Matlab on the same NIR DNG image we process the image from NIR to SAVI

And we can do the same with generating SAVI with our Python code:

SAVI can be very useful in observing very dry or grazed regions where soil and rock is exposed either by climate, animal or human activity.

SAVI has particular interest not only in helping to distinguish low level vegetation from soil and bare soil from rock but it has particular attractiveness in archaeology in distinguishing markings and artificial arrangements in the soil and rock surfaces in a landscape that would otherwise be depicted as largely homogeneous in colour and thus featureless.

This is particularly helpful when mapping regions which may have been used for farming or habitation long ago but have been largely abandoned with rock walls and areas of human habitation receding into the landscape.

In the area of archaelogy it is of course very important to map regions where ancient settlements once stood, and to discover new artifacts.

Also in the field of rewilding, old stone walls that divided farmlands in the past can be an obstacle to migrating wildlife and are an important feature with regards to distribution of species and bottlenecks involving their range in an area and may shape their behavior in feeding as well as predator-prey relations and human interactions.

Matlab Scripts:

All Matlab Files that load DNG need the following loadDNG function file:

function [rawData, tinfo]= loadDNG(dngFilename)   

    if(exist(dngFilename,'file'))

        tinfo = imfinfo(dngFilename);

        t = Tiff(dngFilename,'r');

        offsets = getTag(t,'SubIFD');

        setSubDirectory(t,offsets(1));

        rawData = t.read();

        t.close();

    else

        if(nargin<1 || isempty(dngFilename))

            dngFilename = 'File';

        end

        fprintf(1,'%s could not be found\n',dngFilename);

        rawData = [];

    end

end

Full Poster Script:

%Experimental Vegetation Index Mapping program using DJI Mavic Pro DNG 16-bit images taken using InfraBlue Filter 

%Works with DNG Images only - needs loadDNG script in same directory

%(c)-J. Campbell MuonRay Enterprises Drone-based Vegetation Index Project 2017-2019 

rawData = loadDNG('DJI_0592.DNG'); % load it "functionally" from the command line

%DNG is just a fancy form of TIFF. So we can read DNG with TIFF class in MATLAB

warning off MATLAB:tifflib:TIFFReadDirectory:libraryWarning

%Transpose Matrices for image data because DNG image format is row-major, and Matlab matrix format is column-major.

%rawData = permute(rawData, [2 1 3]);

%Assume Bayer mosaic sensor alignment.

%Seperate to mosaic components.

J1 = rawData(1:2:end, 1:2:end);

J2 = rawData(1:2:end, 2:2:end);

J3 = rawData(2:2:end, 1:2:end);

J4 = rawData(2:2:end, 2:2:end);

figure(1);

J1 = imadjust(J1, [0.09, 0.91], [], 0.45); %Adjust image intensity, and use gamma value of 0.45

J1(:,:,1) = J1(:,:,1)*0.80; %Reduce the overall colour temperature by factor of 0.80

B = im2single(J4(:,:,1)); %VIS Channel Blue - which is used as a magnitude scale in NDVI

R = im2single(J1(:,:,1)); %NIR Chaneel RED which is the IR reflectance of healthy vegeation 

G = im2single(J2(:,:,1)); %VIS Channel Green - which is used as a magnitude scale in ENDVI

L = 0.7; % Soil Reflectance Correction Factor for SAVI Index

%InfraBlue NDVI

InfraBlue = (R - B)./(R + B);

InfraBlue = double(InfraBlue);

%% Stretch NDVI to 0-255 and convert to 8-bit unsigned integer

InfraBlue = floor((InfraBlue + 1) * 32767); % [-1 1] -> [0 256] for 8-bit display range(*128), [0, 65535] for 16-bit display range (*32500)

InfraBlue(InfraBlue < 0) = 0;             % may not really be necessary, just in case & for symmetry

InfraBlue(InfraBlue > 65535) = 65535;         % in case the original value was exactly 1

InfraBlue = uint16(round(InfraBlue));             % change data type from double to uint16

% InfraBlue = uint8(InfraBlue);             % change data type from double to uint8

NDVImag = double(InfraBlue);

%ENDVI - Green Leveraged

ENDVI = ((R+G) - (2*B))./((R+G) + (2*B));

ENDVI = double(ENDVI);

%% Stretch ENDVI to 0-255 and convert to 8-bit unsigned integer

ENDVI = floor((ENDVI + 1) * 32767); % [-1 1] -> [0 256] for 8-bit display range(*128), [0, 65535] for 16-bit display range (*32500)

ENDVI(ENDVI < 0) = 0;             % may not really be necessary, just in case & for symmetry

ENDVI(ENDVI > 65535) = 65535;         % in case the original value was exactly 1

ENDVI = uint16(round(ENDVI));             % change data type from double to uint16

% InfraBlue = uint8(InfraBlue);             % change data type from double to uint8

ENDVImag = double(ENDVI);

%SAVI - Soil Reflectance Leveraged 

%The SAVI is structured similar to the NDVI but with the addition of a

%“soil reflectance correction factor” L.

%L is a constant (related to the slope of the soil-line in a feature-space plot)

%Hence the value of L varies by the amount or cover of green vegetation: in very high vegetation regions, 

%L=0; and in areas with no green vegetation, L=1. Generally, an L=0.5 works

%well in most situations (i.e. mixed vegetation cover)

%So 0.5 (half) is the default value used. When L=0, then SAVI = NDVI.

SAVI = (((R-B)*(1+L))./(R+B+L));

SAVI = double(SAVI);

%% Stretch ENDVI to 0-255 and convert to 8-bit unsigned integer

SAVI = floor((SAVI + 1) * 32767); % [-1 1] -> [0 256] for 8-bit display range(*128), [0, 65535] for 16-bit display range (*32500)

SAVI(SAVI < 0) = 0;             % may not really be necessary, just in case & for symmetry

SAVI(SAVI > 65535) = 65535;         % in case the original value was exactly 1

SAVI = uint16(round(SAVI));             % change data type from double to uint16

% InfraBlue = uint8(InfraBlue);             % change data type from double to uint8

SAVImag = double(SAVI);

% Display them all.

subplot(3, 3, 2);

imshow(rawData);

fontSize = 20;

title('Captured Image', 'FontSize', fontSize)

subplot(3, 3, 4);

imshow(R);

title('Red Channel', 'FontSize', fontSize)

subplot(3, 3, 5);

imshow(G)

title('Green Channel', 'FontSize', fontSize)

subplot(3, 3, 6);

imshow(B);

title('Blue Channel', 'FontSize', fontSize)

subplot(3, 3, 7);

myColorMap = jet(65535); % Whatever you want.

rgbImage = ind2rgb(NDVImag, myColorMap);

imagesc(rgbImage,[0 1])

%imshow(recombinedRGBImage);

title('NDVI Image', 'FontSize', fontSize)

subplot(3, 3, 8);

myColorMap = jet(65535); % Whatever you want.

rgbImage2 = ind2rgb(ENDVImag, myColorMap);

imagesc(rgbImage2,[0 1])

%imshow(recombinedRGBImage);

title('ENDVI Image', 'FontSize', fontSize)

subplot(3, 3, 9);

myColorMap = jet(65535); % Whatever you want.

rgbImage3 = ind2rgb(SAVImag, myColorMap);

imagesc(rgbImage3,[0 1])

%imshow(recombinedRGBImage);

title('SAVI Image', 'FontSize', fontSize)

% Set up figure properties:

% Enlarge figure to full screen.

set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0, 1, 1]);

% Get rid of tool bar and pulldown menus that are along top of figure.

% set(gcf, 'Toolbar', 'none', 'Menu', 'none');

% Give a name to the title bar.

set(gcf, 'Name', 'Demo NDVI Poster', 'NumberTitle', 'Off')

%Save NDVI, ENDVI and SAVI images to PNG format file

imwrite(rgbImage, 'NDVI.png');

imwrite(rgbImage2, 'ENDVI.png');

imwrite(rgbImage3, 'SAVI.png');

Python Script:

"""
Experimental Vegetation Index Mapping program using DJI Mavic 2 Pro
JPEG 16-bit combo images taken using InfraBlue Filter

This script uses rawpy

rawpy is an easy-to-use Python wrapper for the LibRaw library.
rawpy works natively with numpy arrays and supports a lot of options,
including direct access to the unprocessed Bayer data
It also contains some extra functionality for finding and repairing hot/dead pixels.
import rawpy.enhance for this

First, install the LibRaw library on your system.

pip install libraw.py

then install rawpy

pip install rawpy

%(c)-J. Campbell MuonRay Enterprises 2019
This Python script was created using the Spyder Editor

"""

import warnings
warnings.filterwarnings('ignore')
import numpy as np
from matplotlib import pyplot as plt  # For image viewing
import rawpy

#!/usr/bin/python

from matplotlib import colors
from matplotlib import ticker
from matplotlib.colors import LinearSegmentedColormap

#dng reading requires libraw to work


#a nice selection of grayscale colour palettes
cols1 = ['blue', 'green', 'yellow', 'red']
cols2 =  ['gray', 'gray', 'red', 'yellow', 'green']
cols3 = ['gray', 'blue', 'green', 'yellow', 'red']

cols4 = ['black', 'gray', 'blue', 'green', 'yellow', 'red']

def create_colormap(args):
    return LinearSegmentedColormap.from_list(name='custom1', colors=cols3)

#colour bar to match grayscale units
def create_colorbar(fig, image):
        position = fig.add_axes([0.125, 0.19, 0.2, 0.05])
        norm = colors.Normalize(vmin=-1., vmax=1.)
        cbar = plt.colorbar(image,
                            cax=position,
                            orientation='horizontal',
                            norm=norm)
        cbar.ax.tick_params(labelsize=6)
        tick_locator = ticker.MaxNLocator(nbins=3)
        cbar.locator = tick_locator
        cbar.update_ticks()
        cbar.set_label("NDVI", fontsize=10, x=0.5, y=0.5, labelpad=-25)


# Open an image
infile = 'DJI_0592.DNG'

raw = rawpy.imread(infile)

rgb = raw.postprocess(gamma=(1,1), no_auto_bright=True, output_bps=16)

R = rgb[:,:,0]
G = rgb[:,:,1]
B = rgb[:,:,2]
# Get the red band from the rgb image, and open it as a numpy matrix
#NIR = image[:, :, 0]
         
ir = np.asarray(R, float)



# Get one of the IR image bands (all bands should be same)
#blue = image[:, :, 2]

r = np.asarray(B, float)



# Create a numpy matrix of zeros to hold the calculated NDVI values for each pixel
ndvi = np.zeros(r.size)  # The NDVI image will be the same size as the input image

# Calculate NDVI
ndvi = np.true_divide(np.subtract(ir, r), np.add(ir, r))

'''
Extended Normalized Difference Vegetation Index (ENDVI) data gathered for vegetative health
monitoring can be used to provide similar, but spectrally different information as compared
to traditional NDVI data.

Soil background, differing atmospheric conditions and various types
of vegetation can all influence the reflection of visible light somewhat differently.

ENDVI uses Blue and Green visible light instead of the Red only method of the standard NDVI
algorithm. This allows for better isolation of plant health indicators.

A normal healthy plant will reflect both visible Green and NIR light while
absorbing Blue visible light, hence plants appear Green in the Visible Spectrum and
Will Glow Brightly in the Near-Infrared.

In the ENDVI measurement we are calculating the relationship between
high absorption of this Blue light and high reflectance of Green and NIR waves
as an indicator of plant health.

Hence ENDVI analysis may, at times, be able to impart more accurate or reliable information
regarding plant or crop health by additional reference of information in
the green portion of the spectrum with respect to NIR.

'''
#ENDVI = [(NIR + Green) – (2 * Blue)] / [(NIR + Green) + (2 * Blue)]
# Get one of the green image bands (all bands should be same)


g = np.asarray(G, float)

#(NIR + Green)
irg = np.add(ir, g)

# Calculate ENDVI
#ndvi = np.true_divide(np.subtract(irg, 2*r), np.add(irg, 2*r))

#%SAVI - Soil Reflectance Leveraged

#The SAVI is structured similar to the NDVI but with the addition of a
#“soil reflectance correction factor” L.
#L is a constant (related to the slope of the soil-line in a feature-space plot)
#Hence the value of L varies by the amount or cover of green vegetation: in very high vegetation regions,
#L=0; and in areas with no green vegetation, L=1. Generally, an L=0.5 works
#well in most situations (i.e. mixed vegetation cover)
#So 0.5 (half) is the default value used. When L=0, then SAVI = NDVI.



#SAVI = [((R-B)*(1+L))./(R+B+L)];

L=0.5;
one = 1;

#rplusb = np.add(ir, r)
#rminusb = np.subtract(ir, r)
#oneplusL = np.add(one, L)

#ndvi = np.true_divide(np.multiply(rminusb, oneplusL), np.add(rplusb, L))


# Display the results
output_name = 'InfraBlueNDVI3.jpg'

fig, ax = plt.subplots()
image = ax.imshow(ndvi, cmap=create_colormap(colors))
plt.axis('off')
     
create_colorbar(fig, image)
extent = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
     
fig.savefig(output_name, dpi=600, transparent=True, bbox_inches=extent, pad_inches=0)
        # plt.show()