3.1. Radiometer SIC Products
Five radiometer SIC products selected for this study are summarized in
Table 1. Instead of including all available SIC products in our study, we selected products representing different SIC algorithm types. One of the products, the ARTIST (Arctic Radiation and Turbulence Interaction Study) Sea Ice (ASI) SIC algorithm using AMSR2 (ASI-AMSR2) 89 GHz data [
18], represents high-resolution SIC products retrieved with 90 GHz brightness temperature (
) channels in combination with low-resolution weather filters. The OSI-408-a is near-real-time SIC product based on the AMSR2 18.7 and 36.5 GHz
data used in many SIC algorithms [
19]. The OSI-450-a and OSI-458 are climate data records based on SSMIS and AMSR2 data, respectively [
8,
20]. The SICCI-HR-SIC (HR means high resolution) is a resolution-enhancing SSMIS product where SIC with 19.35 and 37.0 GHz
data (the same SIC as in the OSI-450-a product) is pansharpened using a SIC with finer-resolution 91.7 GHz data [
21]. These products have pixel sizes from 3.125 to 25 km. Some of the products are oversampled, like the ASI-AMSR2 3.125 km product. Products with larger pixel sizes are taken to be too coarse for the Baltic Sea, especially for the Bay of Bothnia and Gulf of Finland where sea ice forms every winter.
Descriptions of the SIC retrieval algorithms, e.g., tie point selection and applied weather filters, can be found in the references shown in
Table 1. All products, except the ASI-AMSR2, use the daily dynamic tie points method from Maass and Kaleschke [
6] to correct the land spillover already occurring at the
level, a dynamic open water filter, and the atmospheric correction of
s [
8]. The ASI-AMSR2 is based on fixed tie points and fixed open water filters [
18], and
s are not corrected for atmospheric effects. All SIC products have SIC values between 0 and 100%, i.e., they are truncated to the range between 0% and 100%.
3.2. MODIS OWSI Chart
For the validation of the radiometer SIC products over the Baltic Sea SIC, we derived an open water—sea ice (OWSI) chart using MODIS band 1 reflectance data which have a 250 m spatial resolution. The total number of the MODIS OWSI charts is 62, and they cover February to April 2019. This OWSI chart derivation follows our earlier study for the Barents and Kara Seas [
25]. In the following, the used MODIS datasets and their processing are first described, including the automatic and manual cloud masking of the MODIS reflectance data. Next, the OWSI chart derivation is described.
The MODIS spectrometer measures reflected solar and emitted thermal radiation in 36 spectral bands ranging from the visible to the thermal infrared (TIR) (0.415 to 14.235 μm) [
26]. The MODIS instrument is aboard NASA Terra and Aqua satellites. In this study, we derived an OWSI chart for the Baltic Sea using band 1 reflectance (
) data at a 250 m spatial resolution and the MODIS cloud mask product [
27]. MODIS sea ice surface temperature (IST) data and an RGB image with a 2-1-1 band combination (RGB211) were used in finding cloud-free MODIS swath datasets. The following Terra and Aqua MODIS products were used: the MOD/MYD/02QKM-MODIS/Terra Calibrated Radiances 5-Min L1B Swath 250 m, MOD/MYD/03-MODIS Geolocation Fields 5-Min L1A Swath 1 km, MOD/MYD/35_L2-MODIS Cloud Mask and Spectral Test Results 5-Min L2 Swath 250 m and 1 km, and MOD/MYD/29-MODIS Sea Ice Extent 5-Min L2 Swath 1 km. The MOD/MYD29 product contained IST data, and the cloud mask was derived from the MOD/MYD35_L2 unobstructed field-of-view flag [
28]. In the following, these products are only referred to with the ‘MOD’ prefix.
The MODIS products were rectified to the Mercator projection used for the Finnish Ice Service (FIS) chart using the MODIS Swath Reprojection Tool developed by the NASA LP DAAC (Land Processes Distributed Active Archive Center). In the rectification of the MODIS products, original pixel sizes (250 or 1000 m) were kept, and fixed Mercator grids were used. The size of the rectified data was 1330 km in the northing and 1120 km in the easting. The MODIS sensor scan angle () was calculated from the sensor zenith angle, Earth radius, and satellite orbit altitude. Raw data from the bands 1 and 2 were converted to top-of-the-atmosphere (TOA) sun angle corrected reflectances ( and ). The land mask for the Baltic Sea study area was retrieved from the MOD44W (Land Water Mask Derived L3 Global 250 m) product.
In all MODIS data analyses herein, the retrieval of the MODIS OWSI chart
was limited to a maximum of 40° (the overall maximum is 55°). At a
of 40° the MODIS spatial resolution in the across-track direction decreased by a factor of around two, e.g., from 250 m at the nadir to ~500 m at 40° [
29]. In the MOD29 and MOD35 products, daylight data were defined as data collected with a
less than 85° [
27,
28]. We used here a somewhat smaller limit of 80° as at high sun zenith angles the
term in the TOA reflectance calculation is highly variable, and, thus, the accuracy of the estimated
is critical. This
limitation excluded all the MOD02 data from December 2018 to January 2019 in the Bay of Bothnia. Therefore, we used here only MODIS data acquired over the Baltic Sea from 1 February 2019 onwards to end of April 2019.
The following procedure was used to select mostly cloud-free MOD02 data over the sea ice-covered part of the Baltic Sea. First, all available MODIS IST data were processed over the Baltic Sea from 1 February to 30 April 2019. Visual analysis of the IST charts with the MOD35 cloud mask was conducted to select times for the MOD02 data, which possibly had large cloud-free areas over the sea ice-covered part of the Baltic Sea. Next, RGB211 images with and without s cloud mask were used in the further screening of the MOD02 datasets for the OWSI chart derivation. We ended up with 20 MOD02 swath datasets for February, 21 for March, and 21 for April. Here, a MOD02 dataset was either one MOD02 swath or two swaths stitched together. For some days, there were two datasets (Terra and Aqua) and their temporal difference was mostly 10 or 15 min and at maximum 100 min.
Cloud masking was started by designating pixels that had ‘confident clear’, ‘probably clear’, or ‘probably cloudy’ values in the unobstructed field-of-view quality flag of the MOD35 product as cloud-free. This logic followed the MOD29 IST product generation with the target to increase the number of retrievals balanced against the cloud conservative nature of the cloud mask [
28]. A visual analysis of the resulting cloud mask over the RGB211 image showed that sometimes cloud-free pixels of open water or the thin ice within led to larger open water/thin ice areas being classified as cloudy, and cloud-free sea ice, in general, was classified as cloudy. This is because NSIDC’s Near-real-time Ice and Snow Extent (NISE) product with a 25 km resolution was used to initialize the sea ice background flag for directing the cloud masking algorithm processing flow [
28]. Open water or thin ice when the NISE sea ice flag is set may be classified as cloudy, and sea ice when the flag is not set may also be classified as cloudy. Methods have been developed for the correction of these cloud mask-induced errors in the OWSI classification, e.g., in [
30], but we decided to correct these errors by the manual editing of the cloud mask. For the sea ice-covered Baltic Sea, which has a rather small area, the manual cloud mask editing for one MODIS dataset was not an overly laborious task, and it also served as a quality control step in the OWSI chart processing.
Before the manual cloud masking, all single or groups of two cloudy pixels within the NISE sea ice flag set area were first removed. These were assumed to be erroneous cloud detections over the open water or thin ice. Next, the cloud mask was aggregated into 5 km blocks (5 by 5 1 km pixels), and if there were more than or equal to three cloudy pixels (12% of all 25 pixels), then the whole block was flagged as cloudy. Using 5 by 5 km blocks for the cloud mask, we could better identify large cloud-free areas and discard areas where the cloud mask was tattered. Small cloud-free areas (‘holes’ in the mask) less than ten blocks in size were flagged as cloudy. Finally, the following manual editing procedures could be conducted: filling holes, removing erroneous cloud mask elements, and masking arbitrary polygonal areas as cloudy or clear. The manual editing was conducted using the RGB211 image. An example of the original and manually edited cloud masks is shown in
Figure 1.
We aimed to develop a simple OWSI classification method based only on the threshold
(bandwidth 620–670 nm, red). We only used the
data as it had the best possible MODIS resolution, 250 m. The
data (841–876 nm, near-infrared (NIR)) also had a 250 m resolution, but we evaluated that it was not needed in the OWSI classification. In addition, we have previously determined an
-based OWSI classification for the Barents and Kara Seas [
25]. In the Baltic Sea, MODIS
data have been used to map the sea ice extent in the Gulf of Riga [
31]. In this study, the threshold for the OWSI classification was determined automatically for each swath dataset based on a sampled
histogram which was bimodal, showing OW and SI peaks. More accurate OWSI classification is also possible using other MODIS reflectance bands and shortwave and thermal infrared bands, e.g., as in [
28,
30,
32], but this is at the cost of decreasing the resolution to 500 m or 1 km. In the MOD29 product, for daytime, there is sea ice extent data based on various reflectance tests, including the threshold
, at 1 km resolution [
28]. An algorithm called IceMap500 has been developed to generate 500 m sea ice extent maps with a systematic accuracy over 90% [
30]. Machine learning methods have also been developed to classify OW, SI, and clouds; for example, Jiang et al. [
33] constructed a training sample library and used a Multi-Feature Level Fusion Random Forest classification algorithm that integrated multiple features from the reflectance bands 1, 2, and 7. Their SI chart had a 250 m pixel size, but this required the interpolation of the band 7 data at a 500 m resolution to 250 m.
For the determination of the threshold for the OWSI classification in the Baltic Sea, we manually selected samples in the cloud-masked RGB211 images from February 2019 for OW, thin ice (THI), and snow-covered sea ice (SSI). These samples represented winter conditions in the Baltic Sea. The size of the sample windows for the OW and SSI was 5.25 by 5.25 km (21 by 21 pixels), and it was 2.75 by 2.75 km (11 by 11 pixels) for the THI. The number of sample windows was 210 for the OW, 278 for the THI, and 128 for the SSI. Another set of samples were selected in April 2019, which represented springtime melting sea ice conditions. Here, we selected windows for the OW (219 windows), SSI (35 windows), and melting sea ice (likely mostly snow-free) (MSI, 385 windows, 11 by 11 pixels window size).
The resulting
probability density function (pdf) for the OW is very narrow (see
Figure 2); for example, the 95th percentile of the
is only 0.070 in the winter data. The SSI is well separated from the OW, and the 5th percentile of the
is around 0.55. In the winter data, the THI has a large variation, from 0.16 to 0.51 by the 5th and 95th percentiles. The MSI likewise has a large variation, from 0.073 to 0.42. Setting the
threshold to 0.10, only 0.12% of the OW samples in the winter data are classified as THI, and only 1.0% of the THI is classified as OW. In the spring data, the OW and MSI misclassifications are 0.46% and 4.5%, respectively. Unfortunately, it is not possible to separate fully the OW and very thin ice which both have small
values. In addition, the
data (TOA) here do not have atmospheric corrections, and in the springtime, water turbity close to the coast and river estuaries may raise
. Therefore, it is better to use somewhat too high than too small of a threshold. We decided to use the threshold of 0.10 for all
data as it gives very small OW misclassifications in our manually selected sample data:
This threshold was manually determined. In the resulting OWSI charts, very thin ice (thickness: a few cm) and very wet sea ice may be classified as OW. Fog and thin unmasked clouds over OW can also be misclassified as sea ice. In addition, the mixed-pixel effect affects the OWSI classification as even a small fraction of sea ice, especially that of snow-covered ice, within the 250 m pixel leads to assigning the SI class to the pixel. This would further lead to the overestimation of SIC when the MODIS data are used to calculate SIC at radiometer data resolutions. An example of the resulting OWSI charts in shown in
Figure 3 for the same date and time as the RGB211 images in
Figure 1.
3.3. Sentinel-2 OWSI Chart
Sentinel-2A and Sentinel-2B (S-2A/B) carry the MultiSpectral Instrument (MSI) with 13 different spectral bands from VNIR (visible–near-infrared) to SWIR (shortwave infrared) with band central wavelengths varying from 443 to 2202 (S-2A) / 2186 (S-2B) nm [
34,
35,
36]. The spatial resolution for different bands is either 10, 20, or 60 m. The swath width of MSI is 290 km. The Sentinel-2 mission supports Copernicus land, marine, security, emergency, and climate change services and studies, including the monitoring of vegetation, soil, and water cover, as well as the observation of inland waterways and coastal areas. The MSI data can also be used for the remote sensing of sea ice, e.g., for open water–sea ice (OWSI) classification, but MSI lacks longwave infrared bands for good cloud masking over sea ice, e.g., 3.9 and 11 µm TIR bands, as in MODIS.
There are two end-user MSI products: Level-1C top-of-atmosphere (TOA) reflectance and Level-2A bottom-of-atmosphere (BOA) reflectance. We used here the L1C product with the processing baseline 02.07 version. The L1C product had 100 by 100 km fixed tiles (or granules) in the UTM/WGS84 projection. The L1C product had opaque (dense clouds) and cirrus cloud masks [
37]. The cloud masks were processed with data sampled at a 60 m spatial resolution for all spectral bands. The land masking of the L1C tiles was conducted here with the ESA SNAP software (version 9.0) (fractional land/water mask). The fixed land mask was based on the ESA Global Land Cover Map data. It did not identify all the small islands in the Bay of Bothnia.
We aimed to use the 10 m resolution bands (best possible resolution)—B2 (490 nm; blue), B3 (560 nm; green), B4 (665 nm; red) and B8 (842 nm; NIR)—from the L1C product for the OWSI classification over the Bay of Bothnia. B4 was equal to MODIS band 1 (bandwidth of 620–670 nm) which was used here for the MODIS OWSI classification.
As S-2 cloud detection over sea ice may be inaccurate, we searched S-2 data over the Bay of Bothnia only for the dates when the MODIS data were mostly cloud-free and the sun zenith angle was mostly below 80°. The Bay of Bothnia was covered with five 100 km tiles in the UTM34 projection. Our data search for February-April 2019 resulted in S-2 imagery for 13 dates; see an example in
Figure 4.
S-2 data have been used previously in sea ice remote sensing studies either as a validation dataset or to classify sea ice for further analyses, like for lead statistics, e.g., in [
38,
39,
40]. Ludwig et al. [
38] evaluated the accuracy of their MODIS and AMSR-E data-based SIC product with a 1 km pixel using 79 cloud-free (visually determined) S-2 L1C products acquired mainly over the Arctic FYI. They investigated the
s of the visible spectrum bands 2–4 for open water (OW), thin ice, and thick ice classification. These bands did not show significant differences in the classification, and there was a clear distinction between the thin and thick ice and a clear distinction between the OW and sea ice (SI) if only thick ice was present. They decided to use only the band 4 reflectance (
) for the surface type classification and selected specific classification thresholds (OW vs. thin ice and thin vs. thick ice) for each
image instead of using global thresholds for all images. Muchow et al. [
39] investigated lead-width distribution for the Antarctic sea ice using 20 cloud-free L1C products. Lead detection was conducted by thresholding
. For the threshold determination, they selected manually sample data for OW and the following four ice types: nilas (NI; thickness < 10 cm), gray sea ice (GI, 10–15 cm), gray–white ice (GWI; 15–30 cm) and sea ice covered with snow. The sample data were used to calculate
histograms for each surface type, and then a summation of Gaussian functions was fitted to them. The threshold for each surface type was then determined from the intersection of two curves adjacent to each other. For the lead identification, two different thresholds were used: 0.10 for leads covered with OW and 0.17 for leads covered with OW and NI. The threshold of 0.10 classified 29% of NI as OW and the threshold of 0.17 classified 11% of GWI as NI. Wang et al. [
40] developed a decision tree classification of OW and various ice types in the Liaodong Bay of the Bohai Sea using S-2A/2B data. Using one S-2B image, they selected manually sampled areas for different surface types: OW, NI, GWI, and landfast ice (LFI). Utilizing mean reflectance spectra (
to
), a decision tree classification was determined for these surface types. This classification tree included 10 m and 20 m bands. The total accuracy of the classification was 88% according to visual validation. These previous sea ice classification studies demonstrate that OWSI classification in the Baltic Sea is possible with the S-2 data.
We started the development of the OWSI classification method by first investigating the pdfs of the 10 m reflectances of
(490 nm; blue),
(560 nm; green),
(665 nm; red), and
(842 nm; NIR)), as well as the normalized difference water index (NDWI) in the Bay of Bothnia. The NDWI is calculated as follows [
41]:
The NDWI detects open water features against soil and terrestrial vegetation features (they have a smaller NDWI than OW). It is noted that the normalized difference snow index (NDSI) could not be calculated from the S-2 10 m reflectance data. The pdfs were calculated using imagery of one tile in the northern Bay of Bothnia (tile 34WFT) acquired on 5 March and on 6 April 2019. The first date represented winter conditions (the air temperature
at the Kemi 1 lighthouse was −10 °C), and there were large areas of thin ice (<15 cm) according to the Finnish Ice Service (FIS) ice chart; the second date had spring, early melting conditions (the
was around 0 °C) when there was no thin ice in the Bay of Bothnia. For 5 March, all four 10 m
pdfs show a distinct peak at small
values, representing OW, and another peak at high
values, representing snow-covered SI, mostly LFI; see
Figure 5. There are a few smaller peaks, although they are not very distinct, representing different types of level, deformed, and snow-free thin ice. There are no clear
thresholds between the OW and thin ice. For
and
, the OW peak is somewhat narrower than for
and
. The
pdfs on 6 April also have a OW peak, but they have much a smaller snow-covered SI peak as large areas of LFI snow had melted; see
Figure 5. The snow melting resulted in large fractions of
values in the 0.2–0.3 to 0.6 range. There is now a clear
local minimum after the OW peak, as there was no new thin ice on 6 April. The NDWI pdf shows a clear LFI peak and smaller OW peak but no better discrimination of surface types than the single
pdfs.
Likely, there was not much difference that affected which
band we would use for the OWSI classification, but
and
could have been better than
and
due to their narrower OW peaks. Therefore, we selected
for further analyses as it has also been used previously for lead detection and SIC estimation [
38,
39].
First, we investigated the statistics for the OW on 5 March and 6 April by visually selecting 195 OW rectangles from four S-2 tiles acquired on 5 March and 194 rectangles from five tiles on 6 April. The size of the rectangles was from 0.15 to 870 km2. On 5 March, the tile-based mean varied only from 0.047 to 0.056, the mode was either 0.05 or 0.06 when the pdf bin width was 0.01, and the 99th percentile was from 0.057 to 0.065. On 6 April, the tile based statistics were slightly larger: the mean varied from 0.048 to 0.059, and the 99th percentile varied from 0.056 to 0.078. The mode was again either 0.05 or 0.06.
These OW
statistics suggest that fixed thresholds for the OWSI classification are not feasible. A tile-based threshold is a better approach. This tile-based OW threshold takes into account variation in the sun zenith and sensor zenith angles (TOA reflectance does not correct for atmospheric effects, which depend on the sensor zenith angle, i.e., the path length). In the above rectangle selection, we observed that it was sometimes difficult to visually discriminate between OW and very thin ice, i.e., dark nilas and new ice (including grease ice and slush), which had same kind of ‘dark’ appearance (in level and texture) in the
image. There was also sometimes thin fog over the OW and thin ice, which made their identification more difficult and raised the
. Most notably, there were areas where both the OW and very thin ice had a similar
level (and no texture). The OW and new ice also had very similar S-2 spectra in [
40]. These observations show that full discrimination between OW and very thin ice is not possible, and a
threshold for the OWSI classification will classify some unknown part of thin ice to the OW class. This is likely also the case for OW and melting wet ice in late spring.
Next, we needed to find a way to determine the threshold for each L1C tile. One way to conduct this is to manually sample OW data from a tile and to calculate the 99th (or a lower 95th) percentile of . Unfortunately, this manual sampling is time consuming. We tested fitting a mixture of Gaussian pdfs to the OW and SI data in a tile and used the intersection between the pdf with the smallest mean (representing OW) and the pdf with the second smallest mean as the threshold. The fitting seemed to work fine when there was a local minimum after the OW peak but not when this was not the case. Therefore, we decided to set the threshold manually with the help of the pdf. The manually set threshold varied from 0.06 to 0.11, and it was larger in winter, from 0.09 to 0.11, than in the melting season (12 April onwards), from 0.06 to 0.08. The manual setting was easier in the melting season as new thin ice was not present.
Before the OWSI classification, the opaque (dense clouds) and cirrus cloud masks were combined, and the combined mask was manually edited. This editing was conducted with the help of an RGB image with
-
-
bands (approximating a true color image); see an example in
Figure 4. In the manual cloud masking editing we rather masked too much than too little to cover for certain all cloudy and cloud shadow areas. A visual comparison of the cloud masks of the
and RGB images showed that the opaque cloud mask was sometimes showing erroneously cloudy areas over LFI. In some
images from later in April, there were increased
values over open water at river mouths due to water turbidity. These areas were included in the manual cloud mask as they would otherwise be classified as SI.
An example of the resulting OWSI charts is shown in
Figure 6. Here, the S-2 tile data are the same as in
Figure 4. All tile-based OWSI charts for the same S-2 overpass (max. five tiles) were combined to create an OWSI mosaic over the Bay of Bothnia.
3.5. Co-Location and Comparison
For comparison with the MODIS OWSI chart, the SIC products were rectified to the Mercator projection of the MODIS chart with nearest neighbor sampling. Inside each SIC product pixel, the MODIS SIC was calculated as the fraction of sea ice pixels to the total number of pixels. To reduce the land spillover effect, only those radiometer SIC pixels which had, at maximum, 10% land in the MODIS land mask were included in the comparison.
In similar way, the SIC products were rectified to the UTM34 projection of the S-2 OWSI chart, and then the S-2 SIC was calculated within the SIC product grid cells which had, at maximum, 10% land in the S-2 land mask.
In selecting radiometer data pixels for comparison, the status flag of the SIC retrieval in the OSI-450, OSI-458, and SICCI-HR-SIC products was utilized. This flag showed when an open water filter and a land spill over filter had been applied. We only selected pixels with no filters applied in order to only have ‘true’ retrieved open water pixels in the comparison. There was also a status flag in the OSI-408-a product, but it did not have the needed information. The ASI-AMSR2 product did not have a status flag at all. For these two products, we excluded pixels with a radiometer SIC of exactly zero. These pixels were assumed to be mostly the result of applied filters. Keeping the rejected pixels would have had a large effect on the comparison statistics.
The comparison between the gridded MODIS or S-2 SIC and co-located radiometer SIC was conducted by calculating the mean difference of the radiometer SIC minus the MODIS/S-2 SIC, standard deviation of the difference, and linear regression line and coefficient of determination of the regression ().