A Systematic Approach to Identify Shipping Emissions Using Spatio-Temporally Resolved TROPOMI Data

Abstract

Stringent global regulations aim to reduce nitrogen dioxide (NO₂) emissions from maritime shipping. However, the lack of a global monitoring system makes compliance verification challenging. To address this issue, we propose a systematic approach to monitor shipping emissions using unsupervised clustering techniques on spatio-temporal georeferenced data, specifically NO₂ measurements obtained from the TROPOspheric Monitoring Instrument (TROPOMI) on board the Copernicus Sentinel-5 Precursor satellite. Our method involves partitioning spatio-temporally resolved measurements based on the similarity of NO₂ column levels. We demonstrate the reproducibility of our approach through rigorous testing and validation using data collected from multiple regions and time periods. Our approach improves the spatial correlation coefficients between NO₂ column clusters and shipping traffic frequency. Additionally, we identify a temporal correlation between NO₂ column levels along shipping routes and the global container throughput index. We expect that our approach may serve as a prototype for a tool to identify anthropogenic maritime emissions, distinguishing them from background sources.

1. Introduction

Nitrogen dioxide (NO₂) is a highly reactive trace gas that plays a crucial role in both public health problems and environmental pollution. It is primarily released into the air through fuel combustion [1,2,3]. The presence of NO₂ in the atmosphere exacerbates respiratory symptoms in humans and causes environmental pollution through acid rain [3,4,5]. Shipping is a significant contributor to global air pollution, responsible for about 15% of global NO_x emissions [6,7], particularly in areas near coastlines [8]. As a response, governments and non-governmental organizations are striving to reduce anthropogenic NO₂ emissions from ships. The International Maritime Organization (IMO) has enforced regulations on the emission of nitrogen oxides (NO_x) from newly built vessels to comply with stricter standards [9]. Ships must comply with NO_x emission regulations when entering NO_x emission control areas, and stricter regulations are imposed on newer ships to further reduce marine NO_x emissions as shown in Figure 1.

Governments and international organizations have urgently called for a reduction in shipping emissions [10,11,12]. To achieve this reduction, monitoring of shipping emissions is essential. Recent studies have attempted to track NO₂ ship emissions through remote sensing with the advent of technology such as the TROPOspheric Monitoring Instrument (TROPOMI) on board the Copernicus Sentinel-5 Precursor satellite, which provides high sensitivity and smaller ground sampling [13]. One study has shown the possibility of visually tracing ships based on their exhaust as recorded by TROPOMI [14]. However, this approach requires continual human effort to review numerous remote sensing images daily, making it susceptible to human error and subjectivity. Instead, selectively investigating a region identified as suspicious of particularly high structural bias in terms of NO₂ concentration can be helpful for effective regulation implementation.

To address this issue, we propose a systematic approach to separate shipping signatures from background values in the NO₂ columns. Our approach includes spatio-temporal data preprocessing and unsupervised clustering of NO₂ columns, enhancing the interpretation and visualization of observational TROPOMI data for busy shipping routes. The purpose of this paper is to pinpoint relative hot spots and trends of NO₂ concentration in the shipping route, not the absolute NO₂ emission estimation. We conducted our research during the period from 5 December 2018 to 10 November 2021 on three regions of interest with considerably higher NO₂ emissions [15]: the Mediterranean Sea, the Red Sea, and the Indian Ocean. We verified the validity of our approach by comparing the results with publicly available ship track count estimates [16] and a real-world global logistics index [17] to ensure that the clustering result reflects the shipping patterns and economic trends. Our results can be used to determine whether the amount of NO_x emissions is consistent with the expected level from all the ships in a given shipping route or whether there is a structural bias in the NO_x emissions.

In Section 2, we present a comprehensive review of related works on tracking NO₂ emissions from ships. The current state of the dataset used in this paper is checked in Section 3. In Section 4, we propose the research model and its structure. The methodology of the research model is explored in detail in Section 5. Section 6 presents the results obtained from the methodology applied to the dataset and validates the results by collating them with a real-world global logistics index. Finally, we discuss our findings in Section 7.

2. Related Works

Researchers have attempted to measure ship NO_x emissions by directly analyzing ship plumes using mobile and/or stationary instruments. Berg et al. [18] conducted flight transects across ship plumes to verify the feasibility of airborne optical measurements using reflected skylight from the water surface as the light source. Alföldy et al. [19] established an empirical relationship between engine operating conditions and NO₂ emission factors from ship plumes at the entry to the port of Rotterdam. Pirjola et al. [20] observed a significant decrease in NO_x emissions from ships equipped with selective catalytic reduction, with the help of a mobile laboratory van. Lööv et al. [21] measured NO₂ emissions from ships by extracting ship plumes using a sonde instrument, suggesting that optical remote sensing of NO₂ emissions from ships can provide reliable estimates compared to direct measurements.

After the successful verification of optical remote sensing technology in ship emission measurements, further studies have established a connection between satellite observational data and ship emissions. Richter et al. [22] visualized NO₂ emissions along busy shipping routes, although those observations with poor satellite resolution may be dominated by an aged ship plume. Boersma et al. [23] discovered a noticeable drop in average NO_x emissions in Europe via satellite measurements achieved by the implementation of seafaring ships that operate at slow speed. Georgoulias et al. [14] proposed that one can detect individual ship emissions from the TROPOMI NO₂ column density by comparing them with automatic identification system (AIS) data. Riess et al. [24] visualized the correlation of TROPOMI NO₂ column density with shipping emissions over European shipping lanes, displaying the impact of coronavirus disease 2019 (COVID-19) lockdown measures observed in AIS data.

Despite having the best resolution, TROPOMI measurements are not exempt from spatial dependence. The first law of geography states that “everything is related to everything else, but near things are more related than distant things” [25]. In other words, georeferenced observations tend to have a high correlation with their neighboring observations, known as spatial autocorrelation or spatial heterogeneity. Generally, a spatial dataset is not completely independently distributed but rather displays information redundancy [26]. Georeferenced remote sensing data can have large relative uncertainties in individual measurements due to spatial variations in NO₂ profiles and surface albedo [27]. Therefore, reducing information redundancy before satellite data analysis is important.

To the best of our knowledge, most of the previous analyses using TROPOMI measurements have not properly handled spatial autocorrelation. Most researchers deal with spatial data by oversampling, which is a technique of collecting numerous data from the same geographic coordinates over time and then averaging them. This technique can also be referred to as the process of averaging Level 2 data over a long period to obtain a Level 3 grid that is significantly finer than the Level 2 pixels [28]. Several studies have shown the effectiveness of oversampling in enhancing resolution to resolve complex spatial heterogeneity [29,30,31,32,33]. This may be convenient in creating representative values for given geographic coordinates, but temporal oversampling may have averaged out the daily spatial heterogeneity [34]. Taking into account the first law of geography, a single daily event of NO₂ emission exerts spatial influence on the NO₂ concentration level of neighbors on the same day. In other words, the local spatial autocorrelation is the main driver of spatial heterogeneity on daily TROPOMI measurement. Recently, a study demonstrated an image processing technique that enhances the TROPOMI image quality by reducing background noise using spatial association statistics [35,36]. However, this method has not yet been generalized to other regions and multiple time periods to reflect the real-world business cycle.

In general, machine learning models require well-prepared data for accurate predictions [37]. Therefore, data preprocessing can significantly impact the performance of a predictive model [38]. In this paper, we demonstrate that our methodology can properly identify shipping routes when the processed spatio-temporal TROPOMI data are utilized. The main contributions of our research are as follows:

• We demonstrate an unsupervised method to identify shipping emissions using NO₂ column density data without the true ship location information (e.g., AIS data) provided prior to the analysis.
• We propose guidelines for reflecting spatial and temporal heterogeneity when analyzing shipping emissions using remote sensing data.
• We demonstrate the efficacy of our method in identifying shipping emissions by validating it with real-world data of spatial ship track location and temporal economic trend data.

3. Research Scheme

In this section, we provide an overview of the systematic approach to identify shipping emission activities. To aid comprehension, we illustrate the proposed approach as a schematic representation in Figure 2.

First, we propose using NO₂ column density as a proxy for shipping emissions in the open sea. Since NO₂ is the dominant form of NO_x emitted from aged ship plumes in atmospheric air [22,39,40], we can consider the NO₂ concentration in ship plumes as a proxy for the total amount of ship emissions. Provided there is no stationary anthropogenic emission stronger than shipping, which we assume to be the case, the distinctively high NO₂ column density in the open sea primarily results from cumulative shipping emission activities, namely the cumulative NO₂ concentration along the shipping route. Consequently, we suggest using the average level of cumulative NO₂ column density as a proxy for shipping frequency.

Second, we present a method for handling spatial autocorrelation. As TROPOMI data are inherently georeferenced data, we must deal with spatial dependence for each daily image. Therefore, we need to reduce the interdependence among spatial data, firstly by removing outflow from the land near the coastline, and then quantifying the spatial autocorrelation of the NO₂ column density.

Third, we demonstrate how to partition the NO₂ column density by k-means clustering of spatially processed and temporally oversampled data. The clustering yields labels that the user can annotate with interpretations, such as the shipping route with the highest average NO₂ column density.

Finally, we validate the outcome by visualizing spatially and quantifying how TROPOMI NO₂-based clusters correlate with ship tracks from a publicly available cumulative ship track dataset [16]. We further validate the result by plotting the time series of TROPOMI NO₂ column density and the global economic index to show how the NO₂ column density trend for shipping routes is similar to the real-world container throughput trend [17].

4. Materials

4.1. TROPOMI NO₂ Datasets

In this study, we demonstrate a systematic approach to identify shipping emission activities using remote sensing instruments. Specifically, we utilize the remotely sensed NO₂ column density retrieved by TROPOMI, a satellite instrument that has gained attention due to its unprecedented high resolution and sensitivity.

Compared to its predecessors like Global Ozone Monitoring Experiment [41] and Ozone Monitoring Instrument [42], TROPOMI on board the Sentinel-5 Precursor satellite [13] has a much higher sensitivity and smaller ground sampling technology. With a resolution of $3.5\times 7.0$ ${\mathrm{km}}^2$ near nadir, which has been improved to $3.5\times 5.5$ ${\mathrm{km}}^2$ as of August 2019 [43], and can be further enhanced to $1.0\times 1.0$ ${\mathrm{km}}^2$ by oversampling [29,30,31,32,33], TROPOMI can be used for mapping and tracking shipping activities [14].

From the TROPOMI NO₂ product, we make use of two datasets: (1) the NO₂ tropospheric vertical column density (VCD) and (2) the NO₂ tropospheric slant column density (SCD), retrieved using differential optical absorption spectroscopy [44,45,46]. Specifically, we have used the near real-time (NRTI) Level 2 TROPOMI data and then processed this with our own in-house Level 3 processor. The data are stored in the georeferenced Network Common Data Form format [47], with four dimensions: (1) latitude (unit: decimal degree); (2) longitude (unit: decimal degree); (3) time (unit: day); and (4) NO₂ column density (unit: mol m${}^{-2}$).

To ensure the reliability of the data, we consult the TROPOMI Algorithm Theoretical Basis Document [43] and select NO₂ column density data based on the following criteria: quality assurance value (qa_value) $>0.5$, cloud fraction $<0.2$, and viewing zenith angle $<60$°. Regarding the viewing zenith angle, the value of $<60$° is arbitrary, but this should be enough if the purpose of the filter is to filter our poorer quality data. These criteria were recommended at the time of data retrieval to include good-quality data and moderately filter out data with errors in retrievals, to gather data from a clear sky, and ensure data reliability, respectively. Regarding the qa_value in particular, we used qa_value $>0.5$ as it is a useful criterion where the averaging kernels are used [48]. The proposed method suffices this condition as we implement spatial processing of NO₂ column density data.

It is worth noting that the purpose of this paper is to verify the validity of the statistical analysis on trends and changes of the NO₂ concentrations rather than to estimate the absolute amount of emission. Therefore, we choose to include lower-quality data to demonstrate the validity and robustness of the proposed method against such data. After selecting the quality of measurement, we regrid and average the data. Specifically, we interpolate all irregularly gridded Level 2 data on a regular lattice defined in the World Geodetic System 1984 (WGS84) and time-average the data for the same measurement pixel available with valid data from different points in time on the same day.

To promote the generalizability of the proposed method, we analyzed TROPOMI NO₂ column density data from multiple regions over an extended period.

Table 1. Spatio-temporal coverage of the six different TROPOMI NO₂ datasets used.

Dataset	Latitudinal Coverage	Longitudinal Coverage	Temporal Coverage
Mediterranean Sea—VCD	33.015625${}^{°}$N–44.984375${}^{°}$N	11.98438${}^{°}$W–15.98438${}^{°}$E
Mediterranean Sea—SCD	(384 pixels)	(896 pixels)
Red Sea—VCD	12.515625${}^{°}$N–22.484375${}^{°}$N	37.51563${}^{°}$E–47.48438${}^{°}$E	5 Dec 2018–10 Nov 2021
Red Sea—SCD	(320 pixels)	(320 pixels)	(1072 days)
Indian Ocean—VCD	0.015625${}^{°}$N–14.984375${}^{°}$N	75.01563${}^{°}$E–104.98438${}^{°}$E
Indian Ocean—SCD	(480 pixels)	(960 pixels)

shows the properties of the datasets used in the paper, including the spatial and temporal coverage of each dataset.

In terms of spatial coverage, we selected three of the busiest shipping routes in the world, as illustrated in Figure 3. The dataset covers the following geographic regions with bounding WGS84 coordinates in decimal degrees as shown in Table 1. The equirectangular-regridded NO₂ column density data points have the resolution of a square pixel with a sampling rate of 0.03125 degrees in both latitude and longitude. Thus, a pixel in the regridded map represents the NO₂ column density within a quadrangle of approximately $3.5\times 5.5$ ${\mathrm{km}}^2$ at the equator, where the Earth’s circumference is approximately 40,075 km long. The pixel width of the regridded map decreases as the latitude of the pixel moves away from the equator. Consequently, there is a slight oversampling in the regridded map as the latitude moves away from the equator.

In terms of temporal coverage, we selected the time period covering before and after the outbreak of the COVID-19 pandemic as shown in Table 1. For each day, there is a single NO₂ column density per pixel, derived from the TROPOMI measurement at the 13:30 overpass local solar time. In other words, the dataset has a temporal coverage of a total of 1072 days with one measurement per day.

Formally speaking, the regridded TROPOMI NO₂ column density datasets serve as the original input data X defined in a cubic lattice space of latitude, longitude, and time. We define $X\left(t\right)$ as a square lattice that consists of all original input data values at time t. Within $X\left(t\right)$, we define a single entry of $X\left(t\right)$ as $x_i\left(t\right)\in \mathbb{R}$, referenced over coordinates $i\in {\mathbb{R}}^2$ at time $t\in \mathbb{R}$.

4.2. Spatio-Temporal Validation

To validate the proposed method for identifying shipping emissions, we imported two open-source datasets:

• Spatially, we brought in a ship track count data [16] to verify that our results match real-world ship track data. Each geographic coordinate has vessel location information representing the total counts of ship tracks recorded in the data collection period. We used the normalized form of the total ship track counts per region. Lastly, the normalized counts were regridded on the same coordinates as the TROPOMI NO₂ column density.
• Temporally, we brought in a global container throughput index [17] to demonstrate that our results reflect the global business cycle of shipping volume. This index is an estimate of global container throughput, representing about 60% of the global trade volume via shipment collected from a total of 89 international ports from December 2018 to November 2021. Since our regions of interest lie along busy shipping routes connecting Europe and Asia, the trend of the global trade volume should be reflected in the trend of emission concentrations in these regions.

In general, the spatial validation dataset serves as the true value of the normalized ship track counts $y_i^{spat}\in {\mathbb{R}}^2$ considered over coordinates i. The temporal validation dataset serves as the true value of the global container throughput index $y_i^{temp}\in {\mathbb{R}}^2$ considered over coordinates i.

5. Methods

5.1. Data Preprocessing

5.1.1. Near-Coast Removal

We begin by removing spatially irrelevant NO₂ data from inland emission sources. The NO₂ column densities show relatively large uncertainties near the coastline and localized sources such as power plants [27]. Therefore, we need to remove not only the NO₂ data from inland emission sources but also from those near the coastline, as these are most likely to be outlier values.

Inland data must be removed since the NO₂ signature will be dominated by stationary inland emission sources. To remove data from irrelevant inland emission sources, we first mask the data over land. This can be accomplished by masking with the freely available Natural Earth 1:10 m physical vector map of land polygons, including major islands [49].

Since NO₂ travels in the atmosphere via the wind, the NO₂ data over near-coast oceans are highly influenced by spillover from inland emission sources and must be treated as outliers. Here, we presume the width of the near-coast oceanic region to be about 70 km away from the coastline. The near-coast ocean region is removed using morphological binary erosion [50]. Binary erosion can be defined as the binary morphological operation (⊖) between the original image of an image in a Euclidean space E, operand $A\subset E$, and the shape parameter or the kernel, operand $B\subset E$, such that

$$A\ominus B=\{z\in E|B_z\subseteq A\}$$

(1)

where $z\in E$ denotes a translation vector of B to probe across A, and $B_z$ is a translation of B to z. If $B_z$ can be contained in A, then z is said to belong to the erosion $A\ominus B$. z has a value of one when it belongs to the erosion $A\ominus B$, and zero otherwise. A visual example of the binary erosion operation is shown in Figure A2. To erode the NO₂ column density near the coast, we use a binary Gauss circular kernel, illustrated in Figure A3, with a radius of twenty pixels, or equivalently, a geodetic distance of 70 km. In other words, k and l are set to twenty for the erosion of near-coast pixels. Note that the choice of k and l are arbitrary, and they are fixed to twenty to promote the simplicity of the proposed method. A visual example of near-coast removal is shown in Figure A4. By removing data near the coast, most of the data with large NO₂ column density with almost no ship track are removed, as displayed in Figure A5. After the removal of near-coast pixels, we obtain a NO₂ column density that corresponds better with open-sea ship tracks.

5.1.2. Spatial Data Preprocessing

Most of the time, spatial data contain redundant information [26]. In other words, geographical data do not typically represent independently generated incidents, but instead reflect a sum of spatial spillover influences exerted by neighboring data [51,52]. For example, remotely sensed gas leakage images [53] show that pixels near the emission source tend to have higher readings than those farther away. As mentioned earlier, NO₂ emission events have a spatial spillover effect on the NO₂ levels of neighboring areas. This effect is much larger for local pockets of spatial dependence [54,55] near the anthropogenic NO₂ source [56,57].

Traditional regression models, such as the least squares method, assume independently distributed data, which may not hold true for spatial data [51]. Additionally, the least squares method produces distorted results with outliers, as they dominate the sum of squares calculation [58]. Thus, to resolve the information redundancy of spatial data, we quantify the original NO₂ column density into a local indicator of spatial autocorrelation (LISA). Two popular LISA metrics are Getis-Ord $G_i^{*}$ [54,59] and local Moran’s I [55].

In this study, we use the $G_i^{*}$ statistic because we aim to identify the hot and cold spots of NO₂ concentration to identify the shipping signature. Hot (or cold) spots are the specific spatial clusters of high (or low) values, with zero being the division criterion. Although local Moran’s I is a popular LISA metric that can quantify the statistical significance of spatial heterogeneity, it only represents the significance of the spatial heterogeneity and not whether it is a hot or cold spot. Using local Moran’s I may assign high values to both shipping and the background signature if both are significantly spatially heterogeneous. In contrast, we want to clearly distinguish the shipping signature from the background and not categorize them together with the same statistical significance. Therefore, we select $G_i^{*}$ as the spatial data preprocessing tool to generate both hot and cold spots.

$G_i^{*}$ was originally introduced by Getis and Ord [54] as a statistic quantifying the spatial autocorrelation of a local pocket near a target data point. First, we define the target data point as $x_i\left(t\right)$ with known coordinates i at time t, using the same notation as in Section 4. Then, $G_i^{*}$ measured at $x_i\left(t\right)$ is defined as follows:

$$G_i^{*}\left(x_i\left(t\right)\right)=\frac{{\sum }_j{w}_{ij}\left(d\right)x_j\left(t\right)}{{\sum }_j{x}_j\left(t\right)}$$

(2)

where $w_{ij}\left(d\right)$ denotes the symmetric binary spatial weight matrix between $x_i\left(t\right)$ and $x_j\left(t\right)$, whose elements are one when $x_j\left(t\right)$ is within distance d from data $x_i\left(t\right)$ and zero otherwise (j may equal i). Typically, d is chosen as a squared distance in Euclidean space.

Later in [59], the statistic is revised into a standardized statistic, denoted here as $Std{G}_i^{*}$. $Std{G}_i^{*}$ measured at $x_i\left(t\right)$ is defined as follows:

$$\begin{array}{cc}\hfill Std{G}_i^{*}\left(x_i\left(t\right)\right)& =\frac{{\sum }_j{w}_{ij}\left(d\right)x_j\left(t\right)-\overline{x}\left(t\right)}{s\sqrt{\frac{N{\sum }_j{w}_{ij}{\left(d\right)}^2-{\left({\sum }_j{w}_{ij}\left(d\right)\right)}^2}{N-1}}}\hfill \\ \hfill \mathrm{subject}\mathrm{to}\overline{x}\left(t\right)& =\frac{{\sum }_j{x}_j\left(t\right)}{N}\hfill \\ \hfill s& =\sqrt{\frac{{\sum }_j{x}_j{\left(t\right)}^2}{N}-{\left(\overline{x}\left(t\right)\right)}^2}\hfill \end{array}$$

(3)

where $\overline{x}\left(t\right)$ denotes the sample mean of all measurements in $X\left(t\right)$, s denotes the sample standard deviation of all measurements in $X\left(t\right)$, and N denotes the sample size.

A drawback of Equation (3) is the expensive computation of the weight matrix, especially for remote sensing data. As remote sensing data generally have a large number of data entries, the observations are usually non-zero, unless filtered out by a low quality-assurance value or a high cloud fraction. In conventional spatial data analytics tools such as ArcGIS [60], the spatial weight matrix $w_{ij}\left(d\right)$ for every target region i is constructed for all neighbors within d. Although it utilizes a sparse spatial methodology to save storage space and reduce execution time [61], selecting a high number of neighbors can cause memory issues [60].

Therefore, we modify the calculation of the spatial weight matrix with convolution filtering [62], one of the most widely used image processing algorithms in pattern recognition and computer vision. By replacing the weight matrix $w_{ij}\left(d\right)$ with its convolution, we obtain the same result without the expensive computation of the spatial weight matrix. The finite discrete convolution operation $K*I\in {\mathbb{R}}^2$ between the kernel matrix $K\in {\mathbb{R}}^2$ and the target intensity matrix $I\in {\mathbb{R}}^2$ is defined as follows:

$$\begin{array}{cc}\hfill K\ast I=& \left(\begin{array}{ccc}K[-k,-l]& \cdots & K[-k,l]\\ ⋮& \ddots & ⋮\\ K[k,-l]& \cdots & K[k,l]\end{array}\right)\ast \left(\begin{array}{ccc}I[1,1]& \cdots & I[1,n]\\ ⋮& \ddots & ⋮\\ I[m,1]& \cdots & I[m,n]\end{array}\right)\hfill \\ \hfill =& \sum _{dk=-k}^k\sum _{dl=-l}^{l}K[dk,dl]I[\alpha -dk,\beta -dl]\hfill \end{array}$$

(4)

where k denotes the vertical lattice tick, l denotes the horizontal lattice tick, and $\alpha$ and $\beta$ are the coordinates of I for $\forall \alpha \in [1,m]$ and $\forall \beta \in [1,n]$. We refer to the element of $K*I$ at coordinates i as ${\left(K*I\right)}_i$. When the convolution operation encounters an empty cell from I, such as the cell beyond the boundary of I or values filtered out for quality, the operation regards such an entry as empty and calculates using other available data. If none of the entries are available, then the operation returns the target data point as empty.

Equation (4) represents the mathematical model for performing convolution operations between the kernel matrix and target intensity matrix. The kernel matrix K can be designed to detect specific features in an image, such as edges, corners, or other relevant patterns. The convolution operation convolves the target intensity matrix I with the kernel matrix K to produce an output matrix $K*I$, which highlights the detected features in the image.

Using the aforementioned convolution operation, we propose a modified version of the $Std{G}_i^{*}$ statistic, denoted as $Conv{G}_i^{*}$, which measures the spatial autocorrelation of the local pocket near a target data using convolution filtering. The $Conv{G}_i^{*}$ at coordinates i between the kernel K and input intensity dataset at time t, $X\left(t\right)$, is defined as follows:

$$Conv{G}_i^{*}\left(x_i\left(t\right)\right)=\frac{{\left(K\ast X\left(t\right)\right)}_i-K\ast \left(\overline{x}\left(t\right)J_{m,n}\right)}{s\sqrt{\frac{N\left(J_{1,k}{K}^2{J}_{l,1}\right)-{\left(J_{1,k}K{J}_{l,1}\right)}^2}{\left(N-1\right)}}}$$

(5)

where $J_{m,n}$, $J_{1,k}$, and $J_{l,1}$ denote matrices of size $m\times n$, $1\times k$, and $l\times 1$, respectively. The kernel size parameters are generally selected as $k<m$ and $l<n$ to compute the local spatial autocorrelation. While arbitrary size is allowed for the kernel, we use a square kernel matrix of size $2k+1\times 2l+1$ for convenience. Note that if even numbers of rows and columns are selected, then the kernel size becomes $2k\times 2l$. For $Conv{G}_i^{*}$ computation, we use a binary Gauss circular kernel with a radius of 5 pixels, or equivalently, 17 km geodetic distance. In other words, k and l are set as five for the computation of spatial autocorrelation using convolution. Note that the choice of k and l are arbitrary, and they are fixed to five to promote the simplicity of the proposed method. By such computation, we obtain the spatially resolved $Conv{G}_i^{*}$ statistic of the NO₂ column density, as visualized in Figure A6.

5.1.3. Temporal Data Preprocessing

After resolving the spatial autocorrelation, we can proceed to the next step of temporal data preprocessing to further enhance the resolution. As oversampling is a proven technique to achieve finer resolution [29,30,31,32,33], we can average the NO₂ column density over time. Several previous studies have conducted research on the monthly averaged TROPOMI NO₂ column density. However, we strongly believe that the location of a busy shipping route does not vary month by month, because the shortest shipping path, which is the most economic route [63], is already determined and is likely to remain constant for years. Therefore, we extend the oversampling window over the entire test period as shown in Table 1 to obtain a single shipping route that reflects the shipping patterns over the entire period. Conventional oversampling is equivalent to taking the sample mean of the NO₂ column density at coordinates i over time $t\in [1,T]$.

Instead of conventional oversampling, we propose to first resolve the spatial heterogeneity. Here, we oversample from $Conv{G}_i^{*}\left(x_i\left(t\right)\right)$ instead of the original input data:

$${\overline{x}}_i=\frac{{\sum }_{t}Conv{G}_i^{*}\left(x_i\left(t\right)\right)}{T}$$

(6)

As mentioned earlier, oversampling may exclude daily events with overshooting NO₂ column density. While oversampling is effective in suppressing noise signatures [28], it may underestimate unexpectedly high peaks of NO₂ emission from non-compliant ships. Despite this drawback, oversampling retains its spatial significance throughout the test period as it captures constant ship activities around the busy shipping route. Therefore, the oversampling technique does not compromise our aim of retaining shipping route information from the TROPOMI NO₂ column density.

5.2. Clustering

Partitioning can aid human interpretation by assigning labels to each pixel while considering the column density values of neighboring pixels. As TROPOMI datasets contain continuous numeric data, it is important to distinguish which values are relatively high in comparison to other values in the region of interest. We partition the spatio-temporally resolved TROPOMI datasets using a clustering algorithm to obtain cluster labels for the NO₂ column density. Although the approach is not dependent on any specific clustering algorithm, we demonstrate the approach using the k-means clustering algorithm, an unsupervised classification of grouping patterns into k clusters. The original concept of k-means clustering was introduced by Steinhaus [64], and the full algorithm was developed by researchers such as Lloyd [65] and MacQueen [66] in the mid-20th century. Despite its age, k-means clustering remains one of the most commonly used clustering techniques due to its simplicity and speed [67]. We adopt the k-means clustering algorithm suggested by Lloyd [65] to partition the spatio-temporally resolved NO₂ column density using squared distance as the dissimilarity measure. We define the k-means clustering algorithm as an optimization problem of finding k centroids that minimize the within-cluster sum of squares (WCSS) as follows:

$$\begin{array}{cc}\hfill min& \sum _{i=1}^n\sum _{j=1}^k{a}_{ij}|{\overline{x}}_i-{\mu }_j|\hfill \\ \hfill \mathrm{subject}\mathrm{to}& {\mu }_j=\frac{{\sum }_{i=1}^n{a}_{ij}{\overline{x}}_i}{{\sum }_{i=1}^n{a}_{ij}},\hfill \\ \hfill & \sum _{j=1}^k{a}_{ij}=1,\hfill \\ \hfill & \sum _{i=1}^n{a}_{ij}\ge 1\hfill \\ \hfill \mathrm{for}& i\in [1,n],j\in [1,k]\hfill \end{array}$$

(7)

where ${\mu }_j$ denotes the centroid of the jth cluster $C_j$, and $a_{ij}$ denotes a weight of value one if ${\overline{x}}_i$ belongs to $Cj$ and zero otherwise.

One caveat of using the k-means clustering algorithm is that the centroids of produced clusters are only meaningful when the cluster shape is roughly globular [68]. In our case, we aim to identify shipping emissions, where ships travel along linear routes and emit exhaust while sailing. As a result, the ship exhaust detected by remote sensing tends to follow a linear path and form an elongated shape, rather than a globular one. Therefore, we cannot simply cluster data that exhibit spatial spillover effects.

To properly handle spatially autocorrelated data such as TROPOMI images, we propose quantizing the NO₂ column density as a scalar by collapsing the spatial dimension. This is because the collapsed data no longer depend on the spatial parametric setting, and scalar quantized data can be georeferenced back to their original coordinates. Thus, clustering will preserve the non-parametric spatial structure, such as an elongated shape, even after clustering. We use Lloyd’s algorithm in its original format of scalar quantization, where georeference information is collapsed and all data are sorted in terms of the NO₂ column density level. Our objective is to find a scalar quantizer that partitions the data into k levels. Scalar quantization can be performed by applying k-means clustering to the one-dimensional NO₂ column density.

Another caveat is the dependency of the algorithm on the initial assignment of centroids. This issue is often mitigated by performing multiple trials of k-means clustering and optimizing the results by selecting the best outcome. To this end, we employed k-means++ [69], an improved initialization technique with randomized seeding.

A final caveat is that the performance of k-means clustering highly depends on a predetermined number of clusters. As for the dependency on the predetermined cluster number, it can be solved analytically by choosing an optimal number of clusters. While there are numerous techniques for determining the optimal number of clusters, we demonstrate one of the most common methods: the elbow method [70]. This method chooses the optimal number of clusters where the elbow of the WCSS curve is located. We select the elbow of the WCSS curve of cluster k, denoted as $D_k$, by performing k-means clustering over a set range of k. Formally, we define the distortion as follows:

$$\begin{array}{cc}\hfill D_k& =\sum _{{obs}_i\in C_k}\sum _{{obs}_j\in C_k}|{obs}_i-{obs}_j|\hfill \\ & \simeq 2n_k\sum _{{obs}_i\in C_k}|{obs}_i-{\mu }_k|\hfill \end{array}$$

(8)

where any pair of observations ${obs}_i$ and ${obs}_j$ belongs to the same cluster $C_k$ with $n_k$ observations, and ${\mu }_k$ denotes the centroid of $C_k$. The elbow point is defined as the point where the k-means clustering achieves the best fit or where the variance is explained as much as possible before overfitting occurs. We assume that the elbow point is located at the largest difference between the original curve and its reference line, as shown in Figure A7. In other words, it is a maximization problem of the WCSS difference given a set of the number of clusters. We calculate the optimal number of clusters using the elbow method, which is five for the Mediterranean Sea, four for the Red Sea, and six for the Indian Ocean. The cluster labels are sorted in ascending order of the spatially resolved centroid of NO₂ column density using $Conv{G}_i^{*}$. Finally, the shipping signature can be identified as the cluster with the highest NO₂ column density on average, by integrating spatio-temporal data preprocessing and clustering techniques.

5.3. Quantification of Spatial and Temporal Correlation

For the quantification of spatial and temporal data preprocessing, we use the sample Pearson correlation coefficient r as a model assessment criterion. The bivariate sample Pearson correlation between two data points u and v is defined as follows:

$$r_{u,v}=\frac{{\sum }_i\left(u_i-\overline{u}\right)\left(v_i-\overline{v}\right)}{\sqrt{{\sum }_i{\left(u_i-\overline{u}\right)}^2}\sqrt{{\sum }_i{\left(v_i-\overline{v}\right)}^2}}$$

(9)

This Pearson correlation is used to calculate either spatial correlation or temporal correlation between two data.

6. Results and Discussion

6.1. Validation of Spatial Data Preprocessing

In this section, we validate the proposed methodology by demonstrating the increased correspondence of spatio-temporally preprocessed NO₂ column density with publicly available shipping frequency patterns and economic trends, compared to the original column density and clustered column density without spatio-temporal data preprocessing.

First, we calculate the spatial correlation between the clustering results of the NO₂ column density and the normalized ship track counts. The spatial correlation is calculated by solving Equation (9) between the NO₂ column data, such as the original or clustered data, and the validation data, which are the normalized ship track counts. We assume that the true shipping frequency can be estimated by the ship track counts dataset [16]. Such an assumption is based on a strong belief that the location of the shipping route remains constant for years due to the ship operators’ tendency to travel along the shortest path, as explained in Section 5.1.3.

We will now compare the two clustering results: one before resolving spatial autocorrelation by spatial data preprocessing (Figure 4a–c), and the other after resolving it (Figure 4d–f). The cluster label is a sorted representation of the average NO₂ column level within the cluster. As explained in Section 5.2, the cluster labels are sorted in terms of the NO₂ column level of cluster centroid, for instance yielding from 1 to 5 in the Mediterranean Sea. Each cluster label may have its own interpretation, and we assign meanings to the highest and the lowest as the cluster with the highest NO₂ concentration and the background cluster or the reference cluster. Since our proposed method aims towards the systematic identification of shipping emissions, it will be easier if we could assume that the cluster with the highest NO₂ concentration corresponds to the shipping route cluster.

It is non-trivial to manually assign the shipping route label to any cluster in between, since the number of clusters is subject to change per region and such interpretation is prone to human subjectivity. Before preprocessing, the cluster labels with higher labels correspond to the coordinates with relatively high average NO₂ emission concentration levels compared to other clusters. Since it is not spatially preprocessed, the highest may not correspond to the shipping route, but rather the NO₂ emissions from the land. Such behavior is shown in the boxplots of Figure 5a–c, as the clusters with the highest labels do not always have the highest ship track counts. It is therefore a difficult task before preprocessing to correctly interpret which cluster label is the shipping route.

We propose to overcome this non-triviality in shipping route identification by spatially preprocessing the NO₂ concentration data and relying on the systematic approach to identify the shipping route as the highest cluster label. After spatial data preprocessing, higher cluster labels correspond to higher ship track counts, indicating that the preprocessed clusters better reflect the shipping density pattern. Spatial data preprocessing particularly enhances the meaning of the highest cluster label in the Red Sea and in the Indian Ocean, as the clusters before preprocessing may have been largely influenced by spatial spillover and inland emission activities. The changes from spatial preprocessing are apparent, as the cluster shape of the highest cluster labels visually correspond to the major busy shipping routes in test regions of Figure 4.

Figure 5 demonstrates the effectiveness of the $Conv{G}_i^{*}$ operation by showing the proportional relationship between cluster labels and ship track ratios. The results show that spatial data preprocessing improves the correspondence of NO₂ column density to the shipping pattern, as it preserves the increasing trend of normalized ship track counts. We intend to cluster the NO₂ column density in proportion to ship track counts, so the cluster label with the highest number should represent the largest amount of ship track counts. In particular, the cluster labels of the Red Sea recover their monotonicity after spatial preprocessing, as near-coast removal excludes the land-based emission influence from the coastal regions.

The effectiveness of spatial data preprocessing lies in the design of the k-means algorithm. Since the algorithm calculates WCSS generally based on a distance metric, the k-means clustering algorithm is sensitive to outliers. By removing the information redundancy and the outflow effect from coastal region activities via spatial data preprocessing, we see an overall improvement in clustering performance as shown in

Table 2. Spatial correlation between the global container throughput index and three types of NO₂ tropospheric VCD data: (1) original NO₂ tropospheric VCD data; (2) clustered NO₂ tropospheric VCD data without proposed spatial preprocessing; (3) clustered NO₂ tropospheric VCD data using proposed spatial preprocessing. The highest spatial correlation value per test region is highlighted. The p-values of correlation coefficients are shown in parentheses.

Test Region	Original NO2 Data	Clustered NO2 Data without ${\mathit{ConvG}}_{\mathit{i}}^{*}$	Clustered NO2 Data with ${\mathit{ConvG}}_{\mathit{i}}^{*}$
Mediterranean Sea	0.4479 * (0.0000)	0.4586 * (0.0000)	0.5013 * (0.0000)
Red Sea	0.1508 * (0.0000)	0.2549 * (0.0000)	0.6235 * (0.0000)
Indian Ocean	0.2357 * (0.0000)	0.3194 * (0.0000)	0.3215 * (0.0000)

p-values are in parentheses; * p < 0.01.

. As a result of spatial preprocessing, the clustered NO₂ data better represent the relatively higher NO₂ from shipping in the open sea, both visually and quantitatively.

6.2. Validation of Temporal Data Preprocessing

In addition to verifying spatial correspondence with ship tracks, we validate the proposed methodology by demonstrating how well the shipping route cluster reflects the global business cycle. Demand for shipping transportation is a key determinant of global economic activity [71,72]; therefore, the estimated shipping route should resemble the global business cycle. Researchers have found that global shipping throughput data reflect the global economic recession due to the COVID-19 lockdown measures [73,74,75]. We validate the results by comparing the monthly averages of NO₂ column density in the shipping route over the entire time period as shown in Table 1. Here, we assume that the NO₂ data point $x_i\left(t\right)$ belongs to the shipping route cluster with the highest mean NO₂ column density. To temporally validate our results, we compute the correlation between the NO₂ column density of the shipping route and the global container throughput index for the same time window. The temporal correlation is calculated using Equation (9) between the monthly averaged NO₂ column density within the shipping cluster and the monthly index of global container throughput. The global container throughput data are available online as the Rheinisch-Westfälisches Institut für Wirtschaftsforschung and the Institute of Shipping Economics and Logistics (RWI/ISL) Container Throughput Index [17].

Figure 6 shows the monthly averaged NO₂ column density of the shipping route matching the monthly trend of the global logistics volume. Its deseasonalized version is shown in Figure A8 shows the seasonally adjusted NO₂ column density, using the classical seasonal decomposition [76]. Particularly in early 2020 and early 2021, a sudden drop is observed in both the NO₂ trend and the container throughput trend, followed by recovery back to the previous level.

Table 3. Temporal correlation between the global container throughput index and two types of NO₂ column data: NO₂ tropospheric VCD of shipping route and NO₂ tropospheric SCD of shipping route. Each column density is divided into the original set and seasonally adjusted set. Deseasonalized column density is shown in Figure A8. The highest temporal correlation value per test region is highlighted. The p-values of correlation coefficients are shown in parentheses.

Region	NO2 VCD of Shipping Route		NO2 SCD of Shipping Route
	Original	Adjusted	Original	Adjusted
Medit. Sea	−0.2523 (0.1376)	−0.1857 (0.2781)	0.5567 * (0.0004)	0.4582 * (0.0049)
Red Sea	0.1992 (0.2442)	0.6824 * (0.0000)	0.6003 * (0.0001)	0.7103 * (0.0000)
Indian Ocean	−0.2244 (0.1883)	−0.1915 (0.2632)	0.5544 * (0.0005)	0.5188 * (0.0012)

p-values are in parentheses; * p < 0.01.

shows the temporal correlation between the $x_i$ and $y_i^{temp}$ time series plotted in Figure 6. The yearly differences from 2019 to 2021 are displayed in

Table A2. Average NO₂ tropospheric VCD and SCD per region, per year, and per cluster. Test regions span busy shipping routes, namely the Mediterranean Sea (Medit. Sea), the Red Sea, and the Indian Ocean. Test years span from the year 2019 to the year 2021. Clusters are distinguished as either shipping route (the cluster with the highest ship track counts) or background (the cluster with the lowest ship track counts). Each table entry represents the mean column density value within a cluster plus or minus one standard deviation, preprocessed with the $Conv{G}_i^{*}$ method.

Region	Average NO2 VCD—2019		Average NO2 VCD—2020		Average NO2 VCD—2021
	Shipping Route	Background	Shipping Route	Background	Shipping Route	Background
Medit. Sea	2.741 ± 1.627	1.208 ± 0.779	2.695 ± 1.552	1.235 ± 0.743	2.580 ± 1.461	1.166 ± 0.735
Red Sea	1.661 ± 0.948	1.160 ± 0.571	1.627 ± 0.910	1.168 ± 0.609	1.750 ± 0.898	1.273 ± 0.584
Indian Ocean	1.277 ± 0.714	0.490 ± 0.466	1.183 ± 0.748	0.442 ± 0.493	1.182 ± 0.789	0.453 ± 0.484
	Average NO2 SCD—2019		Average NO2 SCD—2020		Average NO2 SCD—2021
	Shipping Route	Background	Shipping Route	Background	Shipping Route	Background
Medit. Sea	15.942 ± 2.552	14.515 ± 2.426	14.806 ± 2.297	13.677 ± 2.097	15.974 ± 2.374	14.842 ± 2.208
Red Sea	12.140 ± 2.194	11.500 ± 1.969	11.606 ± 1.815	11.017 ± 1.663	12.652 ± 1.942	12.118 ± 1.839
Indian Ocean	9.826 ± 1.731	8.225 ± 1.259	9.425 ± 1.610	8.187 ± 1.398	10.113 ± 1.751	8.727 ± 1.675

with VCD and

Table A3. Relative change of NO₂ tropospheric VCD and SCD per region and per cluster, calculated as a ratio between years. Test regions span busy shipping routes, namely the Mediterranean Sea (Medit. Sea), the Red Sea, and the Indian Ocean. Ratios between years span from the change over 2019–2020, over 2020–2021, and over 2019–2021. In each subcolumn heading, the notation of year ratio (e.g., year A/year B) stands for the average NO₂ column density of year A divided by the average NO₂ column density of year B. For each table entry, red font indicates an increase in the ratio (mean value larger than 1), while blue font indicates a decrease in the ratio (mean value smaller than 1). Clusters are distinguished as either shipping route (the cluster with the highest ship track counts) or background (the cluster with the lowest ship track counts). Each table entry represents the mean column density value within a cluster plus or minus one standard deviation, preprocessed with the $Conv{G}_i^{*}$ method.

Region	NO2 VCD Ratio—2020/2019		NO2 VCD Ratio—2021/2020		NO2 VCD Ratio—2021/2019
	Shipping Route	Background	Shipping Route	Background	Shipping Route	Background
Medit. Sea	0.985 ± 0.047	1.029 ± 0.065	0.954 ± 0.061	0.946 ± 0.047	0.939 ± 0.062	0.972 ± 0.070
Red Sea	0.980 ± 0.038	1.008 ± 0.039	1.077 ± 0.050	1.091 ± 0.059	1.055 ± 0.050	1.099 ± 0.048
Indian Ocean	0.930 ± 0.071	0.908 ± 0.098	0.999 ± 0.081	1.023 ± 0.119	0.925 ± 0.069	0.931 ± 0.102
	NO2 SCD Ratio—2020/2019		NO2 SCD Ratio—2021/2020		NO2 SCD Ratio—2021/2019
	Shipping Route	Background	Shipping Route	Background	Shipping Route	Background
Medit. Sea	0.928 ± 0.012	0.941 ± 0.011	1.079 ± 0.012	1.084 ± 0.015	1.002 ± 0.015	1.021 ± 0.020
Red Sea	0.956 ± 0.008	0.957 ± 0.007	1.090 ± 0.010	1.099 ± 0.008	1.041 ± 0.009	1.052 ± 0.008
Indian Ocean	0.963 ± 0.025	0.996 ± 0.013	1.071 ± 0.025	1.066 ± 0.016	1.031 ± 0.016	1.061 ± 0.015

with SCD.

In addition, we have conducted another experiment using Copernicus Atmosphere Monitoring Service (CAMS) global reanalysis ECMWF Atmospheric Composition 4 (EAC4) dataset of atmospheric composition [77]. Figure A9 shows that CAMS data on the shipping route show a similar trend with NO₂ tropospheric VCD on the shipping route. Note that the trend of CAMS data matches with NO₂ tropospheric VCD, which is used for spatial preprocessing of NO₂ data and yielding clusters. This result shows that the quality of our approach to processing TROPOMI NO₂ data is comparable to CAMS global reanalysis (EAC4) result on the shipping route.

We observe that the NO₂ tropospheric SCD can better represent the global shipping density and economic trends. As shown in Table 3, the temporal correlation between the NO₂ column and the global container throughput index is much better when it is compared with NO₂ tropospheric SCD. While we propose using the NO₂ tropospheric VCD in clustering, we suggest also using the NO₂ tropospheric SCD when matching the global economic index due to the following reasons.

The TROPOMI VCD has been known to be sensitive to the a priori estimate of the NO₂ profile shape [78]. This is because the assumptions made during the a priori calculation of the air mass factor may not hold in regions with high pollution levels [79]. In particular, there is a high level of structural uncertainty in air mass factor calculation over regions with high NO₂ pollution [80].

Filtering out persistently strong background measurements, such as those caused by high solar radiation during summer and inland anthropogenic emissions, is functional and can help to rule out seasonality and overestimated inland emissions. However, providing the shipping activity locations may be appropriate for certain analyses, but not when the goal is to analyze the measurement itself directly.

Furthermore, it is worth noting that since shipping emissions are consistently high, there is a risk of underestimation. Despite underestimated NO₂ concentration in hot spots, we would like to emphasize that our approach is to identify the relative trend of the shipping emission, and not the absolute emission estimation. From the empirical result, we suspect the reason for being able to distinguish the shipping route as the tendency of low NO₂ readings remaining low, and very high readings remaining moderately high, such that it is still enough to distinguish the shipping route from the background.

These results are consistent with a previously verified global NO₂ inventory for seafaring ships, which showed relatively high NO₂ emissions during the summer period compared to the winter period ([81], Figure 4). Therefore, our empirical analysis demonstrates that the monthly trend of the NO₂ column density better matches the trend of the global container throughput index.

7. Discussion

This research proposes a method for identifying shipping NO₂ emissions by clustering spatio-temporally resolved TROPOMI data. The method partitions geographic coordinates in the open sea into clusters and identifies the cluster with the largest NO₂ column density as a shipping route. Experiments were conducted in three selected regions with high shipping density, namely the Mediterranean Sea, the Red Sea, and the Indian Ocean. Each region was spatio-temporally preprocessed and visualized on a respective map. Statistical results show that the proposed method increases the spatial correlation between the TROPOMI NO₂ measurement clusters and the ship track counts. Lastly, the TROPOMI NO₂ measurements of the highest cluster labels show temporal trends that match with the real-world global throughput index, such as the impact of the COVID-19 lockdown measures on the global trade volume. In summary, the proposed method has demonstrated its effectiveness in identifying shipping routes with high NO₂ emissions and in providing insight into the impact of global events on shipping emissions.

Since the purpose of this paper is to identify shipping emissions in the open-sea area, it may require different algorithmic settings to identify emission activities in different environments, e.g., flexible choice of kernel size to filter out near-coast measurements or calculating near-coast shipping emissions with the influence of mixed types of inland emissions. One may need to consider more supporting data when applying this method to coastal area emissions, e.g., wind velocity profiles from land and in situ measurements from coast guard authorities. In addition, the proposed method cannot effectively distinguish the clustered NO₂ data by individual vessel types, because the spatial validation is performed at an aggregate NO₂ column level.

Through this study, we identified shipping emissions by spotting the trends and changes in the NO₂ column level along the shipping route. However, it is yet unknown whether higher $Conv{G}_i^{*}$ NO₂ data are also higher in absolute NO₂ column data because the spatial preprocessing yields a quantification of showing a central tendency towards the global mean value. In other words, $Conv{G}_i^{*}$ standardizes the NO₂ column level which shows the relative position of each data point within its cluster. Therefore, it is uncertain to say that the regions with higher $Conv{G}_i^{*}$ NO₂ levels indeed exhibit higher NO₂ levels in an absolute sense.

We expect that the proposed method may serve as an exemplary technique to identify anthropogenic emission sources, distinguishing them from the background in future research. While this study focused on open-sea regions, our method can be applied to other regions with trace gas emissions, including inland anthropogenic emission sources. Depending on the purpose of the study, one may skip the land-masking procedure and/or near-coast removal to see the emissions in all available regions, which demonstrates the versatility of the proposed method. In addition, we plan to conduct optimization studies of fixed distances such as fixed geodetic distance used for $Conv{G}_i^{*}$ computation. With the help of more supporting scientific data such as the chemical transport model, we can optimize the spatial correlation based on the atmospheric domain knowledge, such as using geometric air mass factor in VCD calculation to consider its uncertainty or validating the data based on comparisons with high-quality ground-based or other measurements [82]. With validated real-world case studies and scientifically supported results from unsupervised clustering, our method will be effective in environmental governance and regulatory compliance activities related to emission monitoring. The proposed method may be used by port authorities or inspectorates to function as a monitoring tool of a designated region to see the relative changes in aggregate emission level.