Bathymetric maps from multi-temporal analysis of Sentinel-2 data: the case study of Limassol, Cyprus

During the last years, various methods and techniques have been used to estimate the seafloor topography with high accuracy. The topography is considered an important component for marine applications since several activities and infrastructures are being carried out at the bottom, column and surface of the sea. In this, study, freely and open distributed optical satellite images were used to obtain bathymetric data for shallow waters, (until 30 m below sea level) based on timeseries analysis of multispectral Sentinel-2 datasets. The ratio transform algorithm was implemented for twelve (12) monthly images covering thus a whole year. The different generated Digital Terrain Models (DTMs) were compared with high resolution LIDAR measurements. The results showed that bathymetry can be obtained from satellite data within a Root Mean Square Error (RMSE) ranging from 1.39 up to 2.56 m while more accurate results were generated during the summer period.

1 Introduction

Bathymetrical data are a crucial component of the marine area since several activities and infrastructures are being carried out at bottom, column of the sea and sea surface. The importance of bathymetry is quite high in areas where there are shallow waters, such as harbors. Nowadays, various remote sensing datasets are used to determine the topography of the seafloor with high accuracy, including active sensors. Traditionally, the active sensors used for estimating bathymetry are the echo sounders (Multi-Beam Echo Sounders, Single-Beam Echo Sounders) or LIDAR devices installed on vessels following a course along specially designed sailing lines that have a specific geometry (Kearns and Breman, 2010). Beyond the fact that these methods are very accurate, they are also very time consuming and costly (McCaffrey, 1981; Sánchez-Carnero et al., 2012). Airborne LIDAR is another method that has been implemented for estimating the bathymetry which revealed to be cost effective compared onboard the vessels providing high area coverage faster (Irish and Lillycrop, 1999; Muirhead and Cracknell, 1986; Ogawa et al., 2017).

During the last years new studies related to Satellite-Derived Bathymetry (SDB) emphasize the potential use of optical satellite remote sensing sensors for bathymetric purposes. Indeed, analytical, semi-analytical and empirical approaches may be implemented for the estimation of bathymetry of up to 30 m depth. The analytical and semi-analytical approaches are based on the physics of light's transmission in water using different parameters of the atmosphere, water column and bottom material. By comparing the three methods, the analytical and semi-analytical methods seem to be more complicated than the empirical one, with more parameters in order to extract bathymetrical data (Gao, 2009). Towards the end of 1970s, it was demonstrated that bathymetry could be estimated over clear shallow water using remote sensing with a multi-band log linear algorithm (Lyzenga, 1978, 1980). Afterwards, this analytical approach has been determined in various ways or with small modifications to derive bathymetry (Conger et al., 2006; Lyzenga et al., 2006; Maritorena et al., 1994; Philpot, 1989). An algorithm for determining the depth of penetration (DOP) zones for every band and calibrating depths presented by Jupp (1988). Also, semi-analytical approaches were implemented, using models and methods for retrieving bathymetry data (Adler-Golden et al., 2005; Lee et al., 1999). Systematic errors are experienced by using the analytical and semi-analytical approaches on the coastal areas containing seagrass because of the reflectance which is lower on the neighboring deep waters. This issue derives most probably from the shading which occurs within the canopy (Dekker et al., 2007). Another method was implemented, named ratio method, for the extraction of bathymetric data over clear shallow water (Stumpf et al., 2003), in which the authors argue that the ratio for two bands at the constant depth will be the same, independently of the difference of bottom albedo, and can be scaled with actual depths.

Over the years, many researches used several optical satellite imageries to extract bathymetrical data. As, remote sensing technologies are continuously upgrading, they are also expanding the horizons of researchers in the field of bathymetry using high resolution satellite images (Said et al., 2017; Setiawan et al., 2013). These multispectral sensors, provide a wide range of wavelength bands which can be applied to these models, especially green and blue bands which can penetrate up to 20 m underneath the sea surface in a clear water condition (Sutanto, 1992).

Despite the medium spatial resolution, open and freely distributed data have also attracted the interest of researchers. Indeed, during the years, different studies carried out to estimate bathymetry using freely distributed medium resolution (30 m pixel size) Landsat-7 ETM+ and Landsat-8 LDCM satellites. The most common algorithms used in these studies where the log linear (Lyzenga, 1978, 1980) and the non-linear models (Stumpf et al., 2003). For instance the non-linear algorithm was applied to Landsat images in the area of the Dutch Wadden Sea (Ameland Inlet, Misra et al., 2018) with a coefficient of determination (R²) of 0.85 and Root-Mean-Square Error (RMSE) 3.12 m; on Thousand island (Pushparaj and Hegde, 2017) with R²=0.90 and RMSE = 0.86 m for depths up to 10 m and at the coastal region of Dakshina Kanada in India, Kastela bay in the middle Adriatic with R² of 0.90 (Setiawan et al., 2017).

The exploitation of freely and open distributed datasets for bathymetric purposes has been supported by the new capabilities of the Sentinel-2 sensors, having a 10 m resolution. Recently, the log linear and non-linear methods were used for mapping the bathymetry, on the South Baltic Sea (Chybicki, 2017). In that study, satellite imageries from the Sentinel-2A sensor and soundings data were obtained by a single-beam echo sounder and were used to estimate SDB. Afterwards, error distribution and quality coefficient were calculated to describe the quality of the retrieved SDB. Further, the retrieved SDB was calibrated with the indicator of quality coefficient. In this way, more accurate results were retrieved from SDB, proposing the implementation of this indicator on mapping bathymetry using remote sensing. The best accuracy reported using the ratio transform algorithm method was with an R²=0.91 and RMSE = 1.08 m. At the same time, Kabiri (2017) used linear, non-linear methods for bathymetric purposes on the Nayband bay, in Iran, using two Sentinel-2 satellite imageries of February and March 2017. The best accuracy was estimated to 1.89 and 2.04 m for months February and March respectively, using the ratio transform algorithm.

More, recently the benefits of the Sentinel-2 imageries using big data cloud platforms were presented by Traganos et al. (2018). The Google Earth Engine, a cloud-based geospatial computing platform, was used to acquire Sentinel-2 imageries for the period between 1 August and 31 December 2016, in four different locations at the eastern Mediterranean basin. Analytical and empirical methods were implemented, including the pseudo-invariant features for normalisation and the median 3×3 filter to smooth the remaining noise, hence increasing SDB accuracies. The best accuracy spotted using the analytical method had a R² of 0.90 and RMSE = 1.67 m.

Though these studies have shown the great potentials of estimating bathymetry using freely and open source satellite data such as Sentinel-2 images, there is still a research to be performed regarding the identification of the best period throughout the year where bathymetric errors are minimized. In addition, bathymetric studies based on optical sensors are rare in the wider Eastern Mediterranean region and therefore extensive research will benefit in the future. Even though satellite remote sensing in this area was used in coastal areas for several purposes such as in marine spatial planning, monitoring archaeological sites, quality of coastal water etc. (Agapiou et al., 2011; Hadjimitsis et al., 2010, 2013; Papoutsa et al., 2011) no bathymetrical studies have been made. Currently through an existing funded project named as “THAL-CHOR 2” (“ΘAΛ-XΩP 2” in Greek), bathymetric maps are being developed through the use of the remote sensing which will be used for the implementation of marine spatial planning in Cyprus.

For the framework of this study, optical Sentinel-2 satellite data were used to extract bathymetric information in Cyprus. The present study is comprised of five parts. Further to the introduction, the datasets and methodology part follows, mentioning the field data of Multi-Beam Echo Sounders and LIDAR instruments. Then, the empirical approach and the methodology used in this study is described. The next part concentrates the results of the applied methodology providing the correlation between observed depths and SDB depths, twelve bathymetric maps. Finally, the paper ends with a discussion and the conclusions.

2 Case study and datasets

2.1 Study site

The study site is located south west of the city of Limassol in Cyprus and it covers 19 km² with a coastline that is about 9.6 km in length (Fig. 1). The site is located between 34^∘41^′38.4${}^{\prime \prime }$–34^∘42^′7.2${}^{\prime \prime }$ S and 33^∘4^′30${}^{\prime \prime }$–33^∘10^′44.4${}^{\prime \prime }$ E (WGS 84, 36N). On the coastlines of the study area, the soft substrate with sandy beaches from shallow waters to a depth of twenty meters is gradually becoming muddy (DoF, 2012). Also, breakwaters can be found on the coasts of Limassol at a distance of 120 m from the shoreline. Cyprus has an intense Mediterranean climate with hot dry summers, rainy, rather changeable, winters and in autumn and spring season the weather conditions are changing rapidly. The winds usually are light to moderate and rarely unstable.

Figure 1 Location of study site in Cyprus.

2.2 Sentinel-2 optical satellite data

In this study, Sentinel-2 imagery and field data were used to analyse and estimate bathymetry. Sentinel-2 mission is part of the European Space Agency (ESA) with high-resolution with multi-spectral imagery. The current research aims to compare the estimated results of bathymetry for each month separately. Therefore twelve processed satellite images (one for each month) were acquired through the Copernicus minor site (ESA, 2013), covering a period from April 2017 to March 2018. Table 1 shows the details of acquired satellite image that were used. The acquired packages are of Level 2A i.e. orthoimages and bottom of atmosphere corrected reflectance product. At this stage, the product was processed on radiometric, geometric, cloud screening and atmospheric corrections. These multispectral images have 12 bands of spectral range at three different resolutions (bands/resolution: b1 – Coastal aerosol/60 m , b2 – Blue/10 m, b3 – Green/10 m, b4 – Red/10 m, b5 – Vegetation Red Edge/20 m, b6 – Vegetation Red Edge/20 m, b7 – Vegetation Red Edge/20 m, b8 – NIR/10 m, b8A – Vegetation Red Edge/20 m, b9 – Water vapour/60 m, b10 – SWIR/60 m, b11 – SWIR/60 m, b12 – SWIR/60 m).

Table 1 Details for the processed Sentinel-2 imageries.

2.3 LIDAR in-situ data

The in-situ data were collected and analyzed by the Department of Land Surveyors of Cyprus (DLS, 2014) – the official chartographic provider of the country – covering water depths up to 50 m for an area of approximately 508 km². The data were acquired in two phases; the first phase was conducted in 2012 when the DLS launched process mapping of the sea with multi-beam echo sounders using a conventional vessel covering 238 km²; the second phase was conducted in 2014 using the airborne LIDAR method covering 270 km². The vertical accuracy as estimated by the DLS is ±0.5 m for depths 0 to 20 m and ±1 m for up to 50 m. Consequently, for this research, the DLS provided a set of 244 random control points of 0 to 22 m depth which have been used for: (a) calibrating the bathymetric results obtain from the multispectral imageries to absolute depths and (b) for training the ratio transform algorithm. Figure 2, shows the overall geographical distribution of the control points within the area of interest.

Figure 2 Field data acquired by the Department of Land Survey in Cyprus.

3 Methodology

Several methods have been implemented for the estimation of the bathymetry using remote sensing. This research, gives emphasis in one of the most commonly used methods (ratio transform algorithm) for deriving depths using satellite imagery, (Stumpf et al., 2003). In this method, the approach of Beer's Law was accepted which considers that light attenuation in the water column increases exponentially as depth increases (Gao, 2009). The log transformation algorithm is related to the natural logarithm for reflectance of two bands and actual depths.

Therefore, the bands used for estimating the bathymetry are the blue (440–540 nm), green (500–600 nm), red (600–700 nm) and near infrared (700–800 nm) bands. The penetration in the blue and green part of the spectrum is higher while the absorption of the electromagnetic radiation is increasing as we move to the red part of the spectrum. The ratio transform algorithm can be apply with bands having different water absorption and can be applied in appropriate wavelengths of any sensor. Since blue and green bands have lower absorption, their ratio should remain the same despite the different bottom albedos at a constant depth. The values of these bands were implemented on the following Eq. (1) to estimate bathymetry.

$$\begin{array}{}\text{(1)}& Z=m\mathrm{1}{\displaystyle \frac{\ln R\left(\mathit{\lambda }b\right)}{\ln R\left(\mathit{\lambda }g\right)}}-m\mathrm{0}\end{array}$$

where Z is the satellite derived bathymetry depth, m0 is the offset for a depth of 0 m, m1 is a coefficient to tune the model to the actual depth, and R(λb) and R(λj) are the remote sensing radiances for optical bands λb (blue band) and λg (green). For this model, the bottom depth is estimated on the basis of light attenuation phenomena, as the attenuation of the incoming shortwave radiation varies spectrally.

Figure 3 Workflow processing steps detecting bathymetry from Sentinel-2 (L2A) satellites images.

The ratio transform algorithm was implemented to the multispectral imageries listed in Table 1. Initially the multispectral images were resampled by a referenced band with a 10 m spatial resolution, using the SNAP software. Upon the selection of the area of interest, the ArcGIS software was used where the shapefile was created and uploaded on the SNAP software to mask satellite imageries. Following, each one of the masked imageries classified the blue, green, red and near infrared band and then separated land and water. After the pre-processing of satellite data, the processed product was uploaded on the ENVI 5.3 software for estimating the bathymetry using the log transformation algorithm. The median filter, a nonlinear digital filtering technique, used for removing the noise from images, was applied to this study. A median filter of 7×7 was used to succeed the reduction of the noise of high frequency. The following figure (Fig. 3), shows the processing steps for deriving bathymetry on the study area.

4 Results

Following the above methodology, twelve bathymetry maps, one for each month, were estimated. The median filter has been implemented to eliminate parameters of deterioration in estimating bathymetry such as turbidity, waves, sun glint and sky glint. Following a time series regression analyses, the RMSE and Mean Absolute Errors (MAE) have been calculated.

4.1 Time series regression analyses of Bathymetrical data

In order to estimate the correspondence of the satellite bathymetric model with the Lidar datasets a regression analysis has been performed for each monthly DTM using the linear model. The overall results are found in Fig. 4. Twelve scatterplots have been designed with the x-axis representing the LIDAR points, while the y-axis represented the estimated bathymetry points derived from the Sentinel-2 images. Based on the linear regression analysis the R² coefficient was also estimated for each case. According to the results, strong correlation is observed in shallow waters (less than 12 m) compared to deeper areas (more than 12 m) . Scatterplots of January and February show that R² is relatively high up to 0.90 while R² for the next couple of months is estimated to be between 0.85 and 0.88. During summer season R² value is increased while June provides the best linear fit with a value of 0.95. Lastly, the R² value for autumn, seems to have a downward trend and is the same with winter's one. The scattering starts from 12 m depth, as also clearly showed for the analysis during November. This is reflected to the R² value which is equal to the 0.81 and it is the less precise value from the results.

Figure 4 Correlation between observed depths and SDB depths obtained for each month.

In order to compare the accuracy of the time series extracted results, two statistical parameters were taken into consideration along with correlation coefficient R², RMSE and MAE. The equations for the estimation are given in the following equations:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{J=\mathrm{1}}^{\mathrm{1}}\mathrm{|}{Z}_{\mathrm{SDB}}-Z_{\mathrm{FB}}\mathrm{|}\text{(3)}& {\displaystyle }& {\displaystyle }\mathrm{MAE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{J=\mathrm{1}}^{\mathrm{1}}(Z_{\mathrm{SDB}}-Z_{\mathrm{FB}}{)}^{\mathrm{2}}}\end{array}$$

where n is the number the field points, Z_SDB is the satellite derived bathymetry depth and Z_FB is the field points depth.

The overall results regarding the RMSE and MAE errors are shown on Table 2. As resulted from the table below, the best results were the ones derived for month June where the bathymetry resulted with RMSE to be 1.39 m and MAE 1.05 m. On the other hand, the less precise results having a RMSE of 2.56 m and a MAE of 1.64 m seem to be for month November.

Table 2 R², RMSE and MAE of the estimated depths.

4.2 Results of bathymetrical data

This section demonstrates and compares the various SDB maps produced from the Sentinel-2 images with the LIDAR datasets. Figure 5 below, shows the bathymetrical results of the study area for each month. Overall, the best accuracy is evident during the summer period with RMSE of 1.39 m and MAE of 1.05 m while the autumn period gives the lowest precision with RMSE of 2.56 m and MAE of 1.64 m. The latest may be due to bad weather conditions and the sea state. It is worth noted that bathymetric maps obtained during the period November December encountered error in the South-East. These errors are probably related to cloud coverage of the satellite images.

Figure 5 Predicting and mapping SDB with ranges 0 to −22 m depth, for the months (a) January, (b) February, (c) March, (d) April, (e) May, (f) June, (g) July, (h) August, (i) September, (j) October, (k) November , (l) December.

As it can be seen from the presentation of the results, the period used for estimating bathymetry effects the overall performance of the results. The records provided by the DLS have a 5-year difference from the downloading of the satellite images which is considered to be a very long timeframe, thus data have limitations on the accuracy. In order to compare the estimated bathymetry with the in situ data, comparison maps between the bathymetric satellite DTMs and LIDAR DTMs were created and can been seen on Fig. 6.

Figure 6 Comparison maps of the estimated bathymetry and in situ data for the months: (a) January, (b) February, (c) March, (d) April, (e) May, (f) June, (g) July, (h) August, (i) September, (j) October, (k) November , (l) December.

Based on these findings it is evident that during winter and autumn periods errors can be more than 5 m. Higher errors of more than 4 m are reported in deeper waters for all seasons with the exception of the summer season. In shallow waters the same performance despite the period selected was observed. In addition, in sea areas with breakwaters the differences range from 1 up to 2.5 m.

Bathymetrical profiles were also created in areas that showed frequently encountered errors. Optically observing the maps with the errors' allocation, it can be pointed out that there are more frequent errors in the southwest of the study area. Bathymetrical profile lines of each month were used in order to interpolate the bathymetrical data. Besides, another profile was created, which is based on the LIDAR control points in order to be compared with those of SDB. The above described method was followed with the use of ArcGIS v10.2 Software for the analysis of the diagrams which have been exported by the software, the user has resulted to the diagram depicted in Fig. 7

Figure 7 Sector profiles of SDB.

A very good correlation between the profiles is perceptible to Fig. 7. For a horizontal distance from the coastline up to 200 m with depths of 5 m, a noise was observed between SDB and the interpolated bathymetry from the field points. This might be due to the wave collapse on the shoreline and the breakwaters which are located 120 m from the shoreline. On the other hand, on the coastal distance of 200 to 400 m in the area with depth of 5–10 m seems to be well correlated except the month of January. Moreover, the correlation of the bathymetric data with distance ranging from 400–800 m seems to be good except from the months of March, April and December. Lastly, the distance ranging from 800–1200 m with 20 m depth has a gradual upward shift of from the fields' points' surface. This might be due to the attenuation of the light as it passes through the water column and the deterioration of solar radiation at various depths in the water. Another reason could be the time difference between satellite images and field points.

5 Discussion

On the current study, twelve images, one for each month, have been used for estimating the bathymetry in the south area of Cyprus. During the analysis, the overall results indicate an error ranging from 1.39 up to 2.56 m. Comparing the bathymetric data, it is evident, summer months provide the best results. Indeed, the best accuracy of bathymetry data is found in Sentinel image acquired on June with a R²=0.94, RMSE = 1.39 m and MAE = 1.05 m. Also, a good correlation between the bathymetric maps of May's to Augusts' (Fig. 5) with accuracies of the RMSE ranging from 1.39 to 1.70 m was observed. In contrast errors of more than 5 m during the autumn and winter months, especially in October and November have been reported in deeper waters. The accuracy of the study is considered acceptable, compared to other studies using the same method and satellite imageries such as those of Chybicki (2017), Kabiri (2017) and Traganos and Reinartz (2017). The correlation between situ data and predicted depths, is not the ideal and this may be due the limitations of accuracy that have LIDAR data. Also, a gradual increase in the accuracy of the results was noted in the months of June and December, months of the winter and summer solstice.

An overall analysis indicates high accuracy of the estimated depths in shallow waters which then reduces in deeper waters. Studying the bathymetrical profiles shows a good correlation between the observed and extracted depths from 5 to 12 m while there were deviations in deep waters. The correlation between the in situ data and results can be found in other studies that show better results in shallow waters (Ehses and Rooney, 2015; Pushparaj and Hegde, 2017).

Further to this study, comparison maps were created. The outcome of the maps shows the difference between the estimated bathymetry with the in-situ data. These comparison maps show good correlation to the shallow waters, with the differences ranging from 0.5 to 1.5 m. It is also worth to notice that breakwaters in the study area might have deteriorated the results as shown in the comparison maps.

6 Conclusions

Bathymetry is an important factor for several applications such as navigation, marine environmental protection, coastal engineering etc. During the last years many studies implemented SDB using empirical, analytical and semi-analytical approaches, aiming to investigate the potential use of optical satellite datasets. A critical aspect highlighted on most of these studies refers to the atmospheric and sea conditions of the area including water column and bottom material characteristics. It is therefore important to analyse an optical satellite image with the best possible conditions (cloud coverage, waves, turbidity etc.). For this reason, the study was focused to the identification of the best period for the estimation of the bathymetry following a multi-temporal approach.

Therefore, the study was focused on the analysis of the multi-temporal Sentinel-2 satellite imageries in different periods using the empirical method ratio transform algorithm. The bathymetrical data were extracted using SDB and statistically analyzed for one year. The overall results show that the bathymetry accuracy tends to fluctuate from period to period significantly and therefore this should be taken into consideration in combination with the weather conditions. The latest is due to the fact that turbidity and waves affect the attenuation of the light observed during bad weather conditions and strong winds.

For the purpose of improving the current study, a comparison timeseries analysis exploiting ESA's Sentinel Scientific Data hub will be performed in the near future. Based on the given high temporal resolution of Sentinel-L2A-B sensors, a multitemporal analysis during the summer season (optimum period for bathymetric purposes as evidence in this paper) will be investigated by authors in the near future. The purpose of that study is to identify fluctuations within a specific season (i.e. summer), following other methods as well, aiming to minimize even more the errors and produce more accurate bathymetric maps. In addition, this will allow isolate the optimum window for exploiting the free and open satellite data.