Calculating NDVI: Part 1¶
This exercise follows on from the previous section. In the final exercise of the previous session, you constructed a notebook to create geomedian composite.
In this section, we will create a new notebook based on the notebook from the previous section. Most of the code will remain unchanged, but we will change the area of interest and time extent. We will also add steps to resample the new dataset and create a geomedian. In the next section, we will calculate and plot NDVI.
We will be using the notebook we created in the previous section, Create a geomedian composite. If you have not already set up a copy of the notebook called
Geomedian_composite.ipynb with the required packages and functions, follow the instructions in previous section. Ensure you have completed all the steps, including loading the Sentinel-2 dataset.
Set up notebook¶
Create a copy of the notebook¶
Before you continue with the next step,
Log in to the Sandbox and open the Training folder.
Make a copy of the
Rename the notebook to
See create a copy of a notebook and rename it for more details.
Clearing the notebook¶
We will need to remove any output from previous runs of the notebook.
Select Kernel -> Restart Kernel and Clear All Outputs….
When prompted, select Restart.
Running the notebook¶
This notebook is still set up to run the Session 3 exercise, so you will need to follow the instructions below to modify it. Work cell by cell and pay attention to what needs to be changed.
Run the first cell, which contains the packages and functions for the analysis. No need to change anything here.
dc = datacube.Datacubecommand, change the app name to
"Calculate_ndvi". It should look like:
dc = datacube.Datacube(app="Calculate_ndvi")
Load the data¶
Change the x and y values to those shown below and run the cell.
x=(-6.1495, -6.1380) y=(13.9182, 13.9111)
Change the time in the
Remove the option
If you completed steps 2, 3 and 4, your load cell should look like
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'], resolution=(-10, 10), group_by='solar_day')
Run the cell. The load should return 71 time steps.
The fifth cell of the notebook contains an
rgb command to plot the loaded data. To match our example below, modify this cell so that it matches the code below:
timesteps = [1, 6, 8] rgb(sentinel_2_ds, bands=['red', 'green', 'blue'], index=timesteps)
This will plot images for the 1st, 6th and 8th timestep of the loaded data (remember that Python starts counting at 0). Your image should match the one below.
You may also like to run this cell a few times, experimenting with different values for the
timesteps parameter. The load command should have returned 71 time steps, meaning the values in your
timesteps list can be anywhere from
Resampling the dataset¶
Resampling is used to create a new set of times at regular intervals. Using the resample method, the data can be arranged in days, months, quarterly (three months) or yearly.
Below gives examples of how the data are grouped.
'nD'- number of days (e.g.
'7D'for seven days)
'nM'- number of months (e.g.
'6M'for six months)
'nY'- number of years (e.g.
'2Y'for two years)
Follow the steps below to resample the dataset time steps to quarterly.
Delete the code for plotting all RGB images:
rgb(sentinel_2_ds, bands=['red', 'green', 'blue'], col='time', size=4)
In the cleared cell, write the following code to resample the data and store it in the
resample_sentinel_2_ds = sentinel_2_ds.resample(time='3MS')
resample_sentinel_2_ds describes how to group the data into quarterly segments. We can now use this to calculate the geomedian for each quarterly segment.
S at the end of
'3MS' is to group the data by start of the month.
Compute the geomedian¶
For this session, instead of calling
xr_geomedian(sentinel_2_ds) on the entire array as we did in the previous exercise, we pass the
xr_geomedian function to
map and apply it separately to each resampled group(
Replace the existing
xr_geomedian code with:
geomedian_resample = resample_sentinel_2_ds.map(xr_geomedian)
We can plot the output geomedians, and see the change in the landscape over the year. Replace the exsiting
rgb code with:
rgb(geomedian_resample, bands=['red', 'green', 'blue'], col="time", col_wrap=4)
Comparing the two datasets¶
Comparing the two dataset can you tell the difference from the results shown below?
Take a look at the dimensions. The raw dataset
sentinel_2_ds has 71 time steps loaded, but after resampling to quarterly as in
geomedian_resample, the time dimension is now 4. This makes sense because the year has been divided into four quarters.
Congratulations! You have successfully modified a notebook to create a quarterly geomedian composite by resampling Sentinel-2 data.
If you’d like to experiment futher, try running the code with different resampling values. Did you learn anything interesting to share with us?
In the next section, we will continue working with this notebook to calculate the NDVI values for each of our quarterly geomedians.