diff --git a/episodes/09-raster-calculations.md b/episodes/09-raster-calculations.md
index 2ca0f32a..608f019d 100644
--- a/episodes/09-raster-calculations.md
+++ b/episodes/09-raster-calculations.md
@@ -10,42 +10,43 @@ exercises: 15
 
 :::objectives
 - Carry out operations with two rasters using Python's built-in math operators.
-- Reclassify a continuous raster to a categorical raster.
 :::
 
-
-
 ## Introduction
 
 We often want to combine values of and perform calculations on rasters to create a new output raster. This episode
 covers how to perform basic math operations using raster datasets. It also illustrates how to match rasters with
-different resolutions so that they can be used in the same calculation. As an example, we will calculate a vegetation
-index over one of the satellite scenes.
+different resolutions so that they can be used in the same calculation. As an example, we will calculate
+[a binary classification mask](https://custom-scripts.sentinel-hub.com/sentinel-2/burned_area_ms/) to identify burned
+area over a satellite scene.
 
-### Normalized Difference Vegetation Index (NDVI)
-Suppose we are interested in monitoring vegetation fluctuations using satellite remote sensors. Scientists have defined
-a vegetation index to quantify the amount of green leaf vegetation using the light reflected in different wavelengths.
-This index, named Normalized Difference Vegetation Index (NDVI), exploits the fact that healthy green leaves strongly
-absorb red visible light while they mostly reflect light in the near infrared (NIR). The NDVI is computed as:
+The classification mask requires the following of [the Sentinel-2 bands](https://gisgeography.com/sentinel-2-bands-combinations/)
+(and derived indices):
+
+* [Normalized difference vegetation index (NDVI)](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index),
+  derived from the **near-infrared (NIR)** and **red** bands:
 
 $$ NDVI = \frac{NIR - red}{NIR + red} $$
 
-where $NIR$ and $red$ label the reflectance values of the corresponding wavelengths. NDVI values range from -1 to +1.
-Values close to one indicate high density of green leaves. Poorly vegetated areas typically have NDVI values close to
-zero. Negative NDVI values often indicate cloud and water bodies.
+* [Normalized difference water index (NDWI)](https://en.wikipedia.org/wiki/Normalized_difference_water_index), derived
+  from the **green** and **NIR** bands:
 
-![Source: Wu C-D, McNeely E, Cedeño-Laurent JG, Pan W-C, Adamkiewicz G, Dominici F, et al. (2014) Linking Student Performance in Massachusetts Elementary Schools with the “Greenness” of School Surroundings Using Remote Sensing. PLoS ONE 9(10): e108548. https://doi.org/10.1371/journal.pone.0108548](fig/E09/PONE-NDVI.jpg){alt="PONE-NDVI image"}
+$$ NDWI = \frac{green - NIR}{green + NIR} $$
 
-:::callout
-## More Resources
-Check out more on NDVI in the NASA Earth Observatory portal:
-[Measuring Vegetation](https://earthobservatory.nasa.gov/features/MeasuringVegetation/measuring_vegetation_2.php).
-:::
+* A custom index derived from two of the **short-wave infrared (SWIR)** bands (with wavelenght ~1600 nm and ~2200 nm,
+  respectively):
+
+$$ INDEX = \frac{SWIR_{16} - SWIR_{22}}{SWIR_{16} + SWIR_{22}}$$
+
+* The **blue**, **NIR**, and **SWIR** (1600 nm) bands.
+
+In the following, we start by computing the NDVI.
 
 ### Load and crop the Data
 For this episode, we will use one of the Sentinel-2 scenes that we have already employed in the previous episodes.
 
 :::callout
+
 ## Introduce the Data
 
 We'll continue from the results of the satellite image search that we have carried out in an exercise from
@@ -56,36 +57,43 @@ raster data using this [link](https://figshare.com/ndownloader/files/36028100).
 file in your current working directory, and extract the archive file by double-clicking on it or by running the
 following command in your terminal `tar -zxvf geospatial-python-raster-dataset.tar.gz`. Use the file `geospatial-python-raster-dataset/search.json`
 (instead of `search.json`) to get started with this lesson.
+
 :::
 
 Let's load the results of our initial imagery search using `pystac`:
 
 ```python
 import pystac
-items = pystac.ItemCollection.from_file("search.json")
+items = pystac.ItemCollection.from_file("rhodes_sentinel-2.json")
 ```
 
-
-We then select the second item, and extract the URIs of the red and NIR bands ("red" and "nir08", respectively):
+We then select the first item, and extract the URLs of the red and NIR bands ("red" and "nir", respectively):
 
 ```python
-red_uri = items[1].assets["red"].href
-nir_uri = items[1].assets["nir08"].href
+item = items[0]
+red_href = item.assets["red"].get_absolute_href()
+nir_href = item.assets["nir"].get_absolute_href()
 ```
 
 Let's load the rasters with `open_rasterio` using the argument `masked=True`.
 
 ```python
 import rioxarray
-red = rioxarray.open_rasterio(red_uri, masked=True)
-nir = rioxarray.open_rasterio(nir_uri, masked=True)
+red = rioxarray.open_rasterio(red_href, masked=True)
+nir = rioxarray.open_rasterio(nir_href, masked=True)
 ```
 
-Let's also restrict our analysis to the same crop field area defined in the previous episode by clipping the rasters
-using a bounding box:
+Let’s also restrict our analysis to the island of Rhodes - we can extract the bounding box from the vector file written
+in an earlier episode (note that we need to match the CRS to the one of the raster files):
 
 ```python
-bbox = (629_000, 5_804_000, 639_000, 5_814_000)
+# determine bounding box of Rhodes, in the projected CRS
+import geopandas
+rhodes = geopandas.read_file('rhodes.gpkg')
+rhodes_reprojected = rhodes.to_crs(red.rio.crs)
+bbox = rhodes_reprojected.total_bounds
+
+# crop the rasters
 red_clip = red.rio.clip_box(*bbox)
 nir_clip = nir.rio.clip_box(*bbox)
 ```
@@ -105,65 +113,49 @@ nir_clip.plot(robust=True)
 
 ![](fig/E09/NIR-band.png){alt="near infra-red band image"}
 
-It is immediately evident how crop fields (rectangular shapes in the central part of the two figures) appear as dark and
-bright spots in the red-visible and NIR wavelengths, respectively, suggesting the presence of leafy crop at the time of
-observation (end of March). The same fields would instead appear as dark spots in the off season.
+The burned area is immediately evident as a dark spot in the NIR wavelength, due to the lack of reflection from the
+vegetation in the scorched area.
 
 ## Raster Math
+
 We can perform raster calculations by subtracting (or adding,
 multiplying, etc.) two rasters. In the geospatial world, we call this
 "raster math", and typically it refers to operations on rasters that
 have the same width and height (including `nodata` pixels).
 We can check the shapes of the two rasters in the following way:
+
 ```python
 print(red_clip.shape, nir_clip.shape)
 ```
 
 ```output
-(1, 1000, 1000) (1, 500, 500)
+(1, 1131, 1207) (1, 1131, 1207)
 ```
 
-Both rasters include a single band, but their width and height do not match.
-We can now use the `reproject_match` function, which both reprojects and clips
-a raster to the CRS and extent of another raster.
-
-```python
-red_clip_matched = red_clip.rio.reproject_match(nir_clip)
-print(red_clip_matched.shape)
-```
-
-```output
-(1, 500, 500)
-```
+The shapes of the two rasters match (and, not shown, the coordinates and the CRSs match too).
 
 Let's now compute the NDVI as a new raster using the formula presented above.
-We'll use `rioxarray` objects so that we can easily plot our result and keep
+We'll use `DataArray` objects so that we can easily plot our result and keep
 track of the metadata.
 
 ```python
-ndvi = (nir_clip - red_clip_matched)/ (nir_clip + red_clip_matched)
+ndvi = (nir_clip - red_clip)/ (nir_clip + red_clip)
 print(ndvi)
 ```
 
 ```output
-<xarray.DataArray (band: 1, y: 500, x: 500)>
-array([[[ 0.7379576 ,  0.77153456,  0.54531944, ...,  0.39254385,
-          0.49227372,  0.4465174 ],
-        [ 0.7024894 ,  0.7074668 ,  0.3903298 , ...,  0.423283  ,
-          0.4706971 ,  0.45964912],
-        [ 0.6557818 ,  0.5610572 ,  0.46742022, ...,  0.4510345 ,
-          0.43815723,  0.6005133 ],
+<xarray.DataArray (band: 1, y: 1131, x: 1207)>
+array([[[0., 0., 0., ..., 0., 0., 0.],
+        [0., 0., 0., ..., 0., 0., 0.],
+        [0., 0., 0., ..., 0., 0., 0.],
         ...,
-        [ 0.02391171,  0.21843003,  0.02479339, ..., -0.50923485,
-         -0.53367877, -0.4955414 ],
-        [ 0.11376493,  0.17681159, -0.1673566 , ..., -0.5221932 ,
-         -0.5271318 , -0.4852753 ],
-        [ 0.45398772, -0.00518135,  0.03346133, ..., -0.5019455 ,
-         -0.4987013 , -0.49081364]]], dtype=float32)
+        [0., 0., 0., ..., 0., 0., 0.],
+        [0., 0., 0., ..., 0., 0., 0.],
+        [0., 0., 0., ..., 0., 0., 0.]]], dtype=float32)
 Coordinates:
   * band         (band) int64 1
-  * x            (x) float64 6.29e+05 6.29e+05 6.29e+05 ... 6.39e+05 6.39e+05
-  * y            (y) float64 5.814e+06 5.814e+06 ... 5.804e+06 5.804e+06
+  * x            (x) float64 5.615e+05 5.616e+05 ... 6.097e+05 6.098e+05
+  * y            (y) float64 4.035e+06 4.035e+06 4.035e+06 ... 3.99e+06 3.99e+06
     spatial_ref  int64 0
 ```
 
@@ -178,7 +170,7 @@ ndvi.plot()
 Notice that the range of values for the output NDVI is between -1 and 1.
 Does this make sense for the selected region?
 
-Maps are great, but it can also be informative to plot histograms of values to better understand the distribution. We can accomplish this using a built-in xarray method we have already been using: `plot`
+Maps are great, but it can also be informative to plot histograms of values to better understand the distribution. We can accomplish this using a built-in xarray method we have already been using: `plot.hist()`
 
 ```python
 ndvi.plot.hist()
@@ -187,196 +179,135 @@ ndvi.plot.hist()
 ![](fig/E09/NDVI-hist.png){alt="NDVI histogram"}
 
 :::challenge
-## Exercise: Explore NDVI Raster Values
 
-It's often a good idea to explore the range of values in a raster dataset just like we might explore a dataset that we collected in the field. The histogram we just made is a good start but there's more we can do to improve our understanding of the data.
+### Exercise: NDWI and custom index to detect burned areas
 
-1. What is the min and maximum value for the NDVI raster (`ndvi`) that we just created? Are there missing values?
-2. Plot a histogram with 50 bins instead of 8. What do you notice that wasn't clear before?
-3. Plot the `ndvi` raster using breaks that make sense for the data.
+Calculate the other two indices required to compute the burned area classification mask, specifically:
+
+* The [normalized difference water index (NDWI)](https://en.wikipedia.org/wiki/Normalized_difference_water_index), derived from the **green** and **NIR** bands (with band ID "green" and "nir", respectively):
+
+$$ NDWI = \frac{green - NIR}{green + NIR} $$
+
+* A custom index derived from the  1600 nm and the 2200 nm **short-wave infrared (SWIR)** bands ( "swir16" and "swir22", respectively):
+
+$$ INDEX = \frac{SWIR_{16} - SWIR_{22}}{SWIR_{16} + SWIR_{22}}$$
+
+What challenge do you foresee in combining the data from the two indices?
 
 ::::solution
 
-1. Recall, if there were nodata values in our raster,
-we would need to filter them out with `.where()`. Since we have loaded the rasters with `masked=True`, missing
-values are already encoded as `np.nan`'s. The `ndvi` array actually includes a single missing value.
 ```python
-print(ndvi.min().values)
-print(ndvi.max().values)
-print(ndvi.isnull().sum().values)
+def get_band_and_clip(item, band_name, bbox):
+    href = item.assets[band_name].get_absolute_href()
+    band = rioxarray.open_rasterio(href, masked=True)
+    return band.rio.clip_box(*bbox)
 ```
 
-```output
--0.99864775
-0.9995788
-1
-```
 
-2. Increasing the number of bins gives us a much clearer view of the distribution. Also, there seem to be very few
-NDVI values larger than ~0.9.
 ```python
-ndvi.plot.hist(bins=50)
+green_clip = get_band_and_clip(item, 'green', bbox)
+swir16_clip = get_band_and_clip(item, 'swir16', bbox)
+swir22_clip = get_band_and_clip(item, 'swir22', bbox)
 ```
 
-![](fig/E09/NDVI-hist-bins.png){alt="NDVI histogram with 50 bins"}
 
-3. We can discretize the color bar by specifying the intervals via the `levels` argument to `plot()`.
-Suppose we want to bin our data in the following intervals:
-* $-1 \le NDVI \lt 0$ for water;
-* $0 \le NDVI \lt 0.2$ for no vegetation;
-* $0.2 \le NDVI \lt 0.7$ for sparse vegetation;
-* $0.7 \le NDVI \lt 1$ for dense vegetation.
+```python
+ndwi = (green_clip - nir_clip)/(green_clip + nir_clip)
+index = (swir16_clip - swir22_clip)/(swir16_clip + swir22_clip)
+```
 
 ```python
-class_bins = (-1, 0., 0.2, 0.7, 1)
-ndvi.plot(levels=class_bins)
+ndwi.plot(robust=True)
 ```
 
-![](fig/E09/NDVI-map-binned.png){alt="binned NDVI map"}
-::::
-:::
+![](fig/E09/NDWI.png){alt="NDWI index"}
 
-Missing values can be interpolated from the values of neighbouring grid cells using the `.interpolate_na` method. We
-then save `ndvi` as a GeoTiff file:
 ```python
-ndvi_nonan = ndvi.interpolate_na(dim="x")
-ndvi_nonan.rio.to_raster("NDVI.tif")
+index.plot(robust=True)
 ```
 
+![](fig/E09/custom-index.png){alt="custom index"}
 
-## Classifying Continuous Rasters in Python
+The challenge in combining the different indices is that the SWIR bands (and thus the derived custom index) have lower resolution:
 
-Now that we have a sense of the distribution of our NDVI raster, we
-can reduce the complexity of our map by classifying it. Classification involves
-assigning each pixel in the raster to a class based on its value. In Python, we
-can accomplish this using the `numpy.digitize` function.
+```python
+ndvi.rio.resolution(), ndwi.rio.resolution(), index.rio.resolution()
+```
 
-First, we define NDVI classes based on a list of values, as defined in the last exercise:
-`[-1, 0., 0.2, 0.7, 1]`. When bins are ordered from
-low to high, as here, `numpy.digitize` assigns classes like so:
+```output
+((10.0, -10.0), (10.0, -10.0), (20.0, -20.0))
+```
+
+::::
 
-![Source: Image created for this lesson ([license](../LICENSE.md))](fig/E09/NDVI-classes.jpg){alt="NDVI classes"}
+:::
 
+In order to combine data from the computed indices, we use the `reproject_match` method, which reprojects, clips and
+match the resolution of a raster using another raster as a template. We use the `index` raster as a template, and match
+`ndvi` and `ndwi` to its resolution and extent:
 
-Note that, by default, each class includes the left but not the right bound. This is not an issue here, since the
-computed range of NDVI values is fully contained in the open interval (-1; 1) (see exercise above).
 
 ```python
-import numpy as np
-import xarray
-
-# Defines the bins for pixel values
-class_bins = (-1, 0., 0.2, 0.7, 1)
-
-# The numpy.digitize function returns an unlabeled array, in this case, a
-# classified array without any metadata. That doesn't work--we need the
-# coordinates and other spatial metadata. We can get around this using
-# xarray.apply_ufunc, which can run the function across the data array while
-# preserving metadata.
-ndvi_classified = xarray.apply_ufunc(
-    np.digitize,
-    ndvi_nonan,
-    class_bins
-)
+ndvi_match = ndvi.rio.reproject_match(index)
+ndwi_match = ndwi.rio.reproject_match(index)
 ```
 
-Let's now visualize the classified NDVI, customizing the plot with proper title and legend. We then export the
-figure in PNG format:
+Finally, we also fetch the blue band data and match this and the NIR band data to the template:
 
 ```python
-import earthpy.plot as ep
-import matplotlib.pyplot as plt
-
-from matplotlib.colors import ListedColormap
-
-# Define color map of the map legend
-ndvi_colors = ["blue", "gray", "green", "darkgreen"]
-ndvi_cmap = ListedColormap(ndvi_colors)
-
-# Define class names for the legend
-category_names = [
-    "Water",
-    "No Vegetation",
-    "Sparse Vegetation",
-    "Dense Vegetation"
-]
-
-# We need to know in what order the legend items should be arranged
-category_indices = list(range(len(category_names)))
-
-# Make the plot
-im = ndvi_classified.plot(cmap=ndvi_cmap, add_colorbar=False)
-plt.title("Classified NDVI")
-# earthpy helps us by drawing a legend given an existing image plot and legend items, plus indices
-ep.draw_legend(im_ax=im, classes=category_indices, titles=category_names)
-
-# Save the figure
-plt.savefig("NDVI_classified.png", bbox_inches="tight", dpi=300)
+blue_clip = get_band_and_clip(item, 'blue', bbox)
+blue_match = blue_clip.rio.reproject_match(index)
+nir_match = nir_clip.rio.reproject_match(index)
 ```
 
-![](fig/E09/NDVI-classified.png){alt="classified NDVI map"}
-
-We can finally export the classified NDVI raster object to a GeoTiff file. The `to_raster()` function
-by default writes the output file to your working directory unless you specify a
-full file path.
+We can now go ahead and compute the binary classification mask for burned areas. Note that we need to convert the unit
+of the Sentinel-2 bands [from digital numbers to reflectance](https://docs.sentinel-hub.com/api/latest/data/sentinel-2-l2a/#units)
+(this is achieved by dividing by 10,000):
 
 ```python
-ndvi_classified.rio.to_raster("NDVI_classified.tif", dtype="int32")
+burned = (
+    (ndvi_match <= 0.3) &
+    (ndwi_match <= 0.1) &
+    ((index + nir_match/10_000) <= 0.1) &
+    ((blue_match/10_000) <= 0.1) &
+    ((swir16_clip/10_000) >= 0.1)
+)
 ```
 
-:::challenge
-## Exercise: Compute the NDVI for the Texel island
-
-Data are often more interesting and powerful when we compare them across various
-locations. Let's compare the computed NDVI map with the one of another region in the same Sentinel-2 scene:
-the [Texel island](https://en.wikipedia.org/wiki/Texel), located in the North Sea.
-
-0. You should have the red- and the NIR-band rasters already loaded (`red` and `nir` variables, respectively).
-1. Crop the two rasters using the following bounding box: `(610000, 5870000, 630000, 5900000)`. Don't forget to check the shape of the data, and make sure the cropped areas have the same CRSs, heights and widths.
-2. Compute the NDVI from the two raster layers and check the max/min values to make sure the data
-is what you expect.
-3. Plot the NDVI map and export the NDVI as a GeoTiff.
-4. Compare the distributions of NDVI values for the two regions investigated.
+The classification mask has a single element along the "band" axis, we can drop this dimension in the following way:
 
-::::solution
-1) We crop the area of interest using `clip_box`:
 ```python
-bbox_texel = (610000, 5870000, 630000, 5900000)
-nir_texel = nir.rio.clip_box(*bbox_texel)
-red_texel = red.rio.clip_box(*bbox_texel)
+burned = burned.squeeze()
 ```
 
-2) Reproject and clip one raster to the extent of the smaller raster using `reproject_match`. The lines of code below
-assign a variable to the reprojected raster and calculate the NDVI.
+Let's now fetch and visualize the true color image of Rhodes, after coloring the pixels identified as burned area in red:
 
 ```python
-red_texel_matched = red_texel.rio.reproject_match(nir_texel)
-ndvi_texel = (nir_texel - red_texel_matched)/ (nir_texel + red_texel_matched)
+visual_href = item.assets['visual'].get_absolute_href()
+visual = rioxarray.open_rasterio(visual_href)
+visual_clip = visual.rio.clip_box(*bbox)
 ```
 
-3) Plot the NDVI and save the raster data as a GeoTIFF file.
 
 ```python
-ndvi_texel.plot()
-ndvi_texel.rio.to_raster("NDVI_Texel.tif")
+# set red channel to max (255), green and blue channels to min (0).
+visual_clip[0] = visual_clip[0].where(~burned, 255)
+visual_clip[1:3] = visual_clip[1:3].where(~burned, 0)
 ```
 
-![](fig/E09/NDVI-map-Texel.png){alt="NDVI map Texel"}
-
-4) Compute the NDVI histogram and compare it with the region that we have previously investigated. Many more grid
-cells have negative NDVI values, since the area of interest includes much more water. Also, NDVI values close to
-zero are more abundant, indicating the presence of bare ground (sand) regions.
 ```python
-ndvi_texel.plot.hist(bins=50)
+visual_clip.plot.imshow()
 ```
 
-![](fig/E09/NDVI-hist-Texel.png){alt="NDVI histogram Texel"}
+![](fig/E09/visual-burned-index.png){alt="RGB image with burned area in red"}
 
-::::
-:::
 
+We can save the burned classification mask to disk after converting booleans to integers:
+
+```python
+burned.rio.to_raster('burned.tif', dtype='int8')
+```
 
 :::keypoints
 - Python's built-in math operators are fast and simple options for raster math.
-- numpy.digitize can be used to classify raster values in order to generate a less complicated map.
 :::
diff --git a/episodes/fig/E09/NDVI-classes.jpg b/episodes/fig/E09/NDVI-classes.jpg
deleted file mode 100644
index a2ec0eb1..00000000
Binary files a/episodes/fig/E09/NDVI-classes.jpg and /dev/null differ
diff --git a/episodes/fig/E09/NDVI-classified.png b/episodes/fig/E09/NDVI-classified.png
deleted file mode 100644
index 90c52e5d..00000000
Binary files a/episodes/fig/E09/NDVI-classified.png and /dev/null differ
diff --git a/episodes/fig/E09/NDVI-hist-Texel.png b/episodes/fig/E09/NDVI-hist-Texel.png
deleted file mode 100644
index b972e94b..00000000
Binary files a/episodes/fig/E09/NDVI-hist-Texel.png and /dev/null differ
diff --git a/episodes/fig/E09/NDVI-hist-bins.png b/episodes/fig/E09/NDVI-hist-bins.png
deleted file mode 100644
index c06e06b9..00000000
Binary files a/episodes/fig/E09/NDVI-hist-bins.png and /dev/null differ
diff --git a/episodes/fig/E09/NDVI-hist.png b/episodes/fig/E09/NDVI-hist.png
index 1ecbcc39..7b16ad62 100644
Binary files a/episodes/fig/E09/NDVI-hist.png and b/episodes/fig/E09/NDVI-hist.png differ
diff --git a/episodes/fig/E09/NDVI-map-Texel.png b/episodes/fig/E09/NDVI-map-Texel.png
deleted file mode 100644
index 68a2ecc6..00000000
Binary files a/episodes/fig/E09/NDVI-map-Texel.png and /dev/null differ
diff --git a/episodes/fig/E09/NDVI-map-binned.png b/episodes/fig/E09/NDVI-map-binned.png
deleted file mode 100644
index 06266ea7..00000000
Binary files a/episodes/fig/E09/NDVI-map-binned.png and /dev/null differ
diff --git a/episodes/fig/E09/NDVI-map.png b/episodes/fig/E09/NDVI-map.png
index 4bde6c21..635cac5d 100644
Binary files a/episodes/fig/E09/NDVI-map.png and b/episodes/fig/E09/NDVI-map.png differ
diff --git a/episodes/fig/E09/NDWI.png b/episodes/fig/E09/NDWI.png
new file mode 100644
index 00000000..fd025c1b
Binary files /dev/null and b/episodes/fig/E09/NDWI.png differ
diff --git a/episodes/fig/E09/NIR-band.png b/episodes/fig/E09/NIR-band.png
index 73658e9e..827d2642 100644
Binary files a/episodes/fig/E09/NIR-band.png and b/episodes/fig/E09/NIR-band.png differ
diff --git a/episodes/fig/E09/PONE-NDVI.jpg b/episodes/fig/E09/PONE-NDVI.jpg
deleted file mode 100644
index 030d8d79..00000000
Binary files a/episodes/fig/E09/PONE-NDVI.jpg and /dev/null differ
diff --git a/episodes/fig/E09/custom-index.png b/episodes/fig/E09/custom-index.png
new file mode 100644
index 00000000..92b805a1
Binary files /dev/null and b/episodes/fig/E09/custom-index.png differ
diff --git a/episodes/fig/E09/red-band.png b/episodes/fig/E09/red-band.png
index fb3a5011..b0ea4aba 100644
Binary files a/episodes/fig/E09/red-band.png and b/episodes/fig/E09/red-band.png differ
diff --git a/episodes/fig/E09/visual-burned-index.png b/episodes/fig/E09/visual-burned-index.png
new file mode 100644
index 00000000..9f9fd6ac
Binary files /dev/null and b/episodes/fig/E09/visual-burned-index.png differ