Module 7
Caribbean Regional Node

To learn how bathymetry can be mapped crudely from digital imagery and the limitations of the results.
 
 

This lesson relates to material covered in Chapter 15 of the Remote Sensing Handbook for Tropical Coastal Management and readers are recommended to consult this for further details of the techniques involved. This lesson describes an empirical approach to mapping bathymetry using Landsat Thematic Mapper imagery of the Caicos Bank.

Depth information from remote sensing has been used to augment existing charts (Bullard, 1983; Pirazolli, 1985), assist in interpreting reef features (Jupp et al., 1985) and map shipping corridors (Benny and Dawson, 1983). However, it has not been used as a primary source of bathymetric data for navigational purposes (e.g. mapping shipping hazards). The major limitations are inadequate spatial resolution and lack of accuracy. Hazards to shipping such as emergent coral outcrops or rocks are frequently be much smaller than the sensor pixel and so will fail to be detected.

Among satellite sensors, the measurement of bathymetry can be expected to be best with Landsat TM data because that sensor detects visible light from a wider portion of the visible spectrum, in more bands, than other satellite sensors. Landsat bands 1, 2 and 3 are all useful in measuring bathymetry: so too is band 4 in very shallow (<1 m), clear water over bright sediment (where bands 1-3 tend to be saturated). Bands 5 and 7 are completely absorbed by even a few cm of water and should therefore be removed from the image before bathymetric calculations are performed.

This lesson will help you understand the importance of water depth in determining the reflectance of objects on the seabed and prepare you for understanding Lesson 5 which deals with how one attempts to compensate for these effects.

Familiarity with Bilko for Windows 2.0 is required to carry out this lesson. In particular, you will need experience of using Formula documents to carry out mathematical manipulations of images. Spreadsheet skills are also essential to benefit fully from the lesson. Instructions for spreadsheet manipulations are given for Excel and should be amended as appropriate for other spreadsheets.

The image you will be using was acquired by Landsat 5 TM on 22nd November 1990 at 14.55 hours Universal Time (expressed as a decimal time and thus equivalent to 14:33 GMT). The Turks & Caicos are on GMT - 5 hours so the overpass would have been at 09:33 local time. You are provided with bands 1 (blue), 2 (green), 3 (red) and 4 (near-infrared) of part of this image as the files TMBATHY1.GIF, TMBATHY2.GIF, TMBATHY3.GIF, and TMBATHY4.GIF. These images are of DN values but have been geometrically corrected. The post-correction pixel size is 33 x 33 m and each image is 800 pixels across and 900 pixels along-track.

Question: 4.1. How many kilometres wide and long is the image? Answers

The fundamental principle behind using remote sensing to map bathymetry is that different wavelengths of light will penetrate water to a varying degree. When light passes through water it becomes attenuated by interaction with the water column. The intensity of light remaining (I_d) after passage length p through water, is given by:

where I₀ = intensity of the incident light and k = attenuation coefficient, which varies with wavelength. Equation 4.1 can be made linear by taking natural logarithms:

Red light attenuates rapidly in water and does not penetrate further than about 5 m in clear water. By contrast, blue light penetrates much further and in clear water the seabed can reflect enough light to be detected by a satellite sensor even when the depth of water approaches 30 m. The depth of penetration is dependent on water turbidity. Suspended sediment particles, phytoplankton and dissolved organic compounds will all effect the depth of penetration (and so limit the range over which optical data may be used to estimate depth) because they scatter and absorb light, thus increasing attenuation.

There are two parts to Jupp's method i) the calculation of depth of penetration zones or DOP zones, (ii) the interpolation of depths within DOP zones.

The coefficient of attenuation k depends on wavelength. Longer wavelength light (red in the visible part of the spectrum) has a higher attenuation coefficient than short wavelengths (blue). Therefore red light is removed from white light passing vertically through water faster than is blue. There will therefore be a depth at which all the light detected by band 3 (visible red, 630-690 nm)) of the Landsat TM sensor has been attenuated (and which will therefore appear dark in a band 3 image). However not all wavelengths will have been attenuated to the same extent - there will still be shorter wavelength light at this depth which is detectable by bands 2 and 1 of the Landsat TM sensor (visible green, 520-600 nm; and visible blue, 450-520 nm respectively). This is most easily visualised if the same area is viewed in images consisting of different bands (Figure 4.1). The shallowest areas are bright even in band 4 but deeper areas (> about 15 m) are illuminated only in band 1 because it is just the shorter wavelengths which penetrate to the bottom sufficiently strongly to be reflected back to the satellite.

Figure 4.1. Band 1 and 4 Landsat TM images (DN) of the area around South Caicos. Deeper areas are illuminated by shorter wavebands. There is deep (30 - 2500 m) water in the bottom right hand corner of the images - this is not completely dark because some light is returned to the sensor over these areas by specular reflection and atmospheric scattering. Band 1 illuminates the deep fore reef areas (as indicated by the arrows). The edge of the reef or `drop-off' can be seen running towards the top right hand corner of the band 1 image. In areas of shallow (<1 m) clear water over white sand band 4 reveals some bathymetric detail (as indicated by the arrows).

Landsat TM Band 1 Landsat TM Band 4

The attenuation of light through water has been modelled by Jerlov (1976) and the maximum depth of penetration for different wavelengths calculated for various water types. Taking Landsat MSS as an example, it is known that green light (band 1, 0.50-0.60 m m) will penetrate to a maximum depth of 15 m in the waters of the Great Barrier Reef, red light (band 2, 0.60-0.70 m m) to 5 m, near infra-red (band 3, 0.70-0.80 m m) to 0.5 m and infra-red (band 4, 0.80-1.1 m m) is fully absorbed (Jupp, 1988).

Thus a "depth of penetration" zone can be created for each waveband (Figure 4.2). One zone for depths where green but not red can penetrate, one for where red but not near infra-red can penetrate, and one for where near-infrared (band 3) but not infra-red (band 4) can penetrate. Thus for the Landsat MSS example over clear Great Barrier Reef water three zones could be determined: 5-15 m, 0.5-5 m and surface to 0.5 m depth.

The method of Jupp (1988) has three critical assumptions:

UTM coordinate referenced field survey depth measurements for 515 sites on the Caicos Bank are provided in the Excel spreadsheet file DEPTHS.XLS. This includes the date, time, UTM coordinates from an ordinary or Differential Global Positioning System (GPS or DGPS), depth in metres (to nearest cm) measured using a hand-held echo-sounder, and DN values recorded in Landsat Thematic Mapper bands 1 to 4 at each coordinate (obtained using ERDAS Imagine 8.2 software).

The Tidal Prediction by the Admiralty Simplified Harmonic Method NP 159A Version 2.0 software of the UK Hydrographic Office was used to calculate tidal heights during each day of field survey so that all field survey depth measurements could be corrected to depths below datum (Lowest Astronomical Tide). [Admiralty tide tables could be used if this or similar software is not available]. The predicted tidal height at the time of the satellite overpass at approximately 09.30 h local time on 22 November 1990 was 0.65 m (interpolated between a predicted height of 0.60 m at 09.00 h and one of 0.70 m at 10.00 h). This height was added to all measured depths below datum to give the depth of water at the time of the satellite overpass.

The first task is to find out what the maximum DN values of deep water pixels are for bands 1-4 of the Landsat TM image. At the same time we will note the minimum values. If the image had been radiometrically and atmospherically corrected, the reflectance from deep water would be close to zero but because it has not although no light is being reflected back from the deep ocean the sensors will be recording significant DN values because of atmospheric scattering, etc. Knowing the highest DN values over deep water provides us with a baseline that we can use to determine the depth at which we are getting higher than background reflectance from the seabed. For each waveband, if DN values are greater than the maximum deep water values then we assume the sensor is detecting light reflected from the seabed.

The next stage is to use UTM coordinate referenced field survey data to estimate the maximum depths of penetration for each waveband. This can be done by inspecting the file DEPTHS.XLS which contains details of measured depths and DN values for 515 field survey sites. The survey sites have been sorted on depth with the deepest depths first. Table 4.1 lists the maximum depths of penetration recorded by Jupp (1988) for bands 1-4 of Landsat TM which gives us a starting point for searching for pixels near the maximum depth of penetration.

Note: The (D)GPS derived UTM coordinates for the depth measurements had positional errors in the order of 18-30 m for GPS coordinates and 2-4 m for DGPS ones. The Landsat TM pixels are 33 x 33 m in size. The reflectances recorded at the sensor are an average for the pixel area with contributions from scattering in the atmosphere. There is thus spatial and radiometric uncertainty in the data which needs to be considered. At the edge of the reef depth may change rapidly and spatial uncertainty in matching pixel DN values to depths is likely be most acute.

For Landsat TM1 you will find that the average depth of those pixels which have DN values > L_{deep max} is slightly less than that for pixels with values of 56-57. Given the small sample size and spatial and other errors this is understandable. The maximum depth of penetration for TM1 (z₁) presumably lies somewhere between these two depths and for our protocol we will take the average of the two depths as being our best estimate.

If none of these conditions apply, then pixels are coded to 0. A few pixels may have higher values in band 2 than band 1, which should not happen in theory. These will also be coded to zero and filtered out of the final depth image.

To implement this decision tree you will need to set up a Formula document with a series of conditional statements. This is primarily a lesson in remote sensing not computing so I have prepared one for you which creates four new images, showing the four DOP zones.

Calculating DOP zones does not assign a depth to each pixel, instead it assigns a pixel to a depth range (e.g. 13.5-20.8 m). Step two of Jupp's method involves interpolating depths for each pixel within each DOP zone. We will only work with DOP zone 2.

In DOP zone 2, the DN value of any submerged pixel in TM band 2 (L₂) can be expressed as:

where L_{2 deep mean} is the average deep water pixel value for TM band 2 calculated in Table 4.2, L_{2 surface} is the average DN value at the sea-surface (i.e. with no attenuation in the water column), k₂ is the attenuation coefficient for TM band 2 wavelengths through the water column, and z is the depth. Thus for a particular bottom type the DN value can vary between a minimum which is the mean deep water DN, and a maximum which is the DN which would be obtained if that bottom type was at the surface and no attenuation was occurring. In between the maximum depth of penetration for TM band 2 (z₂) and the surface, the DN value is purely a function of depth (z), with the rate of decrease in DN with depth being controlled by the attenuation coefficient for band 2 (k₂).

Equation 4.3 can be rewritten as:

If we define X₂ as follows: 

we can get rid of the exponential thus:

Now for a given bottom type in TM band 2,  will be a constant which to simplify the look of Equation 4.5, we will call A₂ to give the following linear regression relationship:

This gives us our basis for interpolating depths in each DOP zone. We will illustrate the rationale behind the interpolation continuing to take DOP zone 2 as an example. If a pixel is only just inside DOP zone 2 as opposed to DOP zone 1 (the zone with the next lowest coefficient of attenuation, Figure 4.2) then it has a value of X_{2 min} given by Equation 4.7 as:

The minimum TM band 2 DN in DOP zone 2 (L_{2 min}) will be one more than L_{deep max} for band 2, thus X_{2 min} is log_e([L_{2 deep max} +1] - L_{2 deep mean}). A₂ and the attenuation coefficient k₂ are specific to TM band 2, and z₂ is the maximum depth of penetration of TM band 2 (c. 13.5 m [Table 4.3] and the depth at which minimum reflectances are obtained in DOP zone 2).

Now consider another pixel which is only just inside DOP zone 2 as opposed to DOP zone 3 (the zone with the next highest coefficient of attenuation, Figure 4.2). This has a value of X_{2 max} = log_e(L_{2 max} - L_{2 deep mean}) defined by the equation:

where values are as for Equation 4.7 except that z₃ marks the upper end of DOP zone 2 (c. 4.2 m and the depth at which maximum reflectances are obtained in DOP zone 2).

The values of z₂ and z₃ are known (Table 4.3), L_{2 min} will generally be L_{2 deep max} + 1 (see Table 4.4 for L_{deep max} values for each band), and L_{2 max} and thus X_{2 max} can be determined from histograms of the appropriate images multiplied by their corresponding DOP zone masks (in this example, TMBATHY2.GIF multiplied by the DOP zone 2 mask image). For the purposes of this lesson we will just carry out the procedure for DOP zone 2.

The L_{i min} and L_{i max} values for DOP zones i = 1, 3 and 4 have already been entered in Table 4.5 using the same method, and k_i and A_i have already been calculated using the two equations below. For each band, we can calculate the attenuation coefficient because we know how much the DN value has changed over a known depth range, thus for TM band 1, the DN value has changed from 58 to 69 between 20.8 and 13.5 m depth (i.e. by 11 units in 7.3 m).

In our example for DOP zone 2, Equations 4.7 and 4.8 form a pair of simultaneous equations which can be solved for k₂ (Equation 4.9). Once this is know Equation 4.7 can be re-arranged to find A₂ (Equation 4.10).

Using the ground measurements of depth listed in the spreadsheet DEPTHS.XLS one can compare the depths predicted by the image at each UTM coordinate with those measured on site using a depth sounder. The results of plotting measured depth against predicted depth for the image made using the modified Jupp's method are presented in Figure 4.3.

This method produced a correlation of 0.82 between predicted and measured depth but required extensive field bathymetric data. Even this method does not produce bathymetric maps suitable for navigation since the average difference between depths predicted from imagery and ground-truthed depths ranged from about 0.8 m in shallow water (<2.5 m deep) to about 2.8 m in deeper water (2.5-20 m deep).

However, if you have understood this exercise in trying to predict depth from satellite imagery you will find it much easier to understand the next lesson which is all about trying to compensate for the effect of depth on bottom reflectance so that one can map submerged habitats.

1. The image is 26.4 km across and 29.7 km along track (north-south).
2. The assumption which is most likely not to be satisfied in the clear waters over the Caicos Bank is that the colour and hence reflective properties of the substrate are constant. The water quality is likely to be worse in rough areas where sediment is stirred up and so the attenuation coefficient for each band may also vary across the image.

3. The depth of penetration (z₁) for TM band 1 must lie somewhere between 18.70 and 25.35 m.
4. The average depth of pixels in this range which have DN values greater than L_{deep max} (57) is 21.7 m.
5. The average depths of pixels in this range which have values less than or equal to L_{deep max} (57) is 19.9 m.

6. These shallow water anomalies (areas which seem to be very shallow water on the image but shouldn't be) are caused by clouds over the area. These clouds are particularly visible in the near-infrared.
7. These deepwater anomalies (areas which seem to be deep water on the image but which common sense says cannot be) are caused by dense seagrass beds which have a very low reflectance right across the spectrum and thus, even though shallower than the depth of penetration for TM band 1, reflect no more light than deep water.
8. 14320 TM2 pixels in DOP zone 2 have values = 17 (L_{2 min}).
9. The highest DN value of the TM2 pixels in DOP zone 2 is 41.
10. The top 0.1% of pixels are highlighted at a DN value of 37, so your estimate of L_{2 min} should be 37.

11. It is difficult to judge fine gradations in colour such as the BATHYBLU.PAL palette displays. So the BATHYCON.PAL palette conveys the most information, allowing you to assign areas to depth zones very easily.
12. The modal depth (in the 20 by 20 block with coordinates 220, 800 as its north-west corner) is 1.6 m. Mode is at a pixel value of 16; it includes 24.8% of pixels in the block.
13. The modal depth (in the 20 by 20 block with coordinates 650, 55 as its north-west corner) is 0.8 m. Mode is at a pixel value of 8; it includes 21.0% of pixels in the block.
14. The modal depth (in the 20 by 20 block with coordinates 395, 480 as its north-west corner) is 18.7 m. Mode is at a pixel value of 187; it includes 28% of pixels in the block.
15. The modal depth (in the 20 by 20 block with coordinates 100, 515 as its north-west corner) is 8.3 m. Mode is at a pixel value of 83; it includes 52.5% of pixels in the block.

Back to main page

Back to Module 7

USDOC | NOAA | NESDIS | CoastWatch