Introduction

For the problems considering water level changes be it flooding or shallowing it is required to have a combined elevation model of underwater and ground terrain. There exists now a special term topo-bathymetry for this kind of models. And the most critical part of the model, which should be most precise here, is its nearshore areas of the elevations in the range of possible water level changes.

While the ground terrain is known well enough, the bathymetry is usually very scarce. When implementing a project on modeling flood zones of the Angara River we have analyzed the available sources of information on bathymetry and came to the conclusion that among all the available data sources on the underwater terrain the pilot charts are the best we have, in spite of that they are very schematic. Therefore we have developed a technique for morphing maps to be able to combine the less precise schematic data with the more precise topographic maps [1]. Further comparison of the data obtained this way with the results of manual depth measurements from boat using echo sounder in some areas of the Bratsk reservoir have shown that the actual measurements can differ significantly from the a-priori data obtained from the pilot charts [2]. Thus, in order to refine the model of the underwater terrain, it is required to find additional sources of information about the depths. Meanwhile, the depth measuring using echo sounder is a very time consuming and expensive task; therefore it will not be possible to process the entire water area of a large water body this way. The multi-beam echo sounders [3] make it possible to capture much wider bottom areas under the vessel, but this kind of devices cannot work effectively in the areas of shallow water, which we consider here.

Another very promising technology is the topo-bathymetric LIDARs [4,5]. The topobathymetric LIDAR uses infra-red laser beams, which doesn't penetrate water, to measure position of water surface and green laser beams, which pass through water, to measure the distances to the bottom points up to the 1.5-3×Secchi depth. The serious devices of this kind require to use airplane and their prices (which are always available only on demand) are prohibitively high (they are comparable to that of the airplane). At the very least, the representatives of Leica Geosystems in Russia don't even try to sell this kind of devices, and we were yet unable to find any company, that can perform this kind of field work.

Satellite remote sensing allows us to obtain the largest possible coverage of the observed territory. The mapping from such satellites as Sentinel and LandSat is carried out regularly, and the resulting data is freely available at the resolution of 10 m per pixel for the primary channels [6]. There exist methods of using remote sensing data to obtain information about the underwater terrain to a certain depth. In [7] an overview of various approaches to solving this problem is given. Among them: the use of stereoscopic satellite images of high accuracy, processing of satellite radar images of various types, inversion of bathymetry from data on water color, restoration of bathymetry from water wave characteristics. Most of these methods require to use the remote sensing data of the kinds, which cannot be obtained free of charge. Among the methods considered here of greatest interest to us are the methods of inverting bathymetry from water color data [8], since they can be applied to freely available data. But for the use of the inverting bathymetry methods extensive calibration and analysis of the accuracy of their results is required, and here the approach proposed in this article may be applied to obtain some elevation data for comparison.

For several years, the Angara-Baikal basin was in the state of shortage of water, which ended in 2019. During this period the water level of Lake Baikal and the Irkutsk reservoir has strongly changed, so that it is possible to find the lines of the coastal terrain level using satellite data, simply by tracking the changes in the actual position of the coastline. Despite the fact that the length of Lake Baikal is about 600 km and tides and seiches (standing waves) appear on it, these phenomena can be safely neglected in a first approximation, since the maximum (Sun + Moon) tide changes the water level by only 32 mm, and seiches by no more than 60 mm [9,10]. The most prominent phenomenon here, which may distort the results, is the water waves. Another object in the Irkutsk region suitable for testing the technology under consideration is the Bratsk reservoir, where the water level changes by up to 4 m per year and the total multiyear range is about 6 m. All the objects under consideration are of very large area; therefore the use of any paid data sources would immediately significantly increase the budget of any project. Therefore, the main goal of our study was to find out what best results can be obtained using the freely available remote sensing data.

The remote sensing images are widely used for tracking changes of shorelines. In [11] LandSat images of Lake Sevan were used for long term (between 1973 and 2015) tracking of the lake shoreline changes. And the article mentions several other research works in the field of water detection including the work [12] where the authors perform a herculean task of tracking changes in global surface water over 32 years at 30 meter/pixel resolution using three million Landsat satellite images and specially designed expert system for pixel classification. For the change assessment of already represented by vector contours shorelines the United States Geological Survey (USGS) has developed a special DSAS tool [13], which is an ArcGIS extension.

The main problem with the open remote sensing data for us is that they have rather low resolution of at most 10 m per pixel, which may be not enough to detect the changes in water level in the range of less than 1 m. And no above mentioned research work considers this problem.

Detection of shorelines

Let us consider the main idea of the approach that makes it possible to use the open satellite images of medium resolution for detection of shorelines with the quality, which is good enough for solving practical tasks.

The traditional approach to shoreline detection

When taking part in implementing a project, which required a combined elevation model of Angara river bed (topo-bathymetry), we have already tried to use shorelines detected on the open satellite data of medium resolution (10 m/pixel) to improve the elevation model. To do it we applied a traditionally recommended approach to finding the coast lines which consists of the following steps:

1) compute the NDWI (Normalized Difference Water Index) raster [14] (for Sentinel-2 images, we use the channels B03 -Green and B08 -Visible and Near Infrared).

NDWI = B03 − B08 B03 + B08;

2) perform binarization of the NDWI raster by some threshold (thresholding); 3) vectorize the binary raster, which leads to obtaining the zigzag lines passing along the boundaries of rectangular pixels; 4) smooth the lines found using some smoothing threshold.

After analyzing the results of the traditional approach we had decided that the resolution of 10 m/pixel of the images is not enough for the problem under consideration and had refused to use these data.

Detection of shorelines with subpixel precision

Now, a few years later we decided to return to the problem of shoreline detection on the open satellite images by taking into account that the binarization step destroys substantial part of the image information. The main idea of the new approach is that instead of vectorization of binary images we can immediately build contour lines of the grayscale NDWI image. The steps of the new algorithm: 2) build contour lines of the NDWI raster with the step of 0.05 between 0.0 and 0.25; 3) determine the best contour line level (usually 0.1).

The Figure 1 allows us to compare the results of the two approaches. In spite of that the smoothed border lines of binary image may look good enough, note that, when the values of the NDWI raster pixels will gradually change, the smoothed binary border lines will remain the same, until a value of some of the pixels will cross the binarization threshold, while the contour lines will immediately reflect the changes. So, the smoothed lines will change in jerks (in a discrete way), and the contour lines will change continuously.

Finding the best contour line level

Next, among the levels of the NDWI contour lines we should find the best one, for which the corresponding contours more precisely match the shorelines. To estimate the contour line quality we compare them to some object, which has a clear boundary with the water in any view direction. A good example of this kind of objects is bridge. And for the bridges (and any other kind of water surrounded objects of stable position and size) we can also use their width to check which contour line is the best: we select the level for which the distance between the contour lines on both sides of the bridge is closer to the bridge width obtained from the corresponding very high resolution image. The Figure 2 shows an example of the use of this approach.

Horizontal alignment of images

After detection of shorelines on several satellite images of the same area we can compare the contours and estimate the resulting quality by their conformity.

The need for shift

As a result we can usually see that the subpixel precision of the contour lines reveals that the accuracy of horizontal positioning of the images is limited by their pixel size. Even for the images corresponding to different water levels we can reveal some inconsistencies, like incorrect order of some contour lines (when a shoreline for the lower water level goes deeper into the ground, than that for the higher water level). I. e. we need to shift one of the two layers by each coordinate in the range of ±pixel size to align them correctly. Figure 3 shows an example of the mismatch and a shift performed manually, which improves the contour matching for the shorelines obtained for the low and the high water levels.

Automatic alignment of contours

Trying to find the best shift for alignment of two vector layers manually is very tedious and timeconsuming task. So we have developed a simple algorithm for finding the shift automatically, which selects the two largest similar contours in the area of intersection of the layers and several times moves the first layer in the direction of the average shift from the points of the contour being moved to their projections on the second contour. We have described the algorithm in more details in [15]. In Figure 4 we can see the original positions of the detected shorelines and the results of the shift performed by the algorithm. Here we use two intersecting images for different but very near (2 days apart) time moments with very close water levels, and the images are even in different projections. And anyway after the shift the detected shorelines become very similar. It means that we have found almost the same shorelines on both images in spite of the differences of their registration details -it makes us more confident about the suggested approach fidelity.

The water bodies of interest

During the last decade here in Irkutsk region we have witnessed rather unusual situation, when the state of water shortage have immediately changed in 2019 by the state of water abundance, which continues thus far. So we are able to find satellite images corresponding to wide range of water levels during relatively small time period. In spite of all the public concern about the level of Baikal, the overall range of the Baikal level changes was rather small, about 1m. So we decided to consider an object with much more prominent water level changes -the Bratsk reservoir. Table 1 contains the ranges of water level change for the Bratsk reservoir by year from the site [16]. From the table we can see that the largest difference between the lowest and the highest levels was observed in 2019.

In fact we will be able to detect the shorelines only when we can see the water surface. And the water will be covered by ice and snow in the late fall, the winter and almost all the spring. So, less than half of the images are suitable for processing. And for this suitable part of the year the main problem is clouds.

Experimental results

Processing images for a tile

After successfully processing several arbitrarily selected Sentinel-2 images, as it was described in the Sections 2 and 3 we decided to select some tile (or granule in term of Sentinel documentation) in the area of the Bratsk reservoir (we use the tile T48VUH) and find out how many images we have in 2019, and process automatically them all. We use only the images of good enough quality with no more than 30% of cloud coverage according to the image attributes (we'll discuss this cloud coverage estimate in more details later). For the tile T48VUH we have 63 images in 2019. Figure 5, a presents all the contour lines for the images. In all the images the colors of the layers represent the corresponding to the image date water level. When considering a particular area it should be noted, that not all the images cover the whole tile: some of them may cover only the smaller left part, and some -the larger right part (it depends on the satellite track (or data strip) of the image). We will consider in more detail a smaller area in the left part (Figure 6), and only 32 images intersect this area (Figure 5,b). And among the images for the left part only 21 of them belong to the "summer" period (Figure 5, c).

The Figure 6 shows a fragment of the images from Figure 5, a. Figure 6, a contains all the contour lines from the 32 images we have for the area in 2019. From Figure 6, a we can see, that most part of the lines in the ground area, which were observed in winter, should be explained by snow areas and not by clouds, because the similar structures were detected in several images. The Figure 6, b represents only the "summer" contour lines (21 layers) -we have found the corresponding time interval manually by looking for the first spring and the last autumn image with the contour lines of good quality. And all the remaining wrong contours in the summer images are explained by clouds, because we can't see any similar objects here. After manually removing the summer layers with the bad contours in the area of interest we have 15 good enough layers, which are shown in Figure 6, c. In fact some of the 15 images may contain bad contours outside the area of interest. And from Figures 6, c and d we can see that the order of the contours is generally correct -the smaller the contour, the lower (more green) its level.

In Figure 6, d we can see the results of alignment of the contours. We perform the alignment by executing the following steps:

1) cut from all the layers the fragments in the area of interest; 2) sort the layers by the corresponding water level in descending order; 3) for each layer starting from the second one align the layer to the previous layer (already aligned) using the algorithm described in the Section 3.2.

Figure 7 allows us to compare the detected contours to the actual shorelines visible on the corresponding images of higher resolution. Here we can see that all the noticeable bends of the shorelines were captured.

Handling clouds

The necessity to manually analyze image quality to detect snow and clouds makes impossible completely automatic processing of the available images. So we have made some attempts to apply algorithms, which detect clouds and snow. We'd even prefer to overestimate the cloud areas to obtain only real contours. Unfortunately the quality of the results of the algorithms was not enough for our purposes. Figure 8 summarizes some our attempts using an example of the image for the already mentioned in the previous chapter Sentinel tile T48VUH for the datetime 20190711T042711. Figure 7, a shows the cloud mask from the file MSK_CLOUDS_B00.gml supplied with the Level 1C images. Here we can see that the mask looks ridiculously small -it have missed almost all the clouds in the area. It even makes us doubt the accuracy of the cloud percent estimates reported for the Sentinel images. The Figure 7, b represents the results obtained by GRASS GIS using the i.sentinel add-on, which is specially designed for processing Sentinel-2 images. Here we first perform atmospheric correction using the i.sentinel.preproc command and then run the i.sentinel.mask command, which detects not only clouds (orange) but also their shadows (yellow). Unfortunately in spite of more than one hour spent on the processing of the single image the GRASS GIS add-on also misses substantial part of clouds, and some of the remaining clouds are located in our area of interest -the shore line. So any attempt to use the image after applying this mask would disrupt the resulting contours.

One of the best packages for cloud detection according to comments of colleagues and reviews from articles is the sen2cor utility [17], which computes Level 2A Sentinel product from Level 1C images. The utility builds a cloud probability mask (in the file MSK_CLDPRB_20m.jp2) and categorizes the image pixels into several classes (water, forest, clouds, cloud shadows, dark areas, and so on). The cloud probability mask is shown in Figure 8, c over the corresponding true color image -here only the zero pixels of the mask are transparent (and the mask pixels are the 8-bit unsigned integer values). The classification result here is indeed much better: almost all the clouds were found. But a closer look at the results reveals that some thin clouds over water are still not detected (Figure 9, a) and, which is much more important for us, that the pixels on the border between water and land are often classified as clouds (Figure 9, a and b). Some pixels of the false clouds on the shore line may have rather high "probability" values of more than 40, and some rather big cloud areas contain only the pixels of low probability (less than 30) so it is not enough to use a threshold of about 30 to filter out the wrong results. And the examination of the results of the pixel classification (Figure 8, d) shows, that the most interesting for us coast areas are often classified as just dark pixels. So by some unfortunate turn of events the quality of the Sen2Cor results in the areas of shorelines prevents us from using them without any sophisticated post-processing.

Conclusion

We have demonstrated that the contour lines of NDWI raster images computed for the Sentinel-2 remote sensing data allow us to detect the shore lines with the subpixel precision. Using the example of a fragment of the Bratsk reservoir, we have shown that it is enough to process the satellite images obtained in one year to obtain the good enough information about the nearshore terrain. We have shown the need to align the contours obtained from several images and have described a way to do it for all the series of available images.

Unfortunately we can't make the process completely automatic, because we have not found a satisfactory solution for the more general problem of cloud detection. And it is still required to develop a technique to combine the contour lines for various water levels into a terrain model, which will take into account the possible errors in the contour lines and will allow us to neutralize the errors by collecting more data.