Please see the following steps presenting the way I used space agencies data for this project
Step 01 Data Selection and Download
We approached to the Copernicus Open Access Hub using the following weblink (https://scihub.copernicus.eu/) to visually inspect and select suitable image for this application. As mentioned we selected Sentinel -2 images for this application. First criteria for image selection was less cloud coverage which can be automatically selected from the platform. Our last two criteria was clear waster and sun glint free image. Using Initial visual inspection I had selected the following images.
https://drive.google.com/file/d/1P0vzWMYOJk1UZ7KLLXQztI_9zLRMtooo/view?usp=sharing
Images are presented with their captured date. However it was almost impossible to have sun glint and cloud free images on the test site. Based on the subjective judgement Sentinel -2 images captured over Saint Martin island on 8th May, 2020 were selected for the following procedure. Its important to mention that the downloaded Sentinel 2 data product from Multispectral Instrument was not atmospherically corrected and it is a Top of Atmosphere Reflectance L1c data product
Step 02:Atmospheric Correction and Image Subset
Various atmospheric correction modules are available. Somme of them are Sen2Cor, Acolite, OPERA, 6s/MODTRAN LUTs etc. In this research we used Sen2Cor processor to generate Level 2A data product. Sen2Cor performs atmospheric, terrain and cirrus correction for Top-Of-Atmosphere Level 1C input data and produce Bottom of Atmosphere Reflectance's L2a data product. The process is implemented on command line interface on SNAP software. After applying atmospheric correction a subset image over the test site is produced using Sentinel B2,B3,B4,B8 spectral band. A true color composition of atmospherically corrected Sentinel level 2a data product should look like the following image.
https://drive.google.com/file/d/1wiI8Z0EroGtnc83sKjulpUtiavNignlZ/view?usp=sharing
Step 03. Land Mask
Land mask layer was created using NIR band. Threshold value for land and water was determined to separate the land from sea. After applying the land mask to RGB bands the sea areas should look like the following image.
https://drive.google.com/file/d/1xZQVl5efGHIfWXxugfmFemPXqIMcgTSn/view?usp=sharing
Step 04. Sun Glint Correction
This type of phenomena occurred in the satellite images when there is an specular reflection of the sun on water surfaces. Corrections are applied to removes sun glint reflected waves by using NIR band. The process improve visual appearance of the data product. Its important to apply such correction as our image is affected by this type of contamination. There are several sun glint removal algorithms available, In this exercise I used the one developed by Headley et al.,(2005). Deglint is implement on spreadsheet using regression analysis.
The equation described by Hedley et al. (2005) for the deglinting of a multispectral image is: R’i = Ri–bi(RNIR–MinNIR)
Where R’i is the deglinted pixel in band i
Ri is the reflectance from visible band i
bi is the regression slope
RNIR is the NIR band value
MinNIR the minimum NIR value of the sample
Theoretically such procedure tries to scale the relationship between sun glint and NIR signals by using one more sample images of the affected area. Linear regression between visible and NIR is performed over deep water. “The slope of the regression then used to predict the brightness in the visible band that would be expected if those pixels had a NIR value of MinNIR”. A graphical presentation of this correction techniques (Headley et al,2005) are presented on the following link.
https://drive.google.com/file/d/1CiyJB4x3_XJ1VCzW5ULmwvAPkc_8FINB/view?usp=sharing
I applied the following statistics on visible bands to produce deglinted image
Blue Band
Reference band: NIR
Min NIR Computed:0.013600
MinNIRUsed:0.013600
Regression slope = 0.579346
Regression R - squared= 0.375671
Green band
Reference band: NIR
Min NIR Computed:0.013600
Regression slope = 0.506637
Regression R - squared= 0.176810
Red band
Reference band: NIR
Min NIR Computed:0.013600
Regression slope = 0.770151
Regression R - squared= 0.667050
https://drive.google.com/file/d/1atblv9IYQqy65Kii0qwu22KJBoTccllD/view?usp=sharing
Images on the left shows original RGB composition and deglinted image on the right
Step 05. Dark Object Substraction
Substracting the atmospherically offset from every pixel in the band can generate the corrected image. As this offset value is different in each band therefore the histogram cut-off point is applied separately on each band.
Step 06. GEBCO 2020 Grid and Satellite Derived Bathymetry
As discussed earlier Multi mission compilation from GEBCO 2020 Grid is used for this study. Minor level stretching were applied on GEBCO 2020 Grid to present the depth variations at the study area. Initially 250 sample points were randmoly selected as shown on the following figure to simulate Ratio based bathymetry model. At later stage based on using visual inspection 50 spot heights were omitted from the analysis due to its inconsistency. Finally 200 sample points were used to drive linear based bathymetry model. Following figure shows GEBCO 2020 Grid with the extracted sample depths. Depth values labeled beside the samples.
https://drive.google.com/file/d/12EThqk7U3W0l-RqqEr4OMtPhK2n2l3aJ/view?usp=sharing
AS long as I had decent size ground truth data sets then its time to simulate the bathymetric model. Empirical model works on the principle that each band has a different absorption level over water and this diversity level theoretically will produce the ratio between band. (Stumpf et al. 2003). With this theory depth can be defined using the following formula.
𝑍= 𝑚1 [ln(𝑛𝑅𝑤(𝜆𝑖))/ln(𝑛𝑅𝑤(𝜆𝑗))] + 𝑚0
where n is a user-set constant which should be chosen to ensure the logarithms are positive and the relationship is maximally linear.
The Rw terms are the atmospherically corrected reflectance's in two bands i and j.
m1 and m0 are constant and calculated by calibrating dataset of points with known depths .
The ratio is applied for the pair of Blue-Green (B2-B3) bands and for n = 1000.
It is important to mention that from GEBCO 2020 Grid we only need to calculate the constant m1 and m0.
We also developed correlative plot using the predicted depth and observed depth.
https://drive.google.com/file/d/1F-bSBqmqOYqX3E_TRzY3x1In7JOku8v1/view?usp=sharing
As shown on the above figure values are more aligned to the X axis or observed depth than the predicted values or y axis. It is also apparent that the R square value is .62 which needed to be improved. It is also visible from the chart that prediction is not following the line after 32 meter since the model works for shallow water bathymetry. The calculated bathymetry using Empirical model should look like the following image.
https://drive.google.com/file/d/1evebNQrYkwnTGjH_l5AklIL-2vzD-uTY/view?usp=sharing
One reason why the ratio model couldn't perform well is because the bottom features are equally visible on both spectral bonds(B2, B3), We just noticed that the blue, greeen, red or even the NIR band can easily distinguish bottom features and water was very clear. Reflectance of NIR and Red band should be fully absorbed on certain water depth however due to shallow condition the bottom is becoming visible. The issue also affected while applying sun glint corrections, subsurface NIR reflectance overcorrect the shallow water areas and may appeared dark. Also observed that there are less variability on the bottom type. The entire bottom had unique type of bright sand. Another reason for low performance is because there are outliers in the calibrated database. A visual inspection was made on the 250 GEBCO sample depths at least 50 outliers found in the database and removed. Some of the outliers could easily be detected as they are having high depth value located close to the shoreline. Similar others extreme values can also be identified by making comparison with RGB composition. As discussed earlier subsurface features are clearly visible and can easily be predicted the high depth and low depth areas therefore any observations doesn't comply with such argument was removed form sample depths. Another way to detect the abnormalities on the input data is observing homogeneity of the depth values. In QGIS depth values were plotted and labeled, with this process regional outliers were easily visible.
Step 07. Developing Linear Model
To get advantage of such environmental settings, We understand that it is more logical to apply linear model which could address the issue more efficiently. Next we developed our own version of bathymetric model using single band (B2). A linear regression is implemented between the reflectance of blue band and corrected GEBCO sample heights databases. This time we found strong correlation with the with the brightness of the blue band and the observed depths. https://drive.google.com/file/d/1J7szMNYFKLTogW0msF7izrksu-_DzvGd/view?usp=sharing
Initially R square value was .74, after applying a Log transformation on the sampled depth R square value increased to .78.
https://drive.google.com/file/d/1Dke97sD7gRu-RB4FEs4blcIGc3tFNyG_/view?usp=sharing
Using the slope and intercept from this regression line, brightness value from blue band is scaled ocean depth. With the linear transformation we produced our primary version of bathymetrioc product and finally a low pass filter with 3x3 AM applied on the bathymetric product and resulted the R square value was .78835. Following image shows GEBCO 220 Grid on the left and Sentinel -2 satellite derived bathymetry on the right.
https://drive.google.com/file/d/1gdWzh8yObONXl2gE2K4fVsmqKY3_Sncy/view?usp=sharing