TROPOMI NO2 Tropospheric Column Data: Regridding to 1 km Grid-Resolution and Assessment of their Consistency with In Situ Surface Observations

This work presents a regridding procedure applied to the nitrogen dioxide (NO2) tropospheric column data, derived from the Copernicus Sentinel 5 Precursor Tropospheric Monitoring Instrument (S5P/TROPOMI). The regridding has been performed to provide a better comparison with punctual surface observations. It will be demonstrated that TROPOMI NO2 tropospheric column data show improved consistency with in situ surface measurements once the satellite retrievals are scaled to 1 km spatial sampling. A geostatistical technique, i.e., the ordinary kriging, has been applied to improve the spatial distribution of Level 2 TROPOMI NO2 data, which is originally sparse and uneven because of gaps introduced by clouds, to a final spatial, regular, sampling of 1 km × 1 km. The analysis has been performed for two study areas, one in the North and the other in the South of Italy, and for May 2018-April 2020, which also covers the period January 2020-April 2020 of COVID-19 diffusion over the Po Valley. The higher spatial sampling NO2 dataset indicated as Level 3 data, allowed us to explore spatial and seasonal data variability, obtaining better information on NO2 sources. In this respect, it will be shown that NO2 concentrations in March 2020 have likely decreased as a consequence of the lockdown because of COVID-19, although the far warmest winter season ever recorded over Europe in 2020 has favored a general NO2 decrease in comparison to the 2019 winter. Moreover, the comparison between NO2 concentrations related to weekdays and weekend days allowed us to show the strong correlation of NO2 emissions with traffic and industrial activities. To assess the quality and capability of TROPOMI NO2 observations, we have studied their relationship and correlation with in situ NO2 concentrations measured at air quality monitoring stations. We have found that the correlation increases when we pass from Level 2 to Level 3 data, showing the importance of regridding the satellite data. In particular, correlation coefficients of Level 3 data, which range between 0.50–0.90 have been found with higher correlation applying to urban, polluted locations and/or cities.


Introduction
Nitrogen dioxide (NO 2 ) and nitrogen oxide (NO), together indicated as nitrogen oxides (NOx = NO + NO 2 ), are significant trace gases in the Earth's atmosphere. They are emitted in the atmosphere because of natural (lightning, fires, soil microbial processes) and anthropogenic (traffic, biomass and fuel combustion) processes and/or activities. In the presence of sunlight, NO is converted into NO 2 through a photochemical process involving ozone (O 3 ) on a time scale of minutes. For this reason, NO 2 can be considered as a robust measure for concentrations of NOx [1]. NO 2 is a trace gas present in both the troposphere and stratosphere. The tropospheric NO 2 is a pollutant, which can severely affect air quality; hence, human health [2]. The gas has a strong reactivity, especially in the presence of OH, and contributes to the formation of tropospheric ozone and aerosol [3][4][5].
In the stratosphere, NO 2 is mainly produced from the oxidation of nitrous oxide (N 2 O). The stratospheric NO 2 enters important catalytic reactions, which can destroy the ozone layer [6] but can also suppress ozone depletion by converting hydrogen compounds into unreactive reservoir species such as nitric acid (HNO 3 ), that constitutes one of the significant aerosol sources. Tropospheric NO 2 is mainly present in polluted regions, where it can account for 50%-90% of the total column amount [7].
Moreover, we know (e.g., [8][9][10][11][12][13][14][15][16]) that the daily evolution of the Planetary Boundary Layer (PBL) height can affect the vertical mixing ration profile of NO 2 and its diurnal cycle. The PBL height can at most reach an altitude of 1-3 km (e.g., [11,12]); usually, it rapidly increases in the morning because of sunlight and tends to stay at or below ≈1000 m at night. Above the PBL, NO 2 concentrations approach zero (e.g., [11]). Apart from short-time scale events such as tropospheric lightning and aviation traffic, the main sources of NO 2 are at or near the ground. As a result of the evolving PBL and the intense photochemical activity of nitrogen dioxide, its concentration near the surface tends to decrease in sunlight.
For this reason, satellite observations need to discriminate between tropospheric and stratospheric nitrogen dioxide to get a better insight into detecting pollution in the lower atmosphere. In this respect, satellite observations are usually performed in the UV-VIS spectral regions to exploit the absorption mode remote sensing, which has weighting functions, which do not drop to zero at the surface. Examples of instrument sensing the near-surface in absorption mode include the Global Ozone Monitoring Experiment (GOME), onboard the European Remote Sensing 2 (ERS-2) satellite [17,18], the Scanning imaging Absorption Spectrometer for Atmosphere Chemistry (SCIAMACHY), aboard the Envisat satellite, launched in 2002 [19]. Today, tropospheric NO 2 is monitored from space with a series of instruments, such as the Ozone Monitoring Instrument or OMI, launched aboard the Earth Observing System (EOS)-Aura in 2004 [20], the GOME-2 instruments onboard EUMETSAT's Meteorological Operational Satellites (MetOp) A, B and C launched in 2006, 2012 and 2018, respectively [21] and the TROPOMI instrument, in orbit aboard the Sentinel 5 Precursor in 2017 [22].
As regards the NO 2 data obtained by the predecessors of TROPOMI, such as the OMI instrument, various downscaling techniques have been introduced to improve their spatial resolution. In [23], for example, a downscaling technique applied to OMI NO 2 vertical column density (VCD) data has been developed, using regional air quality model information and also applying averaging kernel (AK) information to adjust the vertical structure. Results have shown best agreement with ground-based observations when downscaling and AK are applied to OMI NO 2 , with a correlation coefficient of R = 0.89. In [24] an approach to improve the spatial resolution of OMI NO 2 data, by applying local scaling factors from a global three-dimensional model, has been presented. Validation results show that the correlation between OMI-derived surface NO 2 and the ground-based measurements, above all in polluted areas, is significant (R 2 0.86). Moreover, there are several works based on regression techniques. For example, in [25] a land-use regression technique to OMI NO 2 data has been applied. In particular, in situ concentration measurements and information about surrounding land-uses to estimate concentrations for non-measurement locations were used. In [26] a land-use regression technique was also used for estimating within-urban variability in air pollution at a spatial resolution of 100 m. In [27] the spatial resolution of OMI NO 2 was improved through the combined use of a geostatistical technique, the universal kriging, and a land-use regression technique. Other studies, instead, exploit the information relating to the air mass factor, a value needed to convert the slant column measurement into a vertical column amount, to obtain a high resolution of NO 2 data. Between these, in [28] a high spatial resolution model simulation was used, constrained by in situ aircraft observations, to calculate tropospheric air mass factors and to better estimate tropospheric NO 2 data. Remote Sens. 2020, 12, 2212 3 of 23 All the above studies demonstrated that downscaled NO 2 data better capture the spatial variability compared to the original datasets.
In this work, a regridding technique has been, instead, applied to the newly available Sentinel-5P/TROPOMI NO 2 tropospheric column data. In particular, month averaged NO 2 data have been re-gridded from the TROPOMI spatial resolution (7 km × 3.5 km at the sub-satellite point) to the finer resolution of 1 km × 1 km, obtaining Level 3 data, to possibly improve their accuracy and capture the spatial variability to an even more appreciable spatial sampling. The chosen regridding technique is the ordinary kriging, a geostatistical technique that respect to other approaches, provide several advantages, such as the capability to provide both data estimations and associated errors and of giving satisfactory results also in the presence of few input data. Moreover, the main advantage of the kriging with respect to other deterministic techniques lies in the fact that not only information from geometric distribution of the known points is used, but also the knowledge of the spatial structure of the correlation between observations. The results allowed us to analyze the monthly variability of NO 2 concentration concerning two study areas, one in the North and the other in the South of Italy, and also to highlight different levels of NO 2 comparing weekdays and weekend days, and, finally, COVID-19 (Corona Virus Disease 2019) to non-COVID-19 times.
In addition, the correlation between satellite-based NO 2 data and in situ NO 2 concentrations measured by different air quality monitoring stations have also been analyzed. Recent work has already dealt with this issue finding a significant correlation (R 2 = 0.68) between Level 2 TROPOMI NO 2 products and ground-based observations and a weekend NO 2 level smaller than those observed during weekdays by about 30% in Helsinki [29]. OMI NO 2 observations have also been used to examine the NO 2 weekly load over urban cities [30,31].
Both TROPOMI NO 2 Level 2 and Level 3 data have been correlated against ground-based observations. In particular, ten air quality monitoring stations were considered for the area in the North of Italy, while thirty-seven air quality monitoring stations for the area in the South of Italy.
The paper is organized as follows: Section 2 provides a description of the data used for the implementation and the validation of the regridding procedure; it also describes the applied methodology, focusing on the construction of the empirical semi-variogram. Section 3 shows the results obtained from the application of the regridding method to Level 2 TROPOMI NO 2 data and their validation. In Section 4 we further discuss our results and methodology, and finally, conclusions are drawn in Section 5.

TROPOMI NO 2 Observations
TROPOMI is a single payload on board of the Copernicus S5P, one of the European Sentinel satellites dedicated to the monitoring of atmospheric composition, launched on the 13 October 2017. S5P is a near-polar sun-synchronous orbit satellite flying at an altitude of 817 km with an overpass time at ascending node of 13:30 and a repeat cycle of 17 days. TROPOMI operates in a push-broom configuration with a wide swath of about 2600 km, along-track of 7 km, which became 5.6 km from 6 August 2019, and a daily global coverage [32]. The instrument is a four spectrometers system that measures in the ultraviolet (UV), UV-visible (UV-VIS), near-infrared (NIR) and shortwave infrared (SWIR) spectral bands. From these spectral bands, several important air quality and climate-related atmospheric constituents can be retrieved, including NO 2 . In particular, the NO 2 columns retrieval algorithm developed by the Royal Netherlands Meteorological Institute (KNMI) exploits the UV-VIS TROPOMI measurements in the 405-465 nm wavelength range [33]. The TROPOMI NO 2 processing system is based on the DOMINO, an algorithm previously used for the OMI instrument [1,34,35]. The retrieval of the total NO 2 slant column density from the Level-1b radiance and irradiance spectra measured by TROPOMI is based on a Differential Optical Absorption Spectroscopy (DOAS) method [36]. The NO 2 slant column density is separated into a stratospheric and a tropospheric part based on Remote Sens. 2020, 12, 2212 4 of 23 information coming from a data assimilation system based on the TM5-MP chemical transport model. Tropospheric and stratospheric slant column densities are converted into vertical column densities by applying an appropriate air-mass factor (AMF) based on a look-up table of altitude-dependent AMFs and information on the vertical distribution of NO 2 from the TM5-MP chemistry transport model at a resolution of 1 × 1 degree and a time step of 30 min [37]. The altitude-dependent AMF depends on the satellite geometry, terrain height, cloud fraction and height and surface albedo.
In this work, the reprocessed Level 2 NO 2 tropospheric column data for the period May 2018-April 2020 have been used. The reprocessed version is the latest version of data that ensures uniformity of the data processing [1,38]. Two other versions of data are also available: the near-real-time stream, available within 3 h after data acquisition, and the offline stream, available within a few days after acquisition. A useful variable contained in NO 2 products, the quality assurance value (qa_value), a continuous variable ranging from 0 to 1, was also taken into consideration. In particular, a qa_value greater than 0.75, recommended removing cloud scenes covered by snow/ice, errors and problematic retrievals, was used as a pixel filter.
The two target areas are shown in Figures 1 and 2. The first area is the Po Valley, located in the North of Italy ( Figure 1, yellow box), bounded in longitude between 7 • E and 13 • E and in latitude between 44 • N and 46 • N. The second region is located in the South of Italy ( Figure 2, yellow box), bounded in longitude between 13 • E and 19 • E and in latitude between 39 • N and 42 • N.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 23 model at a resolution of 1 × 1 degree and a time step of 30 minutes [37]. The altitude-dependent AMF depends on the satellite geometry, terrain height, cloud fraction and height and surface albedo. In this work, the reprocessed Level 2 NO2 tropospheric column data for the period May 2018 -April 2020 have been used. The reprocessed version is the latest version of data that ensures uniformity of the data processing [1,38]. Two other versions of data are also available: the near-realtime stream, available within 3 hours after data acquisition, and the offline stream, available within a few days after acquisition. A useful variable contained in NO2 products, the quality assurance value ( _ ), a continuous variable ranging from 0 to 1, was also taken into consideration. In particular, a _ greater than 0.75, recommended removing cloud scenes covered by snow/ice, errors and problematic retrievals, was used as a pixel filter.
The two target areas are shown in Figures 1 and 2. The first area is the Po Valley, located in the North of Italy ( Figure 1, yellow box), bounded in longitude between 7° E and 13° E and in latitude between 44° N and 46° N. The second region is located in the South of Italy ( Figure 2, yellow box), bounded in longitude between 13° E and 19° E and in latitude between 39° N and 42° N.

Ground-Based Observations
For inter-comparison and validation, both Level 2 and Level 3 NO 2 TROPOMI data are compared vs. NO 2 observations obtained by air quality monitoring stations distributed on the two target areas. The stations are managed by Agenzia Regionale per la Protezione dell'Ambiente (ARPA) of Basilicata, Campania, Puglia, Lombardia and Emilia-Romagna Regions [39][40][41][42][43]. The network consists of 10 monitoring stations for the study area in the North and 37 in the South of Italy (see Figures 1  and 2), respectively. The surface concentration of NO 2 , measured in µg m 3 , is available on a time slot of 1 h. From these data, hourly values have been extracted, and daily/monthly averages have been computed, as well, for comparison with consistent satellite observations. The method used for the in situ measurement of NO 2 is based on the chemiluminescence reaction between NO and O 3 , which produces excited NO 2, which subsequently returns to its fundamental Remote Sens. 2020, 12, 2212 5 of 23 state by emitting electromagnetic radiation in the UV region. Through a catalytic converter, the NO 2 present in the air sample is reduced to NO. In this way, the total concentration of NOx is obtained, while the measurement of NO 2 is derived from the difference between total oxides and NO alone.
The instruments installed in the various considered air quality monitoring stations are constantly subjected to checks, aimed precisely at ensuring the proper functioning of the equipment over time and the reliability of the product data.

Ground-based Observations
For inter-comparison and validation, both Level 2 and Level 3 NO2 TROPOMI data are compared vs. NO2 observations obtained by air quality monitoring stations distributed on the two target areas. The stations are managed by Agenzia Regionale per la Protezione dell'Ambiente (ARPA) of Basilicata, Campania, Puglia, Lombardia and Emilia-Romagna Regions [39][40][41][42][43]. The network consists of 10 monitoring stations for the study area in the North and 37 in the South of Italy (see Figure 1 and Figure 2), respectively. The surface concentration of NO2, measured in , is available on a time slot of 1 hour. From these data, hourly values have been extracted, and daily/monthly averages have been computed, as well, for comparison with consistent satellite observations. The method used for the in situ measurement of NO2 is based on the chemiluminescence reaction between NO and O3, which produces excited NO2, which subsequently returns to its fundamental state by emitting electromagnetic radiation in the UV region. Through a catalytic converter, the NO2 present in the air sample is reduced to NO. In this way, the total concentration of NOx is obtained, while the measurement of NO2 is derived from the difference between total oxides and NO alone.
The instruments installed in the various considered air quality monitoring stations are constantly subjected to checks, aimed precisely at ensuring the proper functioning of the equipment over time and the reliability of the product data

Averaging and Remapping Operations
The native horizontal spatial resolution of NO2 TROPOMI data is 7 km × 3.5 km, and, of course, based on a single daily satellite overpass, we cannot improve the spatial resolution because we cannot increase information just by acting with mathematical tools. Information can be improved only on the basis of additional physical information. We mean other sources of satellite or correlative data, which is not the case for this study because our methodology relies on satellite data alone and suitable mathematics; no additional physics is used. It should be noted that if the implemented mathematical technique is inappropriate to solve a given problem, it could lead to loss of information [44,45].

Averaging and Remapping Operations
The native horizontal spatial resolution of NO 2 TROPOMI data is 7 km × 3.5 km, and, of course, based on a single daily satellite overpass, we cannot improve the spatial resolution because we cannot increase information just by acting with mathematical tools. Information can be improved only on the basis of additional physical information. We mean other sources of satellite or correlative data, which is not the case for this study because our methodology relies on satellite data alone and suitable mathematics; no additional physics is used. It should be noted that if the implemented mathematical technique is inappropriate to solve a given problem, it could lead to loss of information [44,45].
One way to improve information is to accumulate satellite data from consecutive overpasses for a given spatial region and for a time slot larger than the satellite repeat time. Our accumulation time slot width is one month. Since satellites on polar orbits revisit a given area at different angles and the Field of View (FOV) pattern does not point exactly at the same footprints, by accumulating data on a month makes it possible that a given region of extent less than the basic 7 km × 3.5 km grid mesh is visited more times with overlapping satellite FOVs. This is exemplified in Figure 3, which compares the clear sky footprints of just one TROPOMI overpass for an area of the Po Valley, and the accumulated overpasses over a whole month. It is seen that after one month, the clear sky footprints are much more densely distributed over the area than the single overpass.
The vast data gaps in Figure 3a are because of cloudiness, which makes it impossible to perform the retrieval for NO 2 . In contrast, the mosaic of overpasses accumulated over a monthly basis improves spatial data sampling, hence the colocation with in situ observations. These overlapping observations can be used to recover a map with improved spatial sampling. This can be done by acting with a suitable averaging filter (to smooth the effect of outliers) and remapping the data to a regular finer Remote Sens. 2020, 12, 2212 6 of 23 mesh grid than the original ones. These averaging and remapping operations are performed with our kriging methodology, which we will show in the next section.
After kriging operation, the data are labeled as Level 3, in contrast to the original (Level 2) sparse data. It is understood that, because of averaging, the kriging interpolation will happen at the expense of time resolution. In the presence of persistent processes, i.e., the continuous release of pollutants into the environment, this is not a serious drawback. In effect, air quality safety regulations do not look at a single, isolated event, which goes beyond thresholds for a few minutes, when at a persistent phenomenon that remains above limits for days or months. In this respect, our methodology is aiming at looking at persistent processes, which are those which can cause more risks for the health of human beings and the natural environment.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 23 Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing month makes it possible that a given region of extent less than the basic 7 km × 3.5 km grid mesh is visited more times with overlapping satellite FOVs. This is exemplified in Figure 3, which compares the clear sky footprints of just one TROPOMI overpass for an area of the Po Valley, and the accumulated overpasses over a whole month. It is seen that after one month, the clear sky footprints are much more densely distributed over the area than the single overpass. The vast data gaps in Figure 3a are because of cloudiness, which makes it impossible to perform the retrieval for NO2. In contrast, the mosaic of overpasses accumulated over a monthly basis improves spatial data sampling, hence the colocation with in situ observations. These overlapping observations can be used to recover a map with improved spatial sampling. This can be done by acting with a suitable averaging filter (to smooth the effect of outliers) and remapping the data to a regular finer mesh grid than the original ones. These averaging and remapping operations are performed with our kriging methodology, which we will show in the next section.  After kriging operation, the data are labeled as Level 3, in contrast to the original (Level 2) sparse data. It is understood that, because of averaging, the kriging interpolation will happen at the expense of time resolution. In the presence of persistent processes, i.e., the continuous release of pollutants into the environment, this is not a serious drawback. In effect, air quality safety regulations do not look at a single, isolated event, which goes beyond thresholds for a few minutes, when at a persistent phenomenon that remains above limits for days or months. In this respect, our methodology is aiming at looking at persistent processes, which are those which can cause more risks for the health of human beings and the natural environment.

Regridding Technique
This section briefly presents the ordinary kriging method used to pass from Level 2 NO2 data at a spatial resolution of 7 km × 3.5 km to Level 3 NO2 data at a spatial grid-resolution of 1 km × 1 km. The methodology is quite general and can be applied to any set of data.
The regridding technique, used to estimate a variable of interest at a given spatial location, is based on the spatial autocorrelation of data [46].
Let ( ) be a spatial process, whose values or realizations are known at a set of spatial locations, say, = , … , , where denotes the spatial coordinate, e.g., (longitude, latitude) in a case of 2-D process. Our objective is to estimate the process at a point ∉ = , … , for which no observation is available. The estimate, ( ) of ( ) is obtained according to [46] where λ is a vector of real weights and ( ) represents the process at points where it is known. In order to have an unbiased estimate, the additional condition of weights normalized to unity is assumed, i.e.,

Regridding Technique
This section briefly presents the ordinary kriging method used to pass from Level 2 NO 2 data at a spatial resolution of 7 km × 3.5 km to Level 3 NO 2 data at a spatial grid-resolution of 1 km × 1 km. The methodology is quite general and can be applied to any set of data.
The regridding technique, used to estimate a variable of interest at a given spatial location, is based on the spatial autocorrelation of data [46].
Let Z(s) be a spatial process, whose values or realizations are known at a set of spatial locations, say, s = {s 1 , . . . , s n }, where s j denotes the spatial coordinate, e.g., (longitude, latitude) in a case of 2-D process. Our objective is to estimate the process at a point s 0 s = {s 1 , . . . , s n } for which no observation is available. The estimate,Ẑ(s 0 ) of Z(s 0 ) is obtained according to [46] where λ is a vector of real weights and Z(s) represents the process at points where it is known. In order to have an unbiased estimate, the additional condition of weights normalized to unity is assumed, i.e., The weights must be properly determined to minimize (in the least Squares sense) the estimation of the variance and to ensure the estimator has no bias. In practice, the weights are determined based on the correlation properties of the process Z(s). These properties are described through a functional, called semivariogram, which yields the variance of the differences between pairs of the variable under study at two different sites s i and s j [46,47]. The semivariogram is defined according to: The semivariogram is a continuous function in the origin, but, in practice, it is often observed that γ(0) 0. This situation is known as nugget effect (γ(0) = τ 2 is called nugget), and it is linked to measurement errors or spatial resolution issues, which in turn affect the variogram value at very small distances.
Another important feature of the semivariogram function is the sill, defined as: where σ 2 is known as the partial sill (the sill minus the nugget), and h is the distance between two different observations. In case the sill takes on a finite value, it means that the stochastic process is weakly stationary; also, if this occurs for a finite value of h = h * , then h * is said the range of the variogram. The range quantifies the distance over which two different observations can be considered correlated. Based on the available data we can build up an empirical semivariogram, sayγ. In doing so, a fitting model must be considered to extrapolate the spatial behavior of the observed points to the area of interest. In the literature, there are several theoretical semivariogram models (e.g., linear, exponential, Gaussian, wave and circular), with known analytical properties and physical meaning of parameters.
In our methodology, we use a set of models, which are specified in Figure 3, and that which is best suited to our particular process Z(s) is chosen based on a statistical criterion, which selects the model that best fits the given data.
In this work, the fit criterion described in [48] has been used. This step of the analysis is again exemplified in Figure 4, where for the case shown, the best model is that spherical. It has to be stressed that the fit criterion procedure is applied to case-by-case, which means that for the analysis, we are going to show it has been applied to any single, individual, monthly maps. The weights must be properly determined to minimize (in the least Squares sense) the estimation of the variance and to ensure the estimator has no bias. In practice, the weights are determined based on the correlation properties of the process ( ). These properties are described through a functional, called semivariogram, which yields the variance of the differences between pairs of the variable under study at two different sites and [46,47]. The semivariogram is defined according to: The semivariogram is a continuous function in the origin, but, in practice, it is often observed that (0) 0. This situation is known as nugget effect ( (0) = is called nugget), and it is linked to measurement errors or spatial resolution issues, which in turn affect the variogram value at very small distances.
Another important feature of the semivariogram function is the sill, defined as: where is known as the partial sill (the sill minus the nugget), and ℎ is the distance between two different observations. In case the sill takes on a finite value, it means that the stochastic process is weakly stationary; also, if this occurs for a finite value of ℎ = ℎ * , then ℎ * is said the range of the variogram. The range quantifies the distance over which two different observations can be considered correlated.
Based on the available data we can build up an empirical semivariogram, say . In doing so, a fitting model must be considered to extrapolate the spatial behavior of the observed points to the area of interest. In the literature, there are several theoretical semivariogram models (e.g., linear, exponential, Gaussian, wave and circular), with known analytical properties and physical meaning of parameters.
In our methodology, we use a set of models, which are specified in Figure 3, and that which is best suited to our particular process ( ) is chosen based on a statistical criterion, which selects the model that best fits the given data.
In this work, the fit criterion described in [48] has been used. This step of the analysis is again exemplified in Figure 4, where for the case shown, the best model is that spherical. It has to be stressed that the fit criterion procedure is applied to case-by-case, which means that for the analysis, we are going to show it has been applied to any single, individual, monthly maps.   Once, the semivariogram model has been selected, the weights can be calculated and used for the spatial interpolation.
Finally, before finishing this section, it is also worth noting that the ordinary kriging methodology is a particular case of optimal interpolation and can be looked at as a special implementation of the most general formalism of Kalman filter (e.g., see [49][50][51][52]). This is important because, unlike ad-hoc interpolation methods, the kriging provides a robust and optimal tool for interpolation processes.

Comparing Satellite Data to in Situ Observations
In this section, we will compare in situ observations against Level 2 and Level 3 (krigged data). The final aim is to show that krigged data perform better as far as the correlation with ground-based observations is concerned.
The comparison between satellite data to in situ observations required the space-time colocation of data from the different sources. In particular, each in situ data have been compared with the spatially closest Level 2 and Level 3 TROPOMI value. About temporal colocation, each Level 2 TROPOMI observation is associated with the time of the satellite overpass, and therefore we used the temporally closest in situ measurement.
Toward this objective, the effects of the kriging procedure on Level 2 data is exemplified in Figure 5, which shows the scatterplot of time-space collocated values of TROPOMI NO2 and in situ measurements for the ten validation stations in North Italy (see Figure 1) before and after the kriging operation. In Figure 5a, the scatterplot refers to the observations accumulated over a month but not averaged or regridded, therefore they are instantaneous measurements. After the kriging operation, Figure 5b, the data are resampled monthly at the spatial sampling of 1 km × 1km. It is seen that the linear correlation coefficient increases from 0.50 to 0.85. Once, the semivariogram model has been selected, the weights λ can be calculated and used for the spatial interpolation.
Finally, before finishing this section, it is also worth noting that the ordinary kriging methodology is a particular case of optimal interpolation and can be looked at as a special implementation of the most general formalism of Kalman filter (e.g., see [49][50][51][52]). This is important because, unlike ad-hoc interpolation methods, the kriging provides a robust and optimal tool for interpolation processes.

Comparing Satellite Data to in Situ Observations
In this section, we will compare in situ observations against Level 2 and Level 3 (krigged data). The final aim is to show that krigged data perform better as far as the correlation with ground-based observations is concerned.
The comparison between satellite data to in situ observations required the space-time colocation of data from the different sources. In particular, each in situ data have been compared with the spatially closest Level 2 and Level 3 TROPOMI value. About temporal colocation, each Level 2 TROPOMI observation is associated with the time of the satellite overpass, and therefore we used the temporally closest in situ measurement.
Toward this objective, the effects of the kriging procedure on Level 2 data is exemplified in Figure 5, which shows the scatterplot of time-space collocated values of TROPOMI NO 2 and in situ measurements for the ten validation stations in North Italy (see Figure 1) before and after the kriging operation. In Figure 5a, the scatterplot refers to the observations accumulated over a month but not averaged or regridded, therefore they are instantaneous measurements. After the kriging operation, Figure 5b, the data are resampled monthly at the spatial sampling of 1 km × 1 km. It is seen that the linear correlation coefficient increases from 0.50 to 0.85.
From Figure 5, we see a good linear relationship extends over the range 10-80 µg m 3 . The results are showing that the TROPOMI can sense down to the surface and satellite data are capable of discriminating between low and high pollution load. Furthermore, Level 3 TROPOMI has a far better correlation with in situ observations. In passing, the good relationship we see for this case study shows the relatively high homogeneity of the 10 stations. In effect, they are within the bulk area of the Po Valley, which is affected by more intense and persistent air pollution concerning other regions in Italy. Figure 6 shows the results for the 37 validation stations in Southern Italy. In comparison to Figure 5, the correlation coefficient is smaller for both Level 2 and Level 3 data, but, as before, it increases passing form the Level 2 to 3. The smaller correlation is expected because, unlike the Po Valley target area, that in Southern Italy has a much more complex orography, and is much less homogeneous as far as Remote Sens. 2020, 12, 2212 9 of 23 the air quality conditions. The target area includes large urban areas, such as the city of Naples and countryside, sparsely populated towns, which results in very large variability of air pollution intensity. are showing that the TROPOMI can sense down to the surface and satellite data are capable of discriminating between low and high pollution load. Furthermore, Level 3 TROPOMI has a far better correlation with in situ observations. In passing, the good relationship we see for this case study shows the relatively high homogeneity of the 10 stations. In effect, they are within the bulk area of the Po Valley, which is affected by more intense and persistent air pollution concerning other regions in Italy. Figure 6 shows the results for the 37 validation stations in Southern Italy. In comparison to Figure 5, the correlation coefficient is smaller for both Level 2 and Level 3 data, but, as before, it increases passing form the Level 2 to 3. The smaller correlation is expected because, unlike the Po Valley target area, that in Southern Italy has a much more complex orography, and is much less homogeneous as far as the air quality conditions. The target area includes large urban areas, such as the city of Naples and countryside, sparsely populated towns, which results in very large variability of air pollution intensity.
At a different degree, the results in Figure 5 and Figure 6 show that the kriging technique is capable of improving the spatial sampling and yields a better comparison with in situ observations. However, the results apply to the whole target area. In contrast, it would be better to check the capability of TROPOMI to resolve air pollution at the scale of the single station. Towards this objective, we now come to compare satellite vs. in situ observations at every single station. For the Po Valley, the results are summarized in Table 1, where the correlation between TROPOMI Level 3 and in situ observations is shown. In general, we have found that the correlation improves when we move from Level 2 to Level 3, as exemplified for one of the ten stations listed in Tab. 1, in Figure 7, which shows the results for the Bologna-De Amicis station. The station is located within an urban area densely populated and with intense traffic all year-round. As said, we observe From Figure 5, we see a good linear relationship extends over the range 10-80 . The results are showing that the TROPOMI can sense down to the surface and satellite data are capable of discriminating between low and high pollution load. Furthermore, Level 3 TROPOMI has a far better correlation with in situ observations. In passing, the good relationship we see for this case study shows the relatively high homogeneity of the 10 stations. In effect, they are within the bulk area of the Po Valley, which is affected by more intense and persistent air pollution concerning other regions in Italy. Figure 6 shows the results for the 37 validation stations in Southern Italy. In comparison to Figure 5, the correlation coefficient is smaller for both Level 2 and Level 3 data, but, as before, it increases passing form the Level 2 to 3. The smaller correlation is expected because, unlike the Po Valley target area, that in Southern Italy has a much more complex orography, and is much less homogeneous as far as the air quality conditions. The target area includes large urban areas, such as the city of Naples and countryside, sparsely populated towns, which results in very large variability of air pollution intensity.
At a different degree, the results in Figure 5 and Figure 6 show that the kriging technique is capable of improving the spatial sampling and yields a better comparison with in situ observations. However, the results apply to the whole target area. In contrast, it would be better to check the capability of TROPOMI to resolve air pollution at the scale of the single station. Towards this objective, we now come to compare satellite vs. in situ observations at every single station. For the Po Valley, the results are summarized in Table 1, where the correlation between TROPOMI Level 3 and in situ observations is shown. In general, we have found that the correlation improves when we move from Level 2 to Level 3, as exemplified for one of the ten stations listed in Tab. 1, in Figure 7, which shows the results for the Bologna-De Amicis station. The station is located within an urban area densely populated and with intense traffic all year-round. As said, we observe At a different degree, the results in Figures 5 and 6 show that the kriging technique is capable of improving the spatial sampling and yields a better comparison with in situ observations. However, the results apply to the whole target area. In contrast, it would be better to check the capability of TROPOMI to resolve air pollution at the scale of the single station.
Towards this objective, we now come to compare satellite vs. in situ observations at every single station. For the Po Valley, the results are summarized in Table 1, where the correlation between TROPOMI Level 3 and in situ observations is shown. In general, we have found that the correlation improves when we move from Level 2 to Level 3, as exemplified for one of the ten stations listed in Tab. 1, in Figure 7, which shows the results for the Bologna-De Amicis station. The station is located within an urban area densely populated and with intense traffic all year-round. As said, we observe a poor correlation with Level 2 data. Conversely, we observe a very high correlation, which reaches the value of R 2 = 0.93 when going to Level 3 data. The same good correlation is also seen when we consider monthly time series as done in Figure 8. It is seen that the krigged data are capable to follow the seasonal cycle, showing an increase of NO 2 pollution in wintertime months because of the combined effect of domestic heating and increased traffic.
Overall, the good correlation we have found for the Po Valley is an effect of the better spatial sampling of L3 data, but we think it is also a consequence of the particular, local, weather conditions, which govern air circulation over the valley. The Po river flows in a closed valley surrounded by high mountains, which makes pollution remains trapped, especially in wintertime, when atmospheric conditions prevent dispersal. a poor correlation with Level 2 data. Conversely, we observe a very high correlation, which reaches the value of = 0.93 when going to Level 3 data. The same good correlation is also seen when we consider monthly time series as done in Figure 8. It is seen that the krigged data are capable to follow the seasonal cycle, showing an increase of NO2 pollution in wintertime months because of the combined effect of domestic heating and increased traffic.    Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing the value of = 0.93 when going to Level 3 data. The same good correlation is also seen when we consider monthly time series as done in Figure 8. It is seen that the krigged data are capable to follow the seasonal cycle, showing an increase of NO2 pollution in wintertime months because of the combined effect of domestic heating and increased traffic.    TROPOMI data, as well as satellite observations, estimate tropospheric NO 2 columns rather than ground-level NO 2 concentrations, which, in turn, is normally measured from in situ stations. In addition, satellite data cannot resolve the diurnal cycle because of the polar orbits. The good correlation we observe means that, because of the huge pollution stain over the Po Valley, NO 2 tends to be compressed and trapped within the PBL. We have a persistent phenomenon, which can be measured by satellite on consecutive days. The air pollution persistence is an effect peculiar to the Po Valley and its local weather.
Other Italian regions can have local weather conditions, which can favor the dispersal of pollution. The complex orography of most of the Italian territory can impact on the PBL dynamics. As a result, we can observe vertical and horizontal inhomogeneities of NO 2 mixing ratio profiles.
In effect, Table 2 (South Italy) shows a picture much more complex. In general, higher correlations are found for densely populated cities on the East coast (Adriatic Sea) and inner part (mountains). However, there is not a clear pattern. As a result of the above discussion, the cases of fair correlation (≈0.50) could be related to the particular orography of the given area, which could influence the local weather and climate. As an example, for the city of Naples, which we know is an area of intense air pollution, likely, the local, very complex, orography and the related weather characteristics contribute to increasing the day-to-day variability of air quality and, much important, increase the vertical and horizontal inhomogeneity of NO 2 . In general, using additional physical information, e.g., see [53], the correlation could also be improved for stations showing smaller consistency between the tropospheric column and near-surface measurements. However, this is a task not addressed in this paper, which rather focuses on the direct use of TROPOMI observations per se. Nevertheless, it is fair to recall and caution the reader that the TROPOMI observes NO 2 columns rather than ground-level NO 2 concentrations, which are closely related to potentially harmful effects on human health.
Before closing this section, we stress that also for the second target area, we have that, in general, the correlation improves when we go from Level 2 to Level 3. This is exemplified for two stations in Figures 9-12. The first station refers to the urban area of the Avellino city in the mountains. Figure 9 shows that the correlation clearly improves after the kriging applications to the Level 2 data. In general, using additional physical information, e.g., see [53], the correlation could also be improved for stations showing smaller consistency between the tropospheric column and nearsurface measurements. However, this is a task not addressed in this paper, which rather focuses on the direct use of TROPOMI observations per se. Nevertheless, it is fair to recall and caution the reader that the TROPOMI observes NO2 columns rather than ground-level NO2 concentrations, which are closely related to potentially harmful effects on human health.
Before closing this section, we stress that also for the second target area, we have that, in general, the correlation improves when we go from Level 2 to Level 3. This is exemplified for two stations in Figs. 9 to 12. The first station refers to the urban area of the Avellino city in the mountains. Figure 9 shows that the correlation clearly improves after the kriging applications to the Level 2 data. The good correlation (Level 3) is also reflected in Figure 10, where we show the comparison between monthly time series. TROPOMI data are capable of following and reproduce the annual cycle, with the winter peak due to domestic heating. Figure 11 shows the scatter plot for the second station, which is the Viggiano-Costa Molina station. This station is located in a countryside area some km apart from an oil refinery plant. Still, it is far from urban areas to monitor the supposedly NOx emissions from the industrial plant.
(a) Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing surface measurements. However, this is a task not addressed in this paper, which rather focuses on the direct use of TROPOMI observations per se. Nevertheless, it is fair to recall and caution the reader that the TROPOMI observes NO2 columns rather than ground-level NO2 concentrations, which are closely related to potentially harmful effects on human health. Before closing this section, we stress that also for the second target area, we have that, in general, the correlation improves when we go from Level 2 to Level 3. This is exemplified for two stations in Figs. 9 to 12. The first station refers to the urban area of the Avellino city in the mountains. Figure 9 shows that the correlation clearly improves after the kriging applications to the Level 2 data. The good correlation (Level 3) is also reflected in Figure 10, where we show the comparison between monthly time series. TROPOMI data are capable of following and reproduce the annual cycle, with the winter peak due to domestic heating. Figure 11 shows the scatter plot for the second station, which is the Viggiano-Costa Molina station. This station is located in a countryside area some km apart from an oil refinery plant. Still, it is far from urban areas to monitor the supposedly NOx emissions from the industrial plant. Unlike the former Avellino station, now the correlation coefficient goes to zero in the case of Level 2 data. For this station, we observe no sign of important pollution, as it can be argued by comparing Figure 11 to Figure 9. Nevertheless, when we pass to Level 3 data, it is seen that TROPOMI is still sensitive to these very low values of NO2 concentrations. According to the present Italian regulation the NO2 annual average should not exceed 40 . The monthly averages for Viggiano are smaller by a factor of ten.  Unlike the former Avellino station, now the correlation coefficient goes to zero in the case of Level 2 data. For this station, we observe no sign of important pollution, as it can be argued by comparing Figure 11 to Figure 9. Nevertheless, when we pass to Level 3 data, it is seen that TROPOMI is still sensitive to these very low values of NO2 concentrations. According to the present Italian regulation the NO2 annual average should not exceed 40 . The monthly averages for Viggiano are smaller by a factor of ten. Although relatively small, the correlation found for the Level 3 data is still large enough to correctly reproduce the salient features of the seasonal cycle, as shown in Figure 12. Although relatively small, the correlation found for the Level 3 data is still large enough to correctly reproduce the salient features of the seasonal cycle, as shown in Figure 12.   The good correlation (Level 3) is also reflected in Figure 10, where we show the comparison between monthly time series. TROPOMI data are capable of following and reproduce the annual cycle, with the winter peak due to domestic heating. Figure 11 shows the scatter plot for the second station, which is the Viggiano-Costa Molina station. This station is located in a countryside area some km apart from an oil refinery plant. Still, it is far from urban areas to monitor the supposedly NOx emissions from the industrial plant.
Unlike the former Avellino station, now the correlation coefficient goes to zero in the case of Level 2 data. For this station, we observe no sign of important pollution, as it can be argued by comparing Figure 11 to Figure 9. Nevertheless, when we pass to Level 3 data, it is seen that TROPOMI is still sensitive to these very low values of NO 2 concentrations. According to the present Italian regulation the NO 2 annual average should not exceed 40 µg m 3 . The monthly averages for Viggiano are smaller by a factor of ten.
Although relatively small, the correlation found for the Level 3 data is still large enough to correctly reproduce the salient features of the seasonal cycle, as shown in Figure 12.

Monthly Spatial Map of Level 3 NO 2 Tropospheric Column
With the exceptions and cautions outlined above, the previous section gives us the confidence of the good reliability of TROPOMI, especially for the Po Valley, to follow the lower tropospheric column amount of NO 2 and related evolution and variations. We now pass to illustrate and describe spatial patterns of NO 2 tropospheric column data, as revealed from the TROPOMI instrument.
To begin with, Figure 13 illustrates the monthly NO 2 distribution over the Po Valley area, for the period between May 2018 to April 2019. These maps show persistent NO 2 hotspots over big cities, such as Milano (star symbol), and widespread NO 2 distribution in the colder months, particularly from November 2018 to February 2019, with peaks three times larger than those reached in the months of May-August 2018 and April 2019. This situation underlines the fact that NO 2 concentrations are strictly related to a more massive vehicular traffic (because of school season) and domestic heating, and in addition to less sunlight and colder temperatures, typical of the winter months.
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing such as Milano (star symbol), and widespread NO2 distribution in the colder months, particularly from November 2018 to February 2019, with peaks three times larger than those reached in the months of May-August 2018 and April 2019. This situation underlines the fact that NO2 concentrations are strictly related to a more massive vehicular traffic (because of school season) and domestic heating, and in addition to less sunlight and colder temperatures, typical of the winter months.  Figure 14 illustrates the monthly maps of Level 3 NO2 tropospheric column data related to the study area in the South of Italy. These maps show a persistent NO2 hotspot in correspondence to the city of Naples (star symbol). Moreover, also for this case, we see a seasonal cycle, with a higher level  Figure 14 illustrates the monthly maps of Level 3 NO 2 tropospheric column data related to the study area in the South of Italy. These maps show a persistent NO 2 hotspot in correspondence to the city of Naples (star symbol). Moreover, also for this case, we see a seasonal cycle, with a higher level of NO 2 during the months from October 2018 to February 2019. The comparison between the two target areas (North and South) shows that the densely populated and heavily industrialized Po Valley is affected by larger spatial patterns of NO 2 , with concentration peaks three times larger than those reached in the Southern area.
As said before, the main reasons for the huge pollution stain over the Po Valley are closely connected to intense urban and industrial activities and to orographic aspects, which can determine unfavorable local weather conditions, which make pollution remain trapped, especially in wintertime, and prevent dispersal.
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing is affected by larger spatial patterns of NO2, with concentration peaks three times larger than those reached in the Southern area. As said before, the main reasons for the huge pollution stain over the Po Valley are closely connected to intense urban and industrial activities and to orographic aspects, which can determine unfavorable local weather conditions, which make pollution remain trapped, especially in wintertime, and prevent dispersal.

Weekdays and Weekend Observations
To show the capability of the Level 3 data to detect air quality behavior because of anthropic cyclic activities, in this section, we examine NO2 load on weekdays and weekend times.

Weekdays and Weekend Observations
To show the capability of the Level 3 data to detect air quality behavior because of anthropic cyclic activities, in this section, we examine NO 2 load on weekdays and weekend times.
The analysis of weekly cycles captures lower NO 2 levels during the weekend days compared to the weekdays, as a result of reduced emissions from traffic and industrial activities, which are more intense in urban areas. Differently to previous work [29], we only used two days for weekdays and not all five days from Monday to Friday to have uniformity with the two days of the weekend.
More in detail, 16 weekdays (Wednesday and Thursday) and 16 weekend days for June and July were considered for the two target areas. The analysis has been performed in summer so that the results are not affected by domestic heating. The results are shown in Figure 15 (Southern Italy) and Figure 16 (Po Valley), and for both cases, it is possible to note the higher NO 2 concentrations over urban districts, such as Naples and Milan.
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing More in detail, 16 weekdays (Wednesday and Thursday) and 16 weekend days for June and July were considered for the two target areas. The analysis has been performed in summer so that the results are not affected by domestic heating. The results are shown in Figure 15 (Southern Italy) and Figure 16 (Po Valley), and for both cases, it is possible to note the higher NO2 concentrations over urban districts, such as Naples and Milan.
In passing, we draw attention to a striking feature, which can be seen in Figure 16a. If we draw the diagonal from top left to bottom right, it is seen that the NO2 pollution tends to be higher along that direction, from the bottom right to the middle of the map. In effect, that is the main direction of the heavily trafficked Italian A1 highway.

NO2 Pollution During COVID-19 Epidemic Crisis
While writing this paper, Italy has been hit by the dramatic COVID-19 epidemic. As a result that the science community has drawn attention to the possible interrelationship between COVID-19 diffusion and air quality, we have assessed the evolution of NO2 pollution over the Po Valley in the period Dec. 2019 to Apr. 2020. In this period, COVID-19 has spread and supposedly infected millions of people, leading Italian authorities to declare the lockdown of the Po Valley area, which has been later extended to the whole Nation. We stress that we just put attention on the NO2 evolution per se; in other words, we just show how this evolution has been observed and photographed from space.
For a proper comparison, the monthly maps of NO2 concentration over the Po Valley for the period December to March for both the COVID-19 (that is 2019-2020) and the no COVID (2018-2019) periods are shown in Figs. 17 and 18, respectively. The first striking feature which arises from the Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing were considered for the two target areas. The analysis has been performed in summer so that the results are not affected by domestic heating. The results are shown in Figure 15 (Southern Italy) and Figure 16 (Po Valley), and for both cases, it is possible to note the higher NO2 concentrations over urban districts, such as Naples and Milan. In passing, we draw attention to a striking feature, which can be seen in Figure 16a. If we draw the diagonal from top left to bottom right, it is seen that the NO2 pollution tends to be higher along that direction, from the bottom right to the middle of the map. In effect, that is the main direction of the heavily trafficked Italian A1 highway.

NO2 Pollution During COVID-19 Epidemic Crisis
While writing this paper, Italy has been hit by the dramatic COVID-19 epidemic. As a result that the science community has drawn attention to the possible interrelationship between COVID-19 diffusion and air quality, we have assessed the evolution of NO2 pollution over the Po Valley in the period Dec. 2019 to Apr. 2020. In this period, COVID-19 has spread and supposedly infected millions of people, leading Italian authorities to declare the lockdown of the Po Valley area, which has been later extended to the whole Nation. We stress that we just put attention on the NO2 evolution per se; in other words, we just show how this evolution has been observed and photographed from space.
For a proper comparison, the monthly maps of NO2 concentration over the Po Valley for the period December to March for both the COVID-19 (that is 2019-2020) and the no COVID (2018-2019) periods are shown in Figs. 17 and 18, respectively. The first striking feature which arises from the In passing, we draw attention to a striking feature, which can be seen in Figure 16a. If we draw the diagonal from top left to bottom right, it is seen that the NO 2 pollution tends to be higher along that direction, from the bottom right to the middle of the map. In effect, that is the main direction of the heavily trafficked Italian A1 highway.

NO 2 Pollution During COVID-19 Epidemic Crisis
While writing this paper, Italy has been hit by the dramatic COVID-19 epidemic. As a result that the science community has drawn attention to the possible interrelationship between COVID-19 diffusion and air quality, we have assessed the evolution of NO 2 pollution over the Po Valley in the period Dec. 2019 to Apr. 2020. In this period, COVID-19 has spread and supposedly infected millions of people, leading Italian authorities to declare the lockdown of the Po Valley area, which has been later extended to the whole Nation. We stress that we just put attention on the NO 2 evolution per se; in other words, we just show how this evolution has been observed and photographed from space.
For a proper comparison, the monthly maps of NO 2 concentration over the Po Valley for the period December to March for both the COVID-19 (that is 2019-2020) and the no COVID (2018-2019) periods are shown in Figures 17 and 18, respectively. The first striking feature which arises from the comparison of the corresponding maps for no-COVID-19 ( Figure 17) and COVID-19 ( Figure 18) is the sharp decrease of air pollution in 2020 compared to 2019. This has been most likely the effect of the warmer 2019-20 winter, which has limited the use of heating power for both residential and commercial sectors. It is to be stressed that according to ECMWF, winter 2020 has been the far warmest winter season ever recorded over Europe (e.g., see https://climate.copernicus.eu/boreal-winter-season-1920was-far-warmest-winter-season-ever-recorded-europe-0).
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing commercial sectors. It is to be stressed that according to ECMWF, winter 2020 has been the far warmest winter season ever recorded over Europe (e.g., see https://climate.copernicus.eu/borealwinter-season-1920-was-far-warmest-winter-season-ever-recorded-europe-0). By comparing Figure 17 to Figure 18, we see that only January 2019 and 2020 show comparable NO2 pollution, whereas February 2020 shows a dramatic decrease of NO2 pollution over the whole Po Valley. This month has been supposedly that of maximum virus diffusion among the population.
In March 2020, when the lockdown was active, we see some 50% decrease of NO2 pollution concerning the same month of 2019. The same situation also occurs for April, which is shown in Figure 19. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing warmest winter season ever recorded over Europe (e.g., see https://climate.copernicus.eu/borealwinter-season-1920-was-far-warmest-winter-season-ever-recorded-europe-0). By comparing Figure 17 to Figure 18, we see that only January 2019 and 2020 show comparable NO2 pollution, whereas February 2020 shows a dramatic decrease of NO2 pollution over the whole Po Valley. This month has been supposedly that of maximum virus diffusion among the population.
In March 2020, when the lockdown was active, we see some 50% decrease of NO2 pollution concerning the same month of 2019. The same situation also occurs for April, which is shown in Figure 19.  By comparing Figure 17 to Figure 18, we see that only January 2019 and 2020 show comparable NO 2 pollution, whereas February 2020 shows a dramatic decrease of NO 2 pollution over the whole Po Valley. This month has been supposedly that of maximum virus diffusion among the population.
In March 2020, when the lockdown was active, we see some 50% decrease of NO 2 pollution concerning the same month of 2019. The same situation also occurs for April, which is shown in Figure 19. To sum up, we can say that during the lockdown, NO2 has decreased over the Po Valley. However, this decrease is over-imposed to a general decreasing trend because of the favorable 2020 weather conditions. The general reduction in NO2 pollution during the 2020 winter is also clearly seen at ground stations and therefore is not an artifact of satellite data. This is exemplified in Tab. 3, which shows the percentage variation among 2019 and 2020 winter months as measured at ground stations in Bergamo, Brescia and Milan, respectively. These three cities are those who share the bulk of NO2 concentrations over the Po Valley. Furthermore, the negative variations are confirmed with TROPOMI data, as can be directly checked by the same Table 3.

Discussion
Kriging and related tools are now of widespread use within the context of Geostatistics and the analysis of Geospatial data (e.g., [54][55][56]). Application to the particular area of air quality with a focus on NO2 has been recently reported in the literature [57,58], which emphasizes the use of in situ observations. Our application to TROPOMI is quite new, and it has been designed to improve the spatial sampling of the data, which is inevitably affected by gaps because of cloudy conditions. In To sum up, we can say that during the lockdown, NO 2 has decreased over the Po Valley. However, this decrease is over-imposed to a general decreasing trend because of the favorable 2020 weather conditions. The general reduction in NO 2 pollution during the 2020 winter is also clearly seen at ground stations and therefore is not an artifact of satellite data. This is exemplified in Table 3, which shows the percentage variation among 2019 and 2020 winter months as measured at ground stations in Bergamo, Brescia and Milan, respectively. These three cities are those who share the bulk of NO 2 concentrations over the Po Valley. Furthermore, the negative variations are confirmed with TROPOMI data, as can be directly checked by the same Table 3.

Discussion
Kriging and related tools are now of widespread use within the context of Geostatistics and the analysis of Geospatial data (e.g., [54][55][56]). Application to the particular area of air quality with a focus on NO 2 has been recently reported in the literature [57,58], which emphasizes the use of in situ observations. Our application to TROPOMI is quite new, and it has been designed to improve the spatial sampling of the data, which is inevitably affected by gaps because of cloudy conditions. In fact, we use a mosaic of TROPOMI overpasses, which is used to estimate the semivariogram that governs the data structure at the different spatial scales. Once we have estimated the semivariogram, we perform a regridding of the data on a finer, regular, grid mesh with a step ∆r =1 km. The regridding fills the gaps in the original footprint geometry and improves the colocation with in situ observations. We stress that the step ∆r =1 km should not be confused with the spatial resolution of the original data, which remains fixed to that of the TROPOMI footprint geometry. In fact, in our kriging procedure, there is no explicit de-convolution of the original observations or information from additional sources of data at better spatial resolution. Using regridding we oversample the original TROPOMI footprints to match the position of ground stations better and perform time-averaging of the data. It could be argued that this remapping could be performed in many different ways involving mathematical interpolation tools. The specificity of kriging methods is that they do not only consider the distance between observations, but they also capture the spatial structure in the data by computing the variogram and using it to set the weights to perform the interpolation at un-sampled points in the grid. To show the kriging advantages, we have also produced NO 2 maps for the Po Valley using the technique of the nearest neighbor combined with the ordinary arithmetic mean. In doing so, Level 2 data were accumulated for each month on a regular, fixed, grid which has the size of the TROPOMI original spatial resolution of 7 km × 3.5 km. For each satellite overpass, the nodes of this regular grid were filled with the spatially closest Level 2 value (within a maximum distance that refers to the original Level 2 TROPOMI spatial resolution) and, subsequently, the regularly spatially distributed values were averaged over the month. In this way, we can obtain at each node a value which is representative of the monthly averaged Level 2 NO 2 TROPOMI tropospheric column. The results are shown in Figure 20 for the station of Bologna, so that they can be directly compared to those of Figure 7. It is seen that we improve the correlation with respect to Level 2 data, but not concerning kriging for which we observe a correlation of 0.93 against that of 0.80 shown in Figure 20.
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing governs the data structure at the different spatial scales. Once we have estimated the semivariogram, we perform a regridding of the data on a finer, regular, grid mesh with a step ∆ =1 km. The regridding fills the gaps in the original footprint geometry and improves the colocation with in situ observations. We stress that the step ∆ =1 km should not be confused with the spatial resolution of the original data, which remains fixed to that of the TROPOMI footprint geometry. In fact, in our kriging procedure, there is no explicit de-convolution of the original observations or information from additional sources of data at better spatial resolution. Using regridding we oversample the original TROPOMI footprints to match the position of ground stations better and perform timeaveraging of the data. It could be argued that this remapping could be performed in many different ways involving mathematical interpolation tools. The specificity of kriging methods is that they do not only consider the distance between observations, but they also capture the spatial structure in the data by computing the variogram and using it to set the weights to perform the interpolation at unsampled points in the grid. To show the kriging advantages, we have also produced NO2 maps for the Po Valley using the technique of the nearest neighbor combined with the ordinary arithmetic mean. In doing so, Level 2 data were accumulated for each month on a regular, fixed, grid which has the size of the TROPOMI original spatial resolution of 7 km × 3.5 km. For each satellite overpass, the nodes of this regular grid were filled with the spatially closest Level 2 value (within a maximum distance that refers to the original Level 2 TROPOMI spatial resolution) and, subsequently, the regularly spatially distributed values were averaged over the month. In this way, we can obtain at each node a value which is representative of the monthly averaged Level 2 NO2 TROPOMI tropospheric column. The results are shown in Figure 20 for the station of Bologna, so that they can be directly compared to those of Figure 7. It is seen that we improve the correlation with respect to Level 2 data, but not concerning kriging for which we observe a correlation of 0.93 against that of 0.80 shown in Figure 20. Another important aspect to discuss is that regarding the choice of the grid step ∆ =1 km. Although there is no general prescription on how to choose the grid step, the case study at hand is normally suggesting a suitable pragmatic and effective choice. It is normally suggested (see e.g., [56], section 5.4.2) that ∆ "ℎ * (ℎ * is the range that for our analysis is of the order of 100 km, e.g., see Figure  3) and that the number of couples in any grid-box should range between 6-24 (see e.g., [55], Section 5.4.2). In addition, for our application, it is obvious to require ∆ to be smaller than the average radius of the TROPOMI footprint. Otherwise, as shown in Figure 20, we could not take advantage of the oversampling, which is allowed with our mosaic of monthly scenes. If we take ∆ larger than the TROPOMI footprint, there would be no point in oversampling since we could fill gaps with available clear sky data, if any, for the given scene. Finally, according to [55], at least 6-12 points should be used for each predicted value. In our case, with ∆ =1 km we have that each grid box has ≈ 10 neighbor pairs to be used for interpolation. In this way, for any grid point, considering the four quadrants, there are some 40 neighbors to be considered for interpolation. This density of satellite observations also ensures that we have enough points covering the time window of one month. This Another important aspect to discuss is that regarding the choice of the grid step ∆r =1 km. Although there is no general prescription on how to choose the grid step, the case study at hand is normally suggesting a suitable pragmatic and effective choice. It is normally suggested (see e.g., [56], section 5.4.2) that ∆r"h * (h * is the range that for our analysis is of the order of 100 km, e.g., see Figure 3) and that the number of couples in any grid-box should range between 6-24 (see e.g., [55], Section 5.4.2). In addition, for our application, it is obvious to require ∆r to be smaller than the average radius of the TROPOMI footprint. Otherwise, as shown in Figure 20, we could not take advantage of the oversampling, which is allowed with our mosaic of monthly scenes. If we take ∆r larger than the TROPOMI footprint, there would be no point in oversampling since we could fill gaps with available clear sky data, if any, for the given scene. Finally, according to [55], at least 6-12 points should be used for each predicted value. In our case, with ∆r = 1 km we have that each grid box has ≈ 10 neighbor pairs to be used for interpolation. In this way, for any grid point, considering the four quadrants, there are some 40 neighbors to be considered for interpolation. This density of satellite observations also ensures that we have enough points covering the time window of one month. This allows us to take into account the day-to-day variability, which is needed to have a good estimate of the monthly mean.

Conclusions
In this paper, a geostatistical regridding technique, i.e., the ordinary kriging, was applied to improve the spatial resolution of satellite-based NO 2 tropospheric column estimations from the original Remote Sens. 2020, 12, 2212 20 of 23 coarser resolution (7 km × 3.5 km) to a final finer grid-resolution (1 km × 1 km). The finer resolution allowed us to explore spatial data variability, obtaining more information on the main NO 2 sources.
The analysis of the consistency between satellite and ground-based NO 2 values showed good correlations that demonstrated the capability of the TROPOMI instrument to well reproduce the variability of NO 2 concentrations in the lower troposphere, especially for the Po Valley. This consistency has been tested both for the original Level 2 and for the Level 3 NO 2 data, obtaining an improvement of correlation coefficients when we move from Level 2 to 3.
In effect, once properly filtered and remapped with our kriging procedure, TROPOMI NO 2 concentration showed good consistency with in situ observations. This consistency was assessed through the calculation of the linear correlation coefficients, which reached values larger than 0.90 for severely polluted areas, such as those in the Po Valley. This is an important result in that it can allow us to use TROPOMI data in air quality protocols to check and monitor, e.g., the effect of ad hoc measures to mitigate and lower the intensity of air pollution. This could have a high impact on the European Copernicus Program and boost the use of the European satellite platform.
Having said that, we have also found, for the target area in southern Italy, much more mixed results, which are likely because of an unstable PBL, which could well condition the local weather in a way to increase the vertical and horizontal inhomogeneity of NO 2 . For Southern Italy, there remain cases for which from the satellite data, we can observe pollution spots, whose magnitude and intensity cannot be easily extrapolated from the troposphere to the air close to the surface.
We have also used TROPOMI data to analyze NO 2 pollution during COVID-19 and no-COVID-19 times. We have shown that TROPOMI observations have been sensitive to the lockdown effect on air quality, showing a general decrease of NO 2 concentrations in March and April 2020. However, at the same time TROPOMI, in agreement with in situ observations, has shown a general decrease of lower tropospheric NO 2 to 2019, mostly because of a warmer 2020 winter that has limited the use of power for residential and commercial heating.