Today I wanted to show a simple example of how to calculate NDVI from Sentinel-2A data. This is the first time that I have used Sentinel-2A data, and I was interested in the higher (10 meter) resolution data for looking at crop conditions. The use case for this tutorial is an intensely cultivated region of Turner County, South Dakota, an area known for corn and soy production. I chose this region as a tribute to my Grandmother, Lillian, who passed away this past Friday, October 18th, 2017. My grandfather and grandmother lived and farmed in this area for many, many years. Below is a photo I took this past week while out visiting my great grandfathers farm.

As you can see in the picture above, the grass is still green and many of the trees still have leaves–we have had an unseasonably warm autumn thus far in southeast SD.

Now, back to the satellite imagery. As I mentioned, this was my first experience working with Sentinel-2A data, and also my first time working with JPEG 2000 (.jp2) format. Unfortunately, the version of GDAL that I had installed on my system did not recognize the JPEG 2000 driver, and I had to set up an alternative environment with a more recent version of GDAL in order to process these data. Here is an example of the environment I created specifically for this tutorial (.yml file):

name: cole6

channels:

– conda-forge

– defaults

dependencies:

– python=2.7.12

– gdal=2.1.*

– rasterio

– spyder

If you plan to use the code from this tutorial and have never worked with these data before, I recommend setting up this or a similar environment. If you already use Conda, you can activate set up a new environment with the packages installed as above by typing:

$ conda env create -f cole6.yml

However, I was having trouble installing Matplotlib while keeping this version of GDAL due to dependency conflicts, so for this tutorial I simply did the processing in Python and then used QGIS in order to generate the ‘pretty pictures’ seen below.

To find the Sentinel-2A data used in this tutorial, I went to EarthExplorer. I browsed the area over my grandparents farm until I found a cloud-free image, the most recent of which I found was from September 9th, 2017. I then downloaded and unzipped the file, and used bands 4 (red) and 8 (NIR) to generate the Normalized Difference Vegetation Index, or NDVI. NDVI is a measure of vegetation greenness that can give some indication of the health of vegetation.

The short python script below brings in the red and NIR bands, calculates NDVI, and creates an output GeoTIFF file of the NDVI for this particular Sentinel observation. Here we can see the finished product:

Areas in shades of brown indicate low NDVI, and areas in increasingly dark shades of green indicate higher NDVI. Here the city of Sioux Falls, and in particular the dense commercial districts and downtown areas in the upper right hand corner appear dark brown due to the high concentrations of impervious surfaces, such as buildings and roads. There are a few small lakes in the image as well that appear brown. The riparian areas in the lower left hand corner along the James River Valley that appear in shades of tan are likely grass and rangelands that have mostly dried out by mid-September. However, the rectangular blocks in varying shades of green show the most prevalent, and economically important, land cover of the region: Croplands.

Below, I have zoomed in on my great grandfather and grandfathers farms:

I was amazed by the level of intra-field variability in NDVI that is captured by the 10 meter resolution of the Sentinel-2A data. This information can be useful to farmers as areas within a field that are in lighter shades of tan and green may be suffering from a water or nutrient deficiency.

Link to the source code on git: website/Sentinel_NDVI.py

Below is the code used to generate the images above:

#!/usr/bin/env python2 # -*- coding: utf-8 -*- """ ############################################################################### # How to: Calculate NDVI from Sentinel-2A # ############################################################################### # @author: Cole Krehbiel # # Last Updated: 10-18-17 # ############################################################################### """ # Import libraries import glob import numpy as np from osgeo import gdal # If GDAL doesn't recognize jp2 format, check version</pre> # Set input directory in_dir = '/Users/cole/Desktop/Sentinel/' # Search directory for desired bands red_file = glob.glob(in_dir + '**B04.jp2') # red band nir_file = glob.glob(in_dir + '**B08.jp2') # nir band # Define a function to calculate NDVI using band arrays for red, NIR bands def ndvi(red, nir): return ((nir - red)/(nir + red)) # Open each band using gdal red_link = gdal.Open(red_file[0]) nir_link = gdal.Open(nir_file[0]) # read in each band as array and convert to float for calculations red = red_link.ReadAsArray().astype(np.float) nir = nir_link.ReadAsArray().astype(np.float) # Call the ndvi() function on red, NIR bands ndvi2 = ndvi(red, nir) # Create output filename based on input name outfile_name = red_file[0].split('_B')[0] + '_NDVI.tif' x_pixels = ndvi2.shape[0] # number of pixels in x y_pixels = ndvi2.shape[1] # number of pixels in y # Set up output GeoTIFF driver = gdal.GetDriverByName('GTiff') # Create driver using output filename, x and y pixels, # of bands, and datatype ndvi_data = driver.Create(outfile_name,x_pixels, y_pixels, 1,gdal.GDT_Float32) # Set NDVI array as the 1 output raster band ndvi_data.GetRasterBand(1).WriteArray(ndvi2) # Setting up the coordinate reference system of the output GeoTIFF geotrans=red_link.GetGeoTransform() # Grab input GeoTranform information proj=red_link.GetProjection() # Grab projection information from input file # now set GeoTransform parameters and projection on the output file ndvi_data.SetGeoTransform(geotrans) ndvi_data.SetProjection(proj) ndvi_data.FlushCache() ndvi_data=None ############################################################################### <pre>