<jupyter_start><jupyter_text>Session 4 Solution - Calculate NDVI Part 2 Exercise  Open and run notebook<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline

import sys
import datacube

sys.path.append("../Scripts")
from deafrica_datahandling import load_ard
from deafrica_plotting import rgb
from deafrica_plotting import display_map
from odc.algo import xr_geomedian
dc = datacube.Datacube(app="Calculate_ndvi")<jupyter_output><empty_output><jupyter_text>In the following cell, we set x and y equal to "tuples" (pairs) indicating the longitude and latitude extents of our region of interest, and then display a map showing the region.<jupyter_code>x=(-6.1495, -6.1380)
y=(13.9182, 13.9111)

display_map(x, y)<jupyter_output><empty_output><jupyter_text>Making changes to the load cell Adding nir measurement Here we load data via the ```load_ard``` function.  We specify the instance of the datacube (```dc```), the product that we wish to load (```s2_l2a``` for Sentinel-2), the extents of our query (in x, y, and time), the list of bands to fetch (```['red' ,'green', 'blue', 'nir']```), the spatial resolution of our data (20 meters), and we group nearby results by solar day.<jupyter_code>sentinel_2_ds = load_ard(
    dc=dc,
    products=["s2_l2a"],
    x=x, y=y,
    time=("2019-01", "2019-12"),
    output_crs="EPSG:6933",
    measurements=['red', 'green', 'blue', 'nir'],
    resolution=(-10, 10),
    group_by='solar_day')<jupyter_output>Using pixel quality parameters for Sentinel 2
Finding datasets
    s2_l2a
Applying pixel quality/cloud mask
Loading 71 time steps<jupyter_text>Plot timesteps Here we create a list of time slices to plot, and store it as the variable ```timesteps```.  The values in the list are indices into the list of satellite acquisitions, starting at zero.  The values here, ```[1, 6, 8]```, therefore refer to the second, seventh, and ninth acquisitions.  There should be 71 acquisitions for this particular dataset/extent combination, so any values from 0 to 70 can be used.We then call the ```rgb``` function to produce a series of plots from the data that we loaded, using the bands that we specify (```['red', 'green', 'blue']```), for the time indices stored in ```timesteps```.<jupyter_code>timesteps = timesteps = [1, 6, 8]

rgb(sentinel_2_ds, bands=['red', 'green', 'blue'], index=timesteps, size=5)<jupyter_output><empty_output><jupyter_text>Resampling the dataset Here we create a new variable (```resample_sentinel_2_ds```) which describes a grouping of the data we previously loaded (```sentinel_2_ds```), divided into 3-month (quarterly) intervals.<jupyter_code>resample_sentinel_2_ds = sentinel_2_ds.resample(time='3MS')<jupyter_output><empty_output><jupyter_text>Compute the geomedian Below, we take the data which we just split into quarterly intervals, and we apply the ```xr_geomedian``` function to each interval.  We store the result into the ```geomedian_resample``` variable.<jupyter_code>geomedian_resample = resample_sentinel_2_ds.map(xr_geomedian)
rgb(geomedian_resample, bands=['red', 'green', 'blue'], col="time", col_wrap=4, size=5)<jupyter_output><empty_output><jupyter_text>Comparing the two datasets In the next two cells, we simply print our datasets to screen.  The first cell contains the data that we originally loaded from the Open Data Cube, and the second cell contains the data after resampling it into quarterly geomedian composites.  Note that the first dataset contains many more time slices than the second.<jupyter_code>sentinel_2_ds
geomedian_resample<jupyter_output><empty_output><jupyter_text>Manually calculate an index Calculate NDVI In the following cell, we compute NDVI values from our resampled data using the formula: $$ NDVI = \frac{(NIR - Red) }{ (NIR + Red)} $$To keep our code tidy, we create variables ```nir``` and ```red``` and set them to equal ```geomedian_resample.nir``` and ```geomedian_resample.red``` respectively.<jupyter_code>nir = geomedian_resample.nir
red = geomedian_resample.red

NDVI = (nir - red) / (nir + red)<jupyter_output><empty_output><jupyter_text>Plot NDVI for each geomedian Here we call the ```plot``` function on our NDVI data, which gives us an image for each quarterly interval.<jupyter_code>NDVI.plot(col='time', vmin=-0.50, vmax=0.8, cmap='RdYlGn')<jupyter_output><empty_output><jupyter_text>Plot time series of the NDVI area The following cell generates a plot of the mean NDVI over time.  Note that we specify the dimensions along which to take the mean (```dim=['x', 'y']```).  It is important that we omitted the ```time``` dimension &mdash; if we had included it, our mean would collapse the data along _every_ dimension, returning only a single value (the mean value in x, y, _and_ time).<jupyter_code>NDVI.mean(dim=['x', 'y']).plot(size=6)<jupyter_output><empty_output><jupyter_text>Finally, we decorate the plot from the previous cell by adding a title and a label for our Y axis.<jupyter_code>NDVI.mean(dim=['x', 'y']).plot(size=6)
plt.title('Quarterly Trend in NDVI')
plt.ylabel('Mean NDVI')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11104536489151874, 'fraction_numerical': 0.017159763313609466, 'mean_word_length': 6.052851182197497, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '248349', 'license': 'Apache-2.0', 'n_tokens_mistral': 1725, 'n_tokens_neox': 1614, 'n_words': 588, 'path': 'docs/_static/session_4/Calculate_NDVI_Part_2_solution-annotated.ipynb', 'repo_name': 'digitalearthafrica/notebooks-training', 'size': '5070'}