OGC Testbed-16: Machine Learning Training Data ER

The first challenge was data analysis. Specifically, how to compare two different variables, represented by their respective maps, to estimate their correlation. This is a basic step in feature selection. The process is to select those features that contribute most to the prediction variable or output in a Machine Learning algorithm. This process helps finding relevant features, leaving irrelevant features out of the training data set.

One of the most common methodologies used in data analysis is Mean Squared Error (MSE). This error estimate measures the average of the squares of the errors between two variables. However, as this error estimate only considers local values and ignores shapes of isolines, which is one of the main characteristics in maps, this method is unsuited for map comparison. For example, two maps showing very similar isolines but different absolute values will be considered unrelated by MSE methodology.

In order to overcome this challenge, the participants researched new methodologies in the field of image processing. They determined that Structural Similarity Index (SSIM) was considered suitable for the task. The SSIM index is a method for measuring the similarity between two images. The SSIM index can be viewed as a quality measure of one of the images being compared provided that the other image is regarded as of perfect quality. This is a method for predicting the perceived quality of digital television and cinematic pictures, as well as other kinds of digital images and videos. However, SSIM can be applied to find common patterns and contours between two images and estimate their degree of similarity.

As displayed in the following example, when applying SSIM to comparing two different images the value will be indicative of the similarity of the images. When applying SSIM to maps, the hypothesis is that SSIM will be a better indication than MSE-related metrics for variable correlation or for the accuracy of ML predictions.

The SSIM method was tested with the Web Map Service (WMS) accessible Petawawa data. The Petawawa forest dataset provides a series of maps displaying the values of several measurements, such as volume, height or biomass. The dataset is for a relatively small and compact area but offers up to 26 different measurements with a spatial resolution of 30 meters and the quality of the data is assumed to be very high: that is, close to the real value distribution. Further, the maps come directly as gray scale that can be analyzed directly.

Both methods (MSE and SSIM) were run for all maps. SSIM is an index, topped at 1.0, which indicates identical maps. Values over 0.9 indicate a high similarity (90 percentile), and lower than 0.5 indicate very low similarity (10 percentile). In contrast, MSE is an absolute value, with 0 indicating identical images, values under 500 indicate very few differences (10 percentile), and over 5000 indicate substantial differences (90 percentile).

The interest in using SSIM becomes obvious when analyzing maps that show related variables, such as Lorey’s height, the weighted mean height whereby individual trees are weighted in proportion to their basal area (RF_PRF_LOREYSHT), in relation to Co-dominant - dominant height (RF_PRF_CD_HT) and Top height, the average of the largest 100 trees per ha (RF_PRF_TOPHT). These variables all measure forest heights and should show a strong correlation.

When comparing Lorey’s height in relation to Co-dominant - dominant height, both methods indicate high similarity, with SSIM value at 0.93 (very close to the maximum 1.0) and MSE at 161 (well below the 10th percentile value of 500).

However, when comparing Lorey’s height in relation to Top height, only SSIM, with a value of 0.92, indicates strong correlation. MSE gives a value of 1036, indicating that there might be some correlation. The reason for this misalignment is due to the fact that in spite of both maps showing similar shapes, the absolute values of the variables differ significantly.

Similar issues arise when comparing unrelated variables. For example, when comparing Top height (RF_PRF_TOPHT) in relation to Basal area in the Pole tree size class (10-24cm) (RF_PRF_BAPOLES), both methods indicate very low correlation or none at all, with values of 6746 for MSE and 0.39 for SSIM.

However, when comparing Basal area in relation to Quadratic mean diameter for all trees >9.1cm (RF_PRF_DBHQMERCH_MASKED) only SSIM results, with a value of 0.37, clearly indicate no correlation. The MSE value of 1666 which, although not pointing towards a clear correlation, does not reject the relationship either.

All the previously described tests were carried out using the WMS accessible Petawawa dataset. This is high quality data with one continuous variable per layer and with maps served in grayscale. However, this dataset’s geographic extent is very limited and conclusions therefore do not necessarily extrapolate to countrywide maps.

The CWFIS offers an alternative source of data to analyze the available data and showcase the Structural Similarity Index. The CWFIS maps are countrywide and many refer to discrete variables with their own color code, as ranges of the real-world variable. Moreover, many of these maps are related, such as the Fire Weather Index System or the Fire Behavior Prediction System, and offer a good testbed to analyze correlations.

The first challenge analyzing CWFIS maps is to convert them to a scale that is numerically significant and comparable by an algorithm. Maps in CWFIS have color codes that need to be transformed to a grayscale in which the higher volume according to the color code corresponds to either white or black and vice-versa. This challenge can be overcome as the CFWIS WMS offers a GetLegendGraphic option to associate a color to a value range, which then returns a "png" image with the map legend. This option is useful for visual inspection of the data. Alternatively, a WMS instance could offer a Styled Layer Descriptor (SLD) document. This is an XML file encoding of the legend information in a machine-readable format and accessible using a GetStyles operation. The SLD document can describe rules for rendering, and these rules can contain MinScaleDenominator and MaxScaleDenominator elements specifying the numerical values of the scale ranges. This is a valid method to communicate scale values in a machine-readable format using OGC standards. Both GetLegendGraphic and GetStyles operations are specified in an extension to the WMS standard. This extension is called OpenGIS SLD Specification, which is an optional extension.

Unfortunately, this extension does not seem to be fully implemented in NRCan’s services. The GetLegendGraphic operation works properly and provides a descriptive "png" image that a user can visually interpret. However, when the GetStyles operation is used to retrieve an SLD document, the retrieved document is empty. This means that an ML algorithm does not have a numerical reference to interpret the colors and shades in a map and cannot make proper inferences and regressions. In order to overcome this limitation, the color value is considered the value of the variable.

In the following examples, the lowest value is set to black and the highest to white, while the remaining gray gradation is assigned evenly through the different ranges. The following maps show the gray-scaled maps for Wind Speed, Fine Fuel Moisture Code and Initial Spread Index.

According to the Fire Weather Index structure graph, the combination of Wind Speed (WS) and Fine Fuel Moisture Code (FFMC) maps generates the Initial Spread Index (ISI) map. This derived map when combined with the Buildup Index (BUI) generates the Fire Weather Index (FWI). Applying both MSE and SSIM methods indicates a strong correlation between FFMC and ISI (0.96 according to SSIM and 114 according to MSE) and between ISI and FWI (0.97 according to SSIM and 70 according to MSE). It is even possible to apply a simple regression and reverse-engineer the generation of ISI, as depicted in the following map, with a result pretty close to the original ISI map. Comparing this result to the original ISI map it scores relatively well (0.96 and 80 respectively).

MSE and SSIM methods applied to all CWFIS maps draw consistent correlations, grouping maps by families according to their structural similarity. For example, Head Fire Intensity (HFI) is closely related to Rate of Spread (RoS) and Total Fuel Consumption (TFC) but they are also loosely related to FFMC, ISI or WS. This makes sense as HFI, RoS and TFC belong to the Canadian Forest Fire Behavior Prediction (FBP) System and FFMC, ISI and WS belong to the Canadian Forest Fire Weather Index (FWI) System).

The prototype OGC API - Records implementation was considered a useful API to browse and discover datasets, as well as provide all the required information to implement data retrieval by a ML algorithm.

Additionally, for a dataset to be useful as a source to train a ML algorithm the dataset needs the following attributes:

NRCan’s existing OGC Web Services can be useful for ML applications as long as they comply with the above requirements.

Training data, also known as ‘initial data’ or ‘sample data’ are needed in ML to build mathematical models that can then be automatically detect certain features in the data. High quality training data is essential to produce reliable mathematical models. However, establishing these initial data sets can be incredibly daunting. One of the tasks in OGC Testbed 16 was to determine how well OGC standards can support ML in the process of building algorithms to automatically detect waterbodies in Sentinel-1 scenes.

For training datasets two earlier established waterbody databases were used:

A) WFS3 - Quebec_Rivers_and_Lake_TB15_ML_D104.

B) WMS - Lakes, Rivers and Glaciers in Canada - CanVec Series - Hydrographic Features.

The Quebec Rivers and Lakes database (database A) created as part of the OGC Testbed15 ML activity was not regarded as it did not contain a sufficient number of labeled waterbodies. The labeled hydrological features, one of the 8 themes of the feature classes, included in the CanVec Series (dataset B) are watercourses, water linear flow segments, hydrographic obstacles (falls, rapids, etc.), waterbodies (lakes, watercourses, etc.), permanent snow and ice features, water wells and springs. This NRCan vector product is considered the best available geospatial data source covering Canada. The Hydrographic features theme is considered quality vector geospatial data (current, accurate, and consistent) of Canadian hydrographic phenomena. CanVec aims to offer a geometric description and a set of basic attributes on hydrographic features that comply with international geomatics standards, seamlessly across Canada. According to the metadata provided via the CanVec WMS endpoint, the hydrological features are from:

“the best available data sources covering Canadian territory, offers quality topographical information in vector format, and complies with international geomatics standards. CanVec is a multi-source product coming mainly from the National Topographic Data Base (NTDB), the Mapping the North process conducted by the Canada Center for Mapping and Earth Observation (CCMEO), the Atlas of Canada data, the GeoBase initiative, and the data update using satellite imagery coverage (e.g. Landsat 7, Spot, Radarsat, etc.)”

Although the CanVec hydrological features is a very comprehensive, state-of-the-art hydrological product, using these data as label data for training datasets might be less ideal for several reasons: 1) Merged product from different sources, 2) No date range of which imageries are used to produce this labeled dataset are provided, and 3) Metadata on uncertainty or data accuracy is missing. To be more specific:

Water bodies might be less dynamic over time but Figure 18 illustrates how delineation of water bodies can differ over time or by using different datasets while extracting the same water body. Making additional metadata available through the WMS endpoint that includes what raw data sources were used over which period would most likely always improve the accuracy of the ML model.

Technical details

Both labeled and training datasets were provided as non-overlapping 512 x 512 map tiles. The ML algorithm applied uses raster data only, so data tiles for the vector data (*.mvt or *.geojson) were not needed and therefore not provided through an OGC service. Also, for the waterbody vector data, the styling was removed so the data used does not contain contours anymore. Furthermore, no data values were set to display transparent. The dataset was provided in the same projection and tiling scheme and is currently available through the following URL: https://eratosthenes.pvretano.com/cubewerx/cubeserv/default/ogcapi/mysql_tb16

Derived datasets are authoritative but might be less suitable to use as labels in the ML chain when no adequate metadata is included. This often occurs when labeled data is initially created for another purpose than to train ML models. In general, for labeled data make sure to enrich label datasets with metadata when using OGC WMS standards. Metadata should include: a) Source of raw data, time period of raw data used to generate the label data, information on data accuracy, and b) precisely which raw data files were used (including a link to the data source).

Initially RADARSAT-1 (1995 – 2013) Synthetic Aperture Radar (SAR) was selected to use as a training dataset. However, a geo-reference issue resulting in that a systematic offset was detected after applying the default pre-processing steps used in Sentinel Application Platform (SNAP). These steps were:

Following this procedure for RADARSAT-1 imagery, the resulting raster had an offset. This offset was caused by the poor orbit state vectors for RADARSAT-1. The orbit state vectors are updated every 8 minutes but the image acquisition takes less than 20 seconds. To estimate the satellite positions every 20 seconds by interpolation of the vectors is not accurate enough as this leads to incorrect pixel locations in terrain correction. This is a known issue with the RADARSAT-1 product. With RADARSAT-2 the orbit update period is 2 seconds.

Having identified the problem with pre-processing RADARSAT-1, the OGC Testbed16 participants decided to use Sentinel-1 (S1) SAR data for the remaining part of this project. Sentinel-1 data is considerably more difficult to display than other satellite data in near real time (which are no trivial task either). Each pixel in an image needs to be converted to meaningful physical quantities (the backscatter coefficient) and geopositioned individually. There is no global mapping which can take care of this. Orthorectification is even more complicated as it requires the use of a digital elevation model (DEM), and again, adjusted for each pixel individually. Many service providers do give the option of downloading orthorectified S1 imagery as end users might want to use a different DEM for their specific needs. As such S1 data is only limited available as GeoTIFF. Note that the Sentinel Hub streamed imagery from Amazon Web Services (AWS) is not (yet) orthorectified either but AWS will support this soon.

Although Sentinel-1 data is not new to Sentinel Hub, the fast AWS architecture and cloud optimized GeoTIFFs allow for faster response times compared to the current service running on a different platform with original, non-cloud-optimized files.

In general, serving Sentinel-1 data via a WCS endpoint can be challenging as the data is marked as 16-bit grayscale. As such the data will generate 8-bit JPEG grayscale overviews which will be almost entirely black for Data retrieval, and a contrast-stretching RasterSymbolizer will need to be applied to make the data visible for Map access or else it will look entirely black.

Level 0 products of the side viewing sensor of Sentinel-1 SAR records variability of the surface height of the earth. So a correction through an orthorectfication process needs to take place by using a digital elevation model before being able to host this data as a WMS or WCS endpoint. Therefore Sentinel-1 data was pre-processed by NRCan. The SNAP pre-process steps applied to Sentinel-1 SLC data are:

On the output channels, VH and VV were assigned to channel 1 and channel 2 respectively. Then, for the RGB GeoTIFF displays the following was assigned to the RGB channels:

R: VH (Channel 1)

G: VV (Channel 2)

B: VV - VH (Channel 2 - Channel 1)

In total 89 sceneries from Sentinel-1 SLC data were processed, covering a large part of Canadian territories.

Once processed the properties for the GeoTIFF were:

$ tiffinfo S1B_IW_SLC__1SDV_20170620T143323_20170620T143353_006135_00AC6F.tif

The file information shows that for a non-standard tag 65000 a big XML document in some random format was included. This does not pose any problems except that the document cannot be parsed and the data-capture time is available inside that document. However capture time is also in the filename, although that is not parsed either.

Following the SNAP procedure, a technical flaw in creating a GeoTIFF is that the TIFF specifies a "Samples/Pixel: 2" tag but it doesn’t include an "Extra Samples: 1<unspecified>" tag. The TIFF spec requires the "Extra Samples" tag in this place, but it isn’t a problem since our code and the TIFF library assume the information of this tag.

A practical flaw that does have a big impact is the "Rows/Strip: 8167" element. To understand this, the developer has to understand the bitmap structure of TIFFs. There are three bitmap schemes TIFF supports, single strip images, stripped images and tiled images. For the sake of simplicity tiled images are not discussed here. For single stripped images, all of the bitmap is stored in one large block. A general recommendation is that no one strip should take more than 8 kilobytes uncompressed, which, with black-and-white images limits you to 65,536 pixels in a single strip. For stripped images, horizontal blocks of the image are stored together, and more than one strip is joined vertically to make the entire bitmap. However, if you allocate more strips than needed, then automatically the file size increases significantly, and that occurs in the pre-processing steps. Given the pixel size of the image there should be something like 3 strips. As such, the TIFF library reads scanline-organized images in Strips, which should be around 256 KB in size. However, for these preprocessed images, the Strip size is 8167 scanlines, which equals ~642 MB. This is a significant waste of memory, but worse than that, the TIFF-lib Scanline interface is used to read the image, which re-reads the entire 642 MB of data for each scanline.

Reading this image organized as "Rows/Strip: 8167" takes 7531 seconds (2.1 hours). Reading it organized as "Samples/Pixel: 2" takes 2.3 seconds.

Either files need to be reprocessed instead of using the original source data or the way the program reads Scanline TIFFs has to be changed. Once done, performance should be much faster but still wastes 642 MB of memory. This is quite a lot of machine resource where the RAM is divided between many processors.

When reprocessing, the file properties are:

$ cwimage reprocessed.tif

This is different from the integer version of the images. Assumed in the above example is that the first value is the Amplitude and the second is the Phase angle in radians. These values will need to be adjusted to visualize them.

Providing Sentinel 1 data as an OGC service endpoint can be a cumbersome and not a straightforward process as data is provided in a format that is not supported by any of the OGC standards.

As a default, Sentinel-1 SAR data is not provided as GeoTIFF. To setup an OGC WMS or WCS endpoint or something similar, the data needs to be transformed into the GeoTIFF format. To accomplish this, processing steps are needed to generate a GeoTIFF, including orthorectification of the data. The SAR data provider should consider providing pre-processed data for ML workflows preferable for scenes that are used to create labeled data such that the same digital elevation model is used.

Readers are also invited to review the discussions in the D015 Machine Learning ER surrounding the use of OGC APIs in this task.