Creating Remote Sensing Products with Field Spectroscopy Data

Introduction to Tutorial

We will continue to work with the Absolute_Means_smoothed dataframe generated from the last session to produce remote sensing products like vegetation indices and convolved data.

Tutorial Steps

Part One – Vegetation Indices

With your data now processed and converted into absolute spectral reflectance, we can look more closely at differences between the feature types. We can use dimensionality reduction methods, such as spectral indices, to do this.

A number of spectral indices have been designed to highlight different vegetation properties. In this section, we will use some of them to explore differences between our feature types. A full description of the indices used can be found here

One of the most commonly used vegetation indices is the Normalized Difference Vegetation Index (NDVI). When using Sentinel-2 data, the equation to calculate NDVI is –

$$ NDVI = \dfrac{B_8 - B_4}{B_8 + B_4} $$

where Band 8 is equivalent to 842 nm, and Band 4 is equivalent to 665 nm.

1. Using the inbuilt Python function .iloc, we can calculate NDVI for our leaf means. First, let’s look at what iloc does.

   Absolute_Means_smoothed.iloc[315]
   

2. Here, iloc takes the index number for the Absolute_Leaf_Means dataframe, and reports what is at that value – here, the reflectance of our different leaf groups, and also, as ‘Name’, the wavelength at which those reflectances were measured. You should see that ‘Name’ equals 842.0 nm. Notice the offset between the iloc value given and our wavelength value – 350, which is the wavelength value at which our dataframe begins. With this in mind, we can now proceed to calculate NDVI, first generating a new dataframe to store our output:

   df = []
   Indices = pd.DataFrame(df)
   Indices['NDVI'] = ((Absolute_Means_smoothed.iloc[492] - Absolute_Means_smoothed.iloc[315]) / 
                                  (Absolute_Means_smoothed.iloc[492] + Absolute_Means_smoothed.iloc[315])) 
   print(Indices['NDVI'])
   

   • Question: - The Carotenoid Reflectance Index looks at the the concentrations of carotenoids in vegetation. Carotenoids function in light absorption processes in plants, as well as in protecting plants from the harmful effects of too much light. Higher CRI1 values mean greater carotenoid concentration relative to chlorophyll. The value of this index ranges from 0 to more than 15. The common range for green vegetation is 1 to 12: $$CRI = \dfrac{1}{510} - \dfrac{1}{550}$$ With this information, calculate the CRI for your data set.
     Answer

           
           Indices['CRI'] = ((1 / Absolute_Means_smoothed.iloc[160]) - (1 / Absolute_Means_smoothed.iloc[200]))
           
           

   • Question: - The Plant Senescence Reflectance Index maximizes the sensitivity of the index to the ratio of bulk carotenoids (for example, alpha-carotene and beta-carotene) to chlorophyll. An increase in PSRI indicates increased canopy stress (carotenoid pigment), the onset of canopy senescence, and plant fruit ripening. Applications include vegetation health monitoring, plant physiological stress detection, and crop production and yield analysis. The value of this index ranges from -1 to 1. The common range for green vegetation is -0.1 to 0.2: $$PSRI = \dfrac{680 - 500}{750}$$ With this information, calculate the PSRI for your data set.
     Answer

           
           Indices['PSRI'] = ((Absolute_Means_smoothed.iloc[330] - Absolute_Means_smoothed.iloc[150] ) / ( Absolute_Means_smoothed.iloc[400]))
           
           

   • Question: - The Cellulose Absorption Index is a reflectance measurement that quantifies dried plant material. Strong absorption features present in the 2000nm to 2200nm range indicate strong presence of cellulose. $$CAI = 0.5(2000 + 2200) - 2100$$ Determine the Cellulose Absorption Index for your data set.
     Answer

           
           Indices['CAI'] = (0.5 * (Absolute_Means_smoothed.iloc[1650] + Absolute_Means_smoothed.iloc[1850])) - (Absolute_Means_smoothed.iloc[1750])
           
           

Part Two – Convolution

If your research relates to specific multispectral imaging sensors e.g. Sentinel 2, it can be useful to resample your hyperspectral data so that it matches the bands of your specific sensor. We can do this by convolving the hyperspectral data to a multispectral sensor’s “spectral response function”.

The Field Spectroscopy Facility has a Python package available for swift convolution of hyperspectral data to Sentinel-2 and WorldView multispectral imagery. For this tutorial, we’ll look at how this package works, by convolving the data for one band of Sentinel-2 – Band 4.

Band 4 of the Sentinel-2 sensor is also know as the “Red” band. The centre wavelength – the wavelength at which it has the most sensitivity – is 665 nm.

The band, however, is sensitive to other wavelengths (what is termed as its ‘spectral response feature’), from 646 nm to 684 nm.

In order to accurately compare our hyperspectral data to Band 4, we must convolve the hyperspectral data within the wavelength regions of Band 4’s SRF. This requires knowledge of the varying sensitivity of the band. The .csv file ‘Sentinel 2 SRF.csv’, contained in the folder which you downloaded, contains this information.

1. Read in the Sentinel-2 SRF as a pandas dataframe:

   Sentinel_2_bands = pd.read_csv('Sentinel 2 SRF.csv', index_col=0, header=0)
   

2. Knowing the starting and ending wavelength of the SRF of interest, we can now conduct convolution, which mathematically can be expressed as the result of the trapezoidal integration of the product of the multiplication hyperspectral data over the SRF region by the SRF, divided by the result of the trapezoidal integration of the SRF itself: $$R_v^k=\int_{v_1^k}^{v_2^k} R_h(v) S_b(v) \mathrm{d} v / \int_{v_1^k}^{v_2^k} S_b(v) \mathrm{d} v$$

   Indices['Band 4 - Red'] = (np.trapezoid((Absolute_Means_smoothed.iloc[298:334]).mul(Sentinel_2_bands.iloc[298:334, 3],axis = 0), axis =0)) / (np.trapezoid((Sentinel_2_bands.iloc[298:334, 3]), axis = 0))
   Indices['Band 4 - Red']
   

3. Considering this code in more detail – np.trapezoid is the mathematical function conducted on the data; .iloc[298:334] gives the region to conduct the convolution (iloc[298] is equal to wavelength 646 nm (the narrowest edge of the SRF), and iloc[334] is equal to wavelength 684 nm, the broadest range of the SRF); and the ,3 after 298:334 for the integration of the Sentinel 2 bands relates to the column number in which the SRF for Band 4 is stored (the relation between column number and band in the SRF file is Band Number - 1). The final output gives the reflectance value of the convoluted hyperspectral data.

   • Question: - Repeat the process as outlined above, but for Band 3, ‘Green’. Band 3 has a centre wavelength at 560 nm. The SRF for Band 3 extends from 538 nm to 583 nm. Modifying the code above, create a row in Indices called ‘Band 3 - Green’ that provides the output of the convolution of the hyperspectral data to the Band 3 SRF.
     Answer

           
           Indices['Band 3 - Green'] = (np.trapezoid((Absolute_Means_smoothed.iloc[188:233]).mul(Sentinel_2_bands.iloc[188:233, 2],axis = 0), axis =0)) / (np.trapezoid((Sentinel_2_bands.iloc[188:233, 2]), axis = 0))
           
           

4. It is possible to use the same code to generate the values for each of the Sentinel-2 bands, or indeed for any satellite sensor provided you know the spectral response feature of it. To expedite the process, however, we will use the FieldSpecUtils package, which we installed earlier:

   from FieldSpecUtils import Convolution
   

5. We can now run our convolution as follows:

   Convolution.S2(Absolute_Means_smoothed, Sentinel_2_bands)
   

6. Our hyperspectral data has now been convolved to multispectral bands equivalent to Sentinel-2. We can see how this looks by plotting the data output of the convolution, Ploths_with_convolved_bands.csv:

   collated_convolved = pd.read_csv("Plots_with_convolved_bands.csv", index_col = 0)
   collated_convolved.plot(figsize=(9, 6), legend = True, ylim=(0,1), linestyle='--', marker='o')
   xlabel("Centre Wavelength (nm) of Band with Band Number")
   ylabel("Absolute Reflectance")
   title("Hyperspectral convolved broadband reflectance")
   plt.xticks([0, 1, 2, 3, 4, 5, 6, 7, 8, 9],
           ["490 -- Band 2", "560 -- Band 3", "665 -- Band 4", "705 -- Band 5", "740 -- Band 6",
            "783 -- Band 7", "842 -- Band 8", "865 -- Band 8a", "1610 -- Band 11", "2190 -- Band 12"],
           rotation=20)
   legend(title = "Plot number and vegetation type")
   plt.tight_layout()
   plt.savefig("Convolved_Hyperspectral_Plot_Data_to_S2_Multispectral_Bands.pdf")
   plt.show()
   

Wrapup

In this session, you’ve learned how to utilize field spectrometer data to generate remote sensing products like indices or convolved products. More advanced techniques, such as principal component analysis, can also be conducted using data in this format.

In the next session, we will discuss the process of UAV imagery using ESA SNAP, but for the spectral unmixing session in the final tutorial, you will require your spectral library as generated in the last session.