Mesoscale Ocean Feature Identiﬁcation in the North Aegean Sea with the Use of Sentinel-3 Data

: The identiﬁcation of oceanographic circulation related features is a valuable tool for environmental and ﬁshery management authorities, commercial use and institutional research. Remote sensing techniques are suitable for detection, as in situ measurements are prohibitively costly, spatially sparse and infrequent. Still, these imagery applications require a certain level of technical and theoretical skill making them practically unreachable to the immediate beneﬁciaries. In this paper a new geospatial web service is proposed for providing daily data on mesoscale oceanic feature identiﬁcation in the North Aegean Sea, produced by Sentinel-3 SLSTR Sea Surface Temperature (SST) imagery, to end users. The service encompasses an automated process for: raw data acquisition, interpolation, oceanic feature extraction and publishing through a webGIS application. Level-2 SST data are interpolated through a Co-Kriging algorithm, involving information from short term historical data, in order to retain as much information as possible. A modiﬁed gradient edge detection methodology is then applied to the interpolated products for the mesoscale feature extraction. The resulting datasets are served according to the Open Geospatial Consortium (OGC) standards and are available for visualization, processing and download though a dedicated web portal.


Introduction
Upper layer mesoscale oceanic processes are the primary and most efficient transport of energy across sea regions. Mesoscale phenomena are considered eddies, fronts, cyclonic circulations, meandering currents and other oceanic processes with a spatial radius between 10-100 km [1]. They are responsible for the largest heat exchange between basins and are relevant in both mixing and enriching different water masses with nutrients and mineral traces. Similarly, submesoscale circulations are oceanic features with a smaller radius (1-10 km), but they have an important role in mixing the upper stratified surface layers with subjacent ones [2].
Environmental conditions are highly affected by circulation, both by their spatial and temporal variability, influencing the marine life structure and the fishing industry [3,4]. Nutrients are the driving force of primary production in the ocean and their distribution is mainly affected by circulation. Specifically, mesoscale and submesoscale processes as eddies and upwellings can trap or introduce nutrients at specific locations, respectively making them suitable feeding or spawning grounds for fish populations [5]. Research has also shown that various commercial fish species are dependent on their distance from thermal fronts [6]. Finally, currents and oceanic fronts play a vital role in the interaction between the ocean and the atmosphere and are responsible for enabling or restricting heat transfer across sea regions at different latitudes and longitudes, controlling the conditions of the lower parts of the atmosphere and at extension its circulation [7].
Traditional mapping and monitoring of mesoscale and sub-mesoscale flows is performed by lagrangian drifters [8,9]. These devices produce the most accurate results, as they follow each flow directly giving back a very finely detailed track of their trajectory. They are capable of observing both mesoscale and sub-mesoscale features, but they are limited because each drifter can follow one flow at time. Taking into consideration that they can take weeks to months to complete their survey, they are not suitable for mapping large water bodies at spatial and temporal scales required for understanding fisheries and ecosystem dynamics.
Satellite applications on the other hand, are capable for more frequent measurements over larger geographical regions. Altimetry applications target on measuring the Sea Surface Height and in conjunction with in-situ data, they can estimate the Eulerian velocity of the sea surface flows and the general trajectory of the currents [10][11][12]. However, these are limited due to the time to derive the geostrophic flow field. Applications, relying on ocean color techniques, try to extract information on front location by estimating the difference of sea surface basic parameters values in close distances, and identify sharp edges where potential fronts are located. Such methodologies usually do not rely on complimentary data, and their temporal coverage and resolution are dependent on the sensor from which the data were acquired.
There is an extensive literature of different remote sensing approaches about oceanographic fronts and features detection. The gradient edge detection method is a widespread technique, featured with various alternations and optimizations [13,14]. It is an approach that focuses on distinguishing stiff drops and sharp rises of sea surface physical or biochemical parameters in short distances with statistical means. The Canny edge detector is a similar methodology, where edges are extracted from gradient images after filtering for noise removal and the application of a double threshold [15,16]. Cayula and Cornillon (1992) [17] developed a new histogram algorithm, where fronts are detected within a moving window and as thus it is resilient against high value fluctuations and extremes [18]. Other approaches involve statistical modeling, such as the weighted local likelihood algorithm developed by Hopkins et al. (2010) [19], or the logistic regression used to classify fronts from the results of a Canny edge detector [20]. Finally, very promising results are produced from various machine learning techniques, which can identify oceanic features with very high accuracy [21,22].
While the spatial information of mesoscale oceanographic processes and circulation is vital for several scientific and industry fields, this information is not yet widely available at fine scales and for enclosed sea regions such as the Aegean Sea. The extraction of such information requires knowledge of remote sensing, image processing and time series analysis. More than that, if that information is to be published, those processes need to be enclosed within a web-based service, which in turn requires further technical expertise.
In this paper, a practical application is presented for mesoscale and sub-mesoscale near-real time mapping for the North Aegean Sea. It refers to a fully developed service that aims to give end-users ready to use data, without the need of theoretical and technical knowledge about satellite data handling and oceanographic processes. The service encompasses all the processing stages, from data acquisition and mesoscale and sub-mesoscale feature extraction, to publishing the results as a geospatial web service. The work aims to enhance the oceanic front detection with modern data usage and practicality through a geospatial web service.
The service aims to benefit directly the fisheries management authorities, the fishing industry as well as the scientific community. Daily provided information on oceanic fronts and (sub) mesoscale circulation can help government officials plan out the fleet distribution and restrict or allow access to certain fishing grounds based on the likelihood to host fish populations and the level of the exploitation it has already sustained. Additionally, having access to maps showing how potent an area is for fishing, based on oceanic circulation, and the distance that a front has from the port, professional fishermen, will be able to make more cost effective voyages, by maximizing their catch and minimizing their expenses in fuel and time. Finally, as an open source service, scientific institutions and researchers will have access to a ready-to-use dataset for analysis and implementation into their own research.

Study Area
In order to illustrate the geospatial web service potentials, a well-defined case study area in Greece is used (Figure 1). The study area is located in the North Aegean Sea. The North Aegean is characterized by lower salinity levels (29 psu) compared to the South Aegean (>38.5) [23] and colder waters. The main contributing factor for the low salinity and temperature is the mixing from the outflow of the Dardanelles Strait, where waters from the Black Sea with an annual variation of salinity of 19-22 psu and a mean temperature of 16 • C penetrate the Aegean [24]. Another significant influencing factor for the decreased salinity is the freshwater input from major rivers in the Thermaikos Gulf.
The low salinity and cold waters from the Black Sea influence the Northern Aegean Basin circulation, and they are the main reason they distinguish the Northern from Southern Aegean basin. The Black Sea waters, after entering the Aegean, they are divided into two main masses, the first flows to the North towards the Sea of Thrace and the second flows southern, towards the Sporades islands and then reaches the Southern Aegean through the Kafireas and Mykonos-Ikaria straits [25]. Two major anti-cyclonic formations are present at Limnos and Samothraki islands and a smaller one western of Thasos Island [26]. Another cyclonic formation can be observed at the Sporades archipelagos, which changes seasonally to an anti-cyclonic one. The flow from this formation continues to the South, along the Eastern Coastline of Evoia island, creating a strong and narrow coastal current. and minimizing their expenses in fuel and time. Finally, as an open source service, scientific institutions and researchers will have access to a ready-to-use dataset for analysis and implementation into their own research.

Study Area
In order to illustrate the geospatial web service potentials, a well-defined case study area in Greece is used (Figure 1). The study area is located in the North Aegean Sea. The North Aegean is characterized by lower salinity levels (29 psu) compared to the South Aegean (>38.5) [23] and colder waters. The main contributing factor for the low salinity and temperature is the mixing from the outflow of the Dardanelles Strait, where waters from the Black Sea with an annual variation of salinity of 19-22 psu and a mean temperature of 16 °C penetrate the Aegean [24]. Another significant influencing factor for the decreased salinity is the freshwater input from major rivers in the Thermaikos Gulf.
The low salinity and cold waters from the Black Sea influence the Northern Aegean Basin circulation, and they are the main reason they distinguish the Northern from Southern Aegean basin. The Black Sea waters, after entering the Aegean, they are divided into two main masses, the first flows to the North towards the Sea of Thrace and the second flows southern, towards the Sporades islands and then reaches the Southern Aegean through the Kafireas and Mykonos-Ikaria straits [25]. Two major anti-cyclonic formations are present at Limnos and Samothraki islands and a smaller one western of Thasos Island [26]. Another cyclonic formation can be observed at the Sporades archipelagos, which changes seasonally to an anti-cyclonic one. The flow from this formation continues to the South, along the Eastern Coastline of Evoia island, creating a strong and narrow coastal current. Figure 1. Study Area, the North Aegean Sea, Greece. All major oceanographic processes are mapped, the Black Sea Water driven circulation, the inflow from Southern Aegean and Gulf of Thermaikos and their interactions (adapted from [26,27]). Figure 1. Study Area, the North Aegean Sea, Greece. All major oceanographic processes are mapped, the Black Sea Water driven circulation, the inflow from Southern Aegean and Gulf of Thermaikos and their interactions (adapted from [26,27]).

Data
The dataset used in this study is derived from the Copernicus Sentinel-3 satellite data-stream. Thermal front detection is based on the Sentinel-3 SLSTR sensor which is equipped with 11 bands; S1-S6 are used for cloud detection and vegetation and ice monitoring at 500 m resolution, S7-S9 are used for Sea Surface Temperature (SST) and Land Surface Temperature (LST) estimations and fire monitoring at 1km resolution and lastly F1 and F2 for active monitoring also at 1km resolution. SST data already processed to Level-2 from the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) at a 1km resolution used for further analysis. The data are not interpolated but they provide cloud and pixel quality masks.
Sentinel-3 SLSTR imagery's spatial and temporal properties where the most suitable choice for this application. The high 1km resolution of the SST data is well suited for the differentiation of small enough features, for both mesoscale and submesoscale features. Furthermore, with a swath of 740 km oblique view and 1400 km in nadir view, it covers the whole of the study area without the need of mosaicking. Finally, on the operational side of the service, the sensor has a temporal coverage over the Aegean of 1-2 days, enabling the servicing of daily products.
Level-2 SST data are acquired from the EUMETSATA Copernicus Online Data Access service. The data are stored by the Standard Archive Format for Europe (SAFE), containing the bands (NetCDF) and the manifest header file. A script is run daily searching for newly archived data based on spatial and temporal criteria. If there are data matching the query, their unique product ID is collected and then a URL request is generated for downloading each product, e.g., https: //coda.eumetsat.int/odata/v1/Products(\T1\textquoterightID\T1\textquoteright)/$value.

Methodology
A fully operational service is presented to produce daily mesoscale and sub mesoscale oceanic feature localization maps. The service has three main sections: (a) obtaining the raw satellite data and pre-processing them in order to be used for the detection algorithm; (b) the detection process that produces the maps for the end users; and (c) the implementation of the geospatial web service, that publishes the end results according to the OGC standards ( Figure 2).

Preprocessing-Interpolation
In an effort to investigate the fluctuations of SST, a number of remotely sensed data were used over a number of timespans for the study area. Primary data include unwanted objects such as cloud coverage and land surface cells, which are not relevant to our research question. Thus, there is the need to remove (a) cloud coverage cells and (b) land surface cells as well as (c) extreme values near the coastal areas. The latter are considered to be data "noise" because of the anthropocentric activities nearby, the river outflows and the possibility of low bathymetry. A 5 km buffer zone around coastlines has been excluded. Cloud and land pixels are removed with a masking approach.
After this cell removal activity, there is a need to evaluate SST for the cell locations over sea. This is a typical spatial interpolation problem where some observations are used in the study area aiming to estimate values of a variable on locations with no information [28]. In this paper, these observations are extracted from the same SST data that need to be interpolated and from values derived from data of the previous two days, represented as (T-1) and (T-2). All possible combinations of these datasets are tested for the interpolation modelling and are evaluated according to their accuracy.
For the estimation of values in locations with no information, the well-used Kriging group of spatial interpolation methods are introduced. Kriging is a broad family of interpolation methods with the advantage that it minimizes the estimation variance and error estimation, using a model of spatial covariance acquired by a semivariogram, so it has homoescedastic characteristics [29]. This is a detailed understanding of the spatial variability which comes quite handy in estimating a variable in location with unknown values. The semivariogram is essentially the mean squared difference between pairs of values as a function of increasing separation distance (lagged distance). More details regarding semivariogram is given by Lamorey and Jacobson (1995) [30]. Two variations of Kriging are tested: Ordinary Kriging and Co-Kriging. For both Kriging and Co-Kriging, a fundamental assumption is that the measurement data are normally distributed across the study area both in terms of statistical distribution and geographical distribution [31]. Flow chart describing the structure of the geospatial web service, from data processing to product publishing, where T is the day of the current date.

Preprocessing-Interpolation
In an effort to investigate the fluctuations of SST, a number of remotely sensed data were used over a number of timespans for the study area. Primary data include unwanted objects such as cloud coverage and land surface cells, which are not relevant to our research question. Thus, there is the need to remove (a) cloud coverage cells and (b) land surface cells as well as (c) extreme values near the coastal areas. The latter are considered to be data "noise" because of the anthropocentric activities nearby, the river outflows and the possibility of low bathymetry. A 5 km buffer zone around coastlines has been excluded. Cloud and land pixels are removed with a masking approach.
After this cell removal activity, there is a need to evaluate SST for the cell locations over sea. This is a typical spatial interpolation problem where some observations are used in the study area aiming to estimate values of a variable on locations with no information [28]. In this paper, these observations are extracted from the same SST data that need to be interpolated and from values derived from data of the previous two days, represented as (T-1) and (T-2). All possible combinations of these datasets are tested for the interpolation modelling and are evaluated according to their accuracy.
For the estimation of values in locations with no information, the well-used Kriging group of spatial interpolation methods are introduced. Kriging is a broad family of interpolation methods with the advantage that it minimizes the estimation variance and error estimation, using a model of spatial covariance acquired by a semivariogram, so it has homoescedastic characteristics [29]. This is a detailed understanding of the spatial variability which comes quite handy in estimating a variable in location with unknown values. The semivariogram is essentially the mean squared difference between pairs of values as a function of increasing separation distance (lagged distance). More details regarding semivariogram is given by Lamorey and Jacobson (1995) [30]. Two variations of Kriging Flow chart describing the structure of the geospatial web service, from data processing to product publishing, where T is the day of the current date.
The first method, Ordinary Kriging, implicitly evaluates the mean in a moving neighborhood. It is a weighted average estimation of values while weights are based on Best Linear Unbiased Estimator (BLUE). The second method, Co-Kriging, is used when estimating a single variable from several other variables. There should be a theoretical relationship between the variables of the model. Thus, the covariates must have a statistically strong relationship as well. This approach asks for the spatial covariance model of each extra variable as well as the cross-covariance model of the variables. A couple of Co-Kriging estimations were created with different variables. The two attempts included the following sets of variables: T~(T-1) and T~(T-1) + (T-2).
The evaluation of the performance of Kriging Interpolation was achieved with a classical Data Mining "Leave Out" approach. This is the split of the initial dataset into two group of cases. The first group includes 500 randomly sampled observations which have been excluded from Kriging estimation. This sample serves as the "evaluation group" and it is used only during accuracy assessment of the results. Rest of the cases (n = [8500, 9000]) are used for Kriging estimation.
After estimating values for all locations of the study area, the n = 500 evaluation sample is used for calculating the deviation between estimated values (result of Kriging) against empirical values (evaluation sample data). A group of accuracy metrics are used to compare the two values for the n = 500 observations. These two measures are based on the squared errors, i.e., mean absolute error in Equation (1) as well as the absolute error i.e., root mean squared error in Equation (2).
where e i is the difference between the estimated SST value and the corresponding evaluation data at grid location i. The two measures are scale-dependent, which do not affect the outcome of the measures as sample size (n = 500) is constant across all estimations over time spans.

Oceanic Front Detection
The circulation formation detection method stands on the premise that different water bodies have different physical attributes. It is based on the thermal front gradient detection method [13], which is an edge detection technique for image analysis. A filter is used to identify areas that there is a sudden change of temperature in neighboring pixels. These areas are usually the meet points of different currents that create an eddy, the fringes of a gyre or a cyclonic formation. This category of thermal front detection methods is prone to noise due to inherited low gradient magnitude that characterizes the oceanic flows, but they are very lightweight procedures, making them excellent for operational and automated services.
For the edge detection, the Sobel operator was used [32]. It is a moving window-based operation that compares each pixel with their immediate neighborhood. Peak values are found at pixels where a sudden change in SST occurs. The operator performs the gradient estimation separately for the x and y axis of the image, by calculating the directional gradient of the SST (Equations (3) and (4), respectively). Those directional gradients are then used to calculate the overall magnitude on the image (Equation (5)). The points where the gradient reaches its peak, are marked as potential edges according to the following equations.
where Im n , the SST image of day n, G xn and G xy the SST gradients for the x and y axis, respectively, and G n the overall gradient magnitude. Not all the regions are oceanic fronts though. Pixels can produce some level of gradient magnitude for various reasons not always associated to oceanographic phenomena, such as: sensor noise, shallow areas near the coastline, river outflows and the edges of clouded areas. For that reason, distinction must be made between the possible oceanic fronts and the false positives pixels. Generally, gradient values produced by noise in the image are on the lower range, and those related to coastal processes and clouds at the top extremes. In order to extract only the useful information concerning the oceanic features, a double threshold is performed according to the interquartile ranges of the image histogram. Ranges between the first and fourth quantile are considered pixels with weak gradient magnitude values and those above the fifth pixels related to noise or extremes. These ranges are automatically redefined for each new image and are adapted accordingly.
The weak gradient magnitude that characterizes oceanic processes also has a different effect on thermal front detection, that it occasionally produces discontinuities. Specifically, pixels at the center of a front may have a weak gradient signature and, as thus, they fail to pass the threshold value and split the front into multiple parts. To overcome that effect, a morphological closing technique is applied to the thinned edges. After the threshold is performed, the resulting image is transformed to a binary one, whether each pixel has a strong magnitude signature or not. Next, the closing takes effect and it connects pixels that are in proximity to each other in order to make more coherent objects.
Even though a general positioning of the oceanic features has been established, the localization is not yet accurate enough. At the meeting point of two currents there is usually a mixing zone and if so three water masses are created with different characteristics. The first and second are the two meeting currents with their own separate physical attributes. The third is the mixing zone around the thermal front with characteristics based on the two mixing flows. As such, edges from thermal fronts are not sufficiently strong and the gradient magnitude produces high values around the edge. To compensate for that effect, a skeletonizing technique is applied to reduce the edges to a single pixel line positioned at the center of the edge. The thinning works by comparing each pixel with its neighbors and decides whether to reduce it to a 0 based on a set of rules [33].
It should be stated that the lack of in-situ data renders the evaluation and the accuracy assessment of the results difficult. To combat this limitation, a visual comparison is performed between the produced oceanic features and the known general circulation conditions of the North Aegean sea.

Service Structure
Although research institutions and universities produce a significant amount of quality geospatial data during research activities or educational processes, they usually fail in the proper maintenance and distribution of the provided information [34]. Consequently, valuable scientific datasets are scattered among numerous departmental and laboratory computers, without any provision of exploitation by the wide research community. To overcome this issue and to keep up with the movement of Open Science and the pertinent National Plan of Greece [35], a Spatial Data Infrastructure (SDI) was developed within the Marine Remote Sensing Group (MRSG), University of the Aegean. The promotion of access to geographic information aims to encourage its publication as open data, i.e., freely available to everyone to use and republish as they wish, without restrictions due to copyright, patents, or other control mechanisms [36].
SDIs facilitate the discovery, access, exchange and sharing of geographic information and services among stakeholders from different levels in the spatial data community [37] and they are an essential tool for institutional data dissemination [38]. They also support decision-making processes and planning activities concerning environmental management and nature conservation [39]. An SDI is considered as an amalgam of data, technologies, standards, policies and people. Figure 3 depicts the general architecture of an SDI. In the following paragraphs, the components of the MRSG SDI and its operation to support data organization and diffusion of the mesoscale oceanographic feature identification task are presented.
Three datasets will be published to the service, the finalized skeletonized oceanic feature data, and both the raw and gradient magnitude datasets for reference. The oceanic features are in a vector file format in order to be able to store geographic information for each specific object. These data are stored in a PostgreSQL-PostGIS based database, with geographic object queries support, allowing for spatial and temporal interconnectivity between data from different dates.
All datasets related to the procedures of this paper are documented according to the ISO 19115 Metadata standard for geographic information [40]. In this way, the output datasets are uniformly described in a precise and complete manner and they are discoverable through open interoperability protocols (like OGC CS/W service, www.ogc.org/standards/cat). The GeoNetwork opensource platform (geonetwork-opensource.org) constitutes the metadata catalog of the SDI, as it provides a number of powerful capabilities, such as metadata editing, validating and search functions, monitoring and reporting tools, API and customization facilities and Free and Open Source Software principles support. and planning activities concerning environmental management and nature conservation [39]. An SDI is considered as an amalgam of data, technologies, standards, policies and people. Figure 3 depicts the general architecture of an SDI. In the following paragraphs, the components of the MRSG SDI and its operation to support data organization and diffusion of the mesoscale oceanographic feature identification task are presented. Figure 3. Generalized visualized architecture of an SDI. It integrates the Data layer, that consists of the infrastructure for storing data and metadata, Application layer, where the processing takes place and acts a medium between the user and the data, and End-user layer, the interface for interacting with the application layer and visualizing information.
Three datasets will be published to the service, the finalized skeletonized oceanic feature data, and both the raw and gradient magnitude datasets for reference. The oceanic features are in a vector file format in order to be able to store geographic information for each specific object. These data are stored in a PostgreSQL-PostGIS based database, with geographic object queries support, allowing for spatial and temporal interconnectivity between data from different dates.
All datasets related to the procedures of this paper are documented according to the ISO 19115 Metadata standard for geographic information [40]. In this way, the output datasets are uniformly described in a precise and complete manner and they are discoverable through open interoperability protocols (like OGC CS/W service, www.ogc.org/standards/cat). The GeoNetwork opensource platform (geonetwork-opensource.org) constitutes the metadata catalog of the SDI, as it provides a number of powerful capabilities, such as metadata editing, validating and search functions, monitoring and reporting tools, API and customization facilities and Free and Open Source Software principles support.
The creation of new metadata records is supported by an in-house developed tool that produces XML files to be imported into the GeoNetwork platform. The Purpose metadata attribute is used to discriminate the type of the dataset (i.e., raw SST data, gradient estimations, output skeletonized vector files). Other attributes-such as Title, Abstract, Keywords, Topic category, Lineage, Resource constraints, Data quality information, Metadata details and Contact information-are automatically Figure 3. Generalized visualized architecture of an SDI. It integrates the Data layer, that consists of the infrastructure for storing data and metadata, Application layer, where the processing takes place and acts a medium between the user and the data, and End-user layer, the interface for interacting with the application layer and visualizing information.
The creation of new metadata records is supported by an in-house developed tool that produces XML files to be imported into the GeoNetwork platform. The Purpose metadata attribute is used to discriminate the type of the dataset (i.e., raw SST data, gradient estimations, output skeletonized vector files). Other attributes-such as Title, Abstract, Keywords, Topic category, Lineage, Resource constraints, Data quality information, Metadata details and Contact information-are automatically filled in, according to the Purpose and Date attributes. The Distribution information metadata attributes are used to provide the appropriate links to the data repository and the map server (described subsequently).
For the visualization, cartographic navigation and downloading of the output products, OGC services are utilized (like WMS, WCS, WFS). Specifically, by exploiting the capabilities of the GeoServer open source map server, all datasets are automatically transformed to interactive map layers that end-users may access, display, query and download using a webGIS platform. The Metadata links and Data links attributes of the GeoServer layer editor are used to record the appropriate links to the metadata catalog and the data repository, respectively. End-users of the SDI may have full access to all datasets by activating different layers into a webGIS platform granted by the GET (Geospatial Enabling Technologies, www.getmap.eu) company.
The layers are based on the geographic web services already developed and published into the map server. The webGIS mainly consists of an interactive map and provides user-friendly cartographic operations for visualizing and manipulating thematic maps (like map navigation, overlay, base-map selection, etc.), as well as querying, measuring, printing and downloading operations. The adoption of metadata standards and open interoperability protocols for geospatial information has contributed to the increase in recognizability, accessibility and reusability of the products resulting from our work.

Interpolation Results
For the interpolation, a moving neighborhood of 200 pixels was used. This decision was taken in order to avoid to co-evaluate observations that are behind big land regions, such as islands and the gulfs of Chalkidiki peninsula, where the sea currents and general water conditions can differ vastly between one side of the land and the other. In addition to that, in order to avoid isolated pixels being influenced by further, more populated areas, a distance criterion has also been issued. Specifically, the neighboring pixels were restricted to be in a radius of 200 km. If less than 200 pixels were in that area, then only the subset was used for the estimation.
The accuracy of the Kriging combinations was assessed on 3 SST Level-2 datasets in the month of May, the 7th, 9th and 13th. This set of data was chosen because it contained a combination of both clear weather and heavily clouded imagery, and it was a good opportunity to test the performance of the methodology in variable conditions. For each date, all the combinations of the interpolation were tested, with the MAE and RMSE metrics (Tables 1 and 2, respectively). On the 7th, according to the RMSE index, the best interpolation method is by estimating the SST with the use of the previous day's data. Contradictory, the MAE index had the lowest error value when using only the specific days data. For the 9th, the RMSE index suggested that the best combination to use was the SST~(T-1) + (T-2), but the MAE showed that the SS~(T-1) is the best. Finally, on the 13th of May, both indexes agree that using all three dates for the SST estimation produced the most accurate results. Even though the results show that it should not matter what Kriging model should be used for the estimation, they fail to show the loss in detail. Using only the model SST~1, the estimated values would not encapsulate patterns that are included in previous more complete data, as with, for example, thermal fronts. Mesoscale and sub mesoscale phenomena have a temporal variability spanning days or even months [1]. As such, if a phenomenon is present in a specific date it should have an imprint, however small it is, in the next one. In an image where the interpolation formula considers only the individual dates data, and a large area is missing where a thermal front would be present then that area would be smoothed, with no information about the front. On the other hand, if on the previous day's dataset there was information about the front, then it would be beneficial for the modeling to include this information. Concluding, for this service the co-kringing formula of SST~(T-1) + (T-2) was decided to be used for data interpolation.
An example is given by Figure 4, for the 7th of May 2020 dataset. Here a comparison is drawn between raw data with masked pixels and the interpolation results. Originally large areas with missing information are located in Thermaikos Gulf and in waters southern of Chalkidiki. Figure 4B exhibits the interpolated values using the formula SST~(T-1) + (T-2), meaning the estimation of the SST values of the 7th in conjunction with the information of the 6th and the 5th of the month. An example is given by Figure 4, for the 7th of May 2020 dataset. Here a comparison is drawn between raw data with masked pixels and the interpolation results. Originally large areas with missing information are located in Thermaikos Gulf and in waters southern of Chalkidiki. Figure 4B exhibits the interpolated values using the formula SST ~ (T-1) + (T-2), meaning the estimation of the SST values of the 7th in conjunction with the information of the 6th and the 5th of the month.

Oceanic Front Detection Results
A detail presentation of the produced products is given in Figures 5 and 6, for the dates of 7th and 13th of May 2020, respectively. Starting from the gradient estimation images (Figures 5A and  6A), circulation formations are already visible. Major mesoscale phenomena of a spatial radius between of 20 to 50 km in the form of gyres in the Thracian Sea and the general direction of the Black Sea water is tracked, passing by Limnos and Imbros islands heading towards Sporades. In the Thermaikos gulf, submesoscale features are detected of a scale < 10 km in diameter, such as the river plume (East) and two gyres in the center of the gulf, created by the entering eastern currents. Other formations of the same scale are found in the Southern waters of the study area, specifically one

Oceanic Front Detection Results
A detail presentation of the produced products is given in Figures 5 and 6, for the dates of 7th and 13th of May 2020, respectively. Starting from the gradient estimation images (Figures 5A and 6A), circulation formations are already visible. Major mesoscale phenomena of a spatial radius between of 20 to 50 km in the form of gyres in the Thracian Sea and the general direction of the Black Sea water is tracked, passing by Limnos and Imbros islands heading towards Sporades. In the Thermaikos gulf, submesoscale features are detected of a scale < 10 km in diameter, such as the river plume (East) and two gyres in the center of the gulf, created by the entering eastern currents. Other formations of the same scale are found in the Southern waters of the study area, specifically one eastern of Skyros and another one western of Lesvos. Smaller eddies and meanders (scales of 3-10 km) are found in all the parts of the study area.
Generally, higher magnitude values are not associated to the spatial scale of the circulations (Figures 5B and 6B). The highest values are found in the areas dominated by either the river plumes, closer to the coastline and along the current created by the Black Sea Waters entering from the Dardanelles strait. Swirls, eddies and cyclonic flows are dominated by moderate magnitude signatures.
After applying the threshold, most of the noise is removed. Naturally, as it presented the bigger impact, information from the major formations is retained intact. Problems occur at lower scale phenomena, as the majority of them get fragmented. The gyres in Thermaikos, Skyros and Lesvos are visible, but some pixels failed to pass the threshold and they appear split in multiple parts. Other objects, that where visible but had low gradient magnitude values, either were rejected whole or only some pixels remained, and they appear as noise (Figures 5C and 6C).
The skeletonized data simplify the visualization of the resulting thermal fronts, but they introduce additional noise (Figures 5D and 6D). Especially in more complex objects, the thinning procedure tries to capture the general shape as well as the uneven extensions of them resulting, in many cases, in convoluted lines. Moreover, this technique further degrades the fragmented features and their identification without prior knowledge of the general circulation of the Northern Aegean is harder. For that reason, an overlay of the skeletonized data on either the original interpolated SST image or the gradient magnitude would be more useful for the average user to make photo-interpretation decisions. Especially the gradient magnitude products provide continuous information on steep SST changes indicating the position of both strong and weak circulation formations used in conjunction with the skeletonized data, users can track the fronts' position and where they are fragmented, the underlaid information can help them visualize these gaps.

Geospatial Web Service
The geospatial web service application can be accessed from the URL: http://mrsgsrv.lesvos.aegean.gr/sdi/?lang=EN. Search and information extraction options involve both easy to use and advanced options. The former lets users visualize metadata and information of specific data by simply selecting an object. The advanced option involves the issuing of queries, as well as

Geospatial Web Service
The geospatial web service application can be accessed from the URL: http://mrsg-srv.lesvos. aegean.gr/sdi/?lang=EN. Search and information extraction options involve both easy to use and advanced options. The former lets users visualize metadata and information of specific data by simply selecting an object. The advanced option involves the issuing of queries, as well as combining information from different information levels. For this search option, even though it runs on the database level as SQL queries, users do not need to be familiar with the SQL language as it is built in a simple web interface. The only arguments that should be stated are the information levels that should be used and the geographic relation between them.
In Figure 7, a preview is presented of the webGIS application interface. At the top toolbar, users can alternate between data viewing and file reading, as well as gain access to the user manual. The left panel is responsible for manipulating the visualized data information. All available layers can be selected or unselected, depending on whether they should be included in the map, from the left panel. In addition, the legend tab shows the scalebar of raster data and the feature labels, and the selected tab prints information on chosen areas, pixels and features. Directly above the interactive map, there are tools available for navigation, measuring and basic processing as well as printing the visible layers for rapid data retrieval. From the Search Results (down below), the query results are printed, and option is given for scrolling between different dates for visualizing the respective data automatically in the interactive map. Data access for downloading is given directly from the webGIS interface. Vector and raster datasets are available in different, universally used formats, like GeoTiff, ESRI shapefile, KML, etc. Furthermore, by utilizing the OGC protocols WCS and WFS, users can download data based on specific spatial and temporal criteria. The GetCoverage and GetFeature operations of the WCS and WFS protocols, respectively, enable the retrieval of raw data in a desired format along with their metadata. These calls can be used in order to automate daily data retrieval or for incorporating the datasets to scripts or services. An example of the service endpoint is given by http://mrsgsrv.lesvos.aegean.gr:8080/geoserver/ows?service=wms&version=1.3.0&request=GetCapabilities&na mespace=ocean_fronts.

Conclusions
In this paper, a new geospatial web service implementation is proposed, for providing spatial information on mesososcale and submesoscale processes in the North Aegean Sea. Daily oceanic feature identification maps are produced reliably with a gradient edge detection technique. These processes are encapsulated in a web-based SDI, and the final results are serviced to the end-users ready for use, either for visualization and quick decision purposes or for direct downloading.
The widespread availability of Sentinel-3 satellite data has enabled the extraction of more robust, fine and accurate information. Both temporal and spatial resolutions of the SLSTR sensor have proven that methodologies can be applied for smaller object identification and with smaller periodic cycles. These data are also ready to use and built for convenience, providing already processed data at Level-2, equipped with both cloud and land masks and pixel quality flags.
The proposed methodology combines reliability with simplicity. It is based in a process which is well studied from previous works, modernized, and optimized for current technologies and data. Capable of identifying features at scales bellow 5km in diameter, it has proven that can be applied to both small-and large-scale oceanographic applications. It is not without its flaws, however, as it is sensitive to high gradient magnitude signals it fails to capture features at more homogenous waters. The overlaid data are the vectorized detected oceanic fronts for the 7th of May 2020 and the corresponding SST gradient magnitude image. Users can add their own layers to the application, from their or other sources, provided that those datasets are already published as services (WMS, WFS, WCS, etc.). This feature gives the opportunity for immediate comparison of different datasets for validation or complimentary analysis.
Data access for downloading is given directly from the webGIS interface. Vector and raster datasets are available in different, universally used formats, like GeoTiff, ESRI shapefile, KML, etc. Furthermore, by utilizing the OGC protocols WCS and WFS, users can download data based on specific spatial and temporal criteria. The GetCoverage and GetFeature operations of the WCS and WFS protocols, respectively, enable the retrieval of raw data in a desired format along with their metadata. These calls can be used in order to automate daily data retrieval or for incorporating the datasets to scripts or services. An example of the service endpoint is given by http://mrsg-srv.lesvos.aegean.gr:8080/ geoserver/ows?service=wms&version=1.3.0&request=GetCapabilities&namespace=ocean_fronts.

Conclusions
In this paper, a new geospatial web service implementation is proposed, for providing spatial information on mesososcale and submesoscale processes in the North Aegean Sea. Daily oceanic feature identification maps are produced reliably with a gradient edge detection technique. These processes are encapsulated in a web-based SDI, and the final results are serviced to the end-users ready for use, either for visualization and quick decision purposes or for direct downloading.
The widespread availability of Sentinel-3 satellite data has enabled the extraction of more robust, fine and accurate information. Both temporal and spatial resolutions of the SLSTR sensor have proven that methodologies can be applied for smaller object identification and with smaller periodic cycles. These data are also ready to use and built for convenience, providing already processed data at Level-2, equipped with both cloud and land masks and pixel quality flags.
The proposed methodology combines reliability with simplicity. It is based in a process which is well studied from previous works, modernized, and optimized for current technologies and data. Capable of identifying features at scales bellow 5km in diameter, it has proven that can be applied to both small-and large-scale oceanographic applications. It is not without its flaws, however, as it is sensitive to high gradient magnitude signals it fails to capture features at more homogenous waters. Nonetheless, by not being a trained model but rather a series of image manipulation techniques, it offers great repeatability and correction of such defects should be simpler. Furthermore, as a non-data dependent methodology, it enables immediate application in other marine regions.
The proposed oceanic feature identification processing chain has produced promising results for operational services. Most of the important and impactful mesoscale and submesoscale circulation formations of the North Aegean have been successfully mapped. These data can be readily available through the webGIS interface for visualization and processing or they can be downloaded through open source means for use in various research fields. This is a powerful way to organize and distribute oceanographic data for daily use.
Moving forward, this service is planned to be updated and enriched with additional products as well as expanded for full coverage of the Greek Seas. The first major step towards the improvement of the (sub)mesoscale identification products is the inclusion of Sentinel-3's OLCI data, as chlorophyll fronts play a vital role in the fish population distribution. Data from the OLCI sensor can provide chlorophyll-a concentration estimations at a spatial resolution of 300 m, enabling the possible detection of even smaller scale processes. Furthermore, the SLSTR's swath completely overlaps OLCI's and both datasets can be used in conjunction, detecting both thermal and ocean color fronts, integrated in one product. These products, along with the fine-tuned existing ones, will be used to map the oceanographic processes in the North and South Aegean Sea daily, as well as in the Ionian Sea. Secondly, a statistical quality assessment of the results must also be performed with the use of in-situ measurements, so that any shortcomings can be pin pointed and corrected accordingly. Finally, as one of the goals of the geospatial web service is to work as a decision-making tool for maritime authorities and the fishing industry, a set of marine navigation tools will be incorporated, such as route planning and weather and wind forecast display.
Author Contributions: The conceptualization of study and the structuring of the methodology was performed by S.S. and K.T. The processing chain for extracting mesoscale features from remote sensing data was developed by S.S. D.K. designed and applied the kriging interpolation technique and M.V. and S.S. built the web geospatial service. All authors have read and agreed to the published version of the manuscript.