A novel gridding algorithm to create regional trace gas maps from satellite observations

Introduction Conclusions References


Introduction
Satellite retrieved global and regional maps of trace gases are important in many areas of atmospheric research.They provide details about the spatial and temporal distribution of atmospheric components, such as nitrogen dioxide (NO 2 ) or sulphur dioxide (SO 2 ), which cannot be obtained by ground-based networks or field campaigns alone.Satellite maps have been used to study emissions and transport of trace gases from natural and anthropogenic sources on both a regional and global scale (e.g.Carn et al., 2007).Top-down emission inventories are used to constrain emissions in chemistry transport models (e.g.Martin et al., 2003;Choi et al., 2008).Furthermore, trends and cycles of trace gas concentrations can be identified over regions on different time scales (e.g.Beirle et al., 2003;Richter et al., 2005).Validation studies, comparing space-and ground-based measurements or model simulation, require or can benefit from high-resolution maps to identify the regional distribution of trace gases (e.g.Krotkov et al., 2008;Wenig et al., 2008;Chan et al., 2012).In recent years, the increase of the spatial resolution of satellite instruments made trace gas products more suitable for applications on a regional scale.Examples are the Empa OMI NO 2 product (Boersma et al., 2007;Zhou et al., 2009Zhou et al., , 2010) ) and the Berkeley High-Resolution (BEHR) product (Russell et al., 2011).These products are retrieved from the Ozone Monitoring Instrument (OMI), which has ground pixel sizes from 13 × 24 km 2 at nadir to 30 × 160 km 2 at the edge of the swath (Levelt et al., 2006).
Vertical column densities of trace gases are typically expressed in the instrument's frame of reference using across-track and along-track position (Level 2 product).To visualize and compare the data with other data sets the Level 2 product is projected on a longitude-latitude grid (Level 3 product) using a suitable gridding algorithm.On a regional scale grids may have resolutions between 1 × 1 km well within the limitation set by the instrument's resolution.It can also be assumed that the distribution is continuous within an orbit of measurements.Current gridding algorithms typically assume that measurement values are constant within the pixel boundaries.These algorithms produce non-continuous maps and can underestimate maximum values.Within one orbit, consecutive measurements can have overlapping ground pixels due to the instrument slit function.Therefore, Wenig et al. (2008) average pixels based on measurement error and pixel size.While their approach conserves the mean value of the measurement, averaging overlapping pixels can further reduce maximum values.A different approach is suggested by Brunner et al. (2007) who tested image reconstruction algorithms used in computed tomography on idealized test data.However, their method targets the reconstruction of monthly or annual maps using a large number of overlapping pixels.
In this paper, we present a novel gridding algorithm which was developed to be more suitable to create high-resolution trace gas maps.The algorithm is designed for OMI, but can be employed easily to similar satellite instruments, for example, the TROPOspheric Monitoring Instrument (TROPOMI), planned for launch in 2014 (Veefkind et al., 2012).TROPOMI is expected to have smaller measurement errors than OMI, due to higher spatial resolution (7 × 7 km 2 at nadir) and better cloud correction.
The new gridding algorithm approximates the trace gas distribution by a continuous differentiable, parabolic spline surface, defined on the lattice of tiled pixels in a single orbit.Our parabolic spline method (PSM) is based on an algorithm developed by Kobza and Mlcak (1994), who construct a spline surface from known mean values on a rectangular lattice.The use of polynomial splines for interpolation is well established (e.g.Lancaster and Salkauskas, 1986;de Boor, 2001).Polynomial approximations which conserve mean values within lattice cells are called histopolations, because they are often used for smoothing histograms (Späth, 1995).OMI column densities can be interpreted as weighted averages over the area of a ground pixel which spans several lattice cells.The spatial sensitivity of the sensor over ground is described by the instrument function.PSM derives the mean values within each cell of the lattice considering this Introduction

Conclusions References
Tables Figures

Back Close
Full instrument function.The inverse problem of computing the mean values from the instrument function is ill-posed in the sense that the retrieved solutions are very sensitive to measurement noise and show a tendency to oscillate around the true mean values if not regularised.A second-order difference matrix is used to regularise the inverse problem.This formulation is similar to the smoothing spline method used to interpolate noisy measurements (de Boor, 2001).This paper is organized as follows: the measurement principle of OMI is described in Sect.2.1 and requirements of gridding algorithms are stated in Sect.2.2.The gridding algorithm described by Wenig et al. (2008) is recapitulated in Sect.2.3.The parabolic spline method (PSM) is explained in Sect.2.4 and its performance is evaluated using test scenarios in Sect.3. The algorithms are applied to produce maps of OMI NO 2 column densities in Sect. 4. Finally, Sect. 5 contains discussions and conclusions.

Measurement principles
In this section, we describe the measurement principle of the Ozone Monitoring Instrument (OMI), which is a passive, nadir-viewing imaging spectrometer (Levelt et al., 2006).OMI derives trace gas column densities and aerosol parameters from Earth's reflectance spectra in the near ultraviolet and visible wavelength range.The reflectance spectrum is scanned along the sunlit part of about 14.5 sun-synchronous orbits per day, resulting in daily global coverage in a moderate spatial resolution.The scanner applies the push broom principle: a charged-coupled device (CCD) array is used to perform simultaneously n (= 60) measurements in across-track direction resolving a swath width of 2600 km.In along-track direction continuous measurements are averaged every t (= 2) seconds, leading to m (≈ 1600) measurements rows in a single orbit.In each orbit, single measurements ρi,j are associated with a two-dimensional n × m lattice with across-track position i = 0, . .., n − 1 and along-track position j = 0, . .., m − 1. Introduction

Conclusions References
Tables Figures

Back Close
Full A ground pixel is the area at the earth's surface from which a certain amount of energy is received by the instrument's sensor during a single measurement.Its boundaries are typically described by the four corners of a quadrangle.Geodetic longitude and latitude for the centre of each pixel are included in OMI's Level 1b data product.The coordinates of the four corners of each ground pixel are provided by the "Ground Pixel Corners Coordinates and Sizes" product (Kurosu and Celarier, 2010), which includes two types of pixel corners: -Tiled pixel product: adjacent pixels share corners in across-and along-track direction.In along-track direction tiled pixels are always 13 km long.
-Overlapping product (FoV75): adjacent pixels share corners only in across-track direction, while in along-track direction pixels are overlapping, such that the pixel area contains 75 % of the received energy as specified by the instrument function (see below).
The tiled pixel product is useful to create clean maps, where each pixel is recognizable, while the overlapping product preserves the spatial sensitivity of the instrument better and is assumed to be more accurate for derived products, for example emission inventories.The pixel sizes in across-and along-track direction for the two pixel products are shown in Fig. 1.The pixel size of the overlapping pixel increases towards the edge of the swath, because of the wider instrument function.At the edge of the swath, overlapping pixels have about twice the length of tiled pixels.Three examples of pixel shapes of both products are shown in Fig. 2. Other products may provide different pixel products: in the Dutch OMI NO 2 (DOMINO) product (Boersma et al., 2011;Dirksen et al., 2011) ground pixels are included as a corner product which is similar to NASA's tiled product.
For each pixel, trace gas column densities ρi,j are retrieved from the measurements with measurement error εi,j .It is convenient to use across-track and along-track directions as x and y axes, respectively.In this coordinate system, the measured column density ρi,j can be expressed in terms of the true distribution ρ(x, y ) by Introduction

Conclusions References
Tables Figures

Back Close
Full ρ(x, y) W i,j (x, y) dA i,j + εi,j where A i,j is the area over which radiation is received by the sensor.The instrument function W i,j describes the spatial sensitivity of the instrument.For a push broom scanner, we can neglect the x dependency and write the instrument function as where is the convolution operator.The first term on the right hand side is the instrument's slit function and the second term describes the movement of the sensor above ground (Kurosu and Celarier, 2010).The coefficient c is chosen such that the full width at half maximum (FWHM) of the slit function at the earth's surface corresponds to an instrument's opening angle of 1 • .The FWHM depends on the distance between surface and instrument and hence the width of the instrument function increases towards the edge of the swath (Fig. 3).Since the FWHM changes are small within a pixel, the FWHM can be approximated by the mean value of each pixel.This approximation allows a computationally efficient gridding algorithm using one-dimensional splines.The second factor on the right hand side of Eq. ( 2) is the normalised boxcar function which describes the averaging in along flight direction in the interval y i,j ± 1 2 vt, where y i,j is the pixel centre, v the instrument's speed above ground and t the averaging time.The instrument function is normalized such that the integral is one.In the case of OMI, the nadir port moves roughly 13 km in along flight direction during each measurement.

Gridding algorithm
A gridding algorithm aims to obtain the best reconstruction of the true distribution of trace gas column densities ρ(x, y) from the measurement values ρi,j .The reconstructed distribution f (x, y ) is projected on a latitude-longitude grid with indices u and v , so that the value at a grid cell f u,v is given by Introduction

Conclusions References
Tables Figures

Back Close
Full where A u,v is the area of the grid cell.The resolution of regional maps for trace gas columns vary between 0.01 and 0.05 • , corresponding to 1 and 5 km at 25 • latitude.
Since on this scale A u,v is much smaller than the ground pixel area A i,j , f u,v can be approximated by the value of f (x, y) at the centre of the grid cell.
Specific requirements for a gridding algorithm may depend on the particular application.We want our gridding algorithm to have the following properties: -The weighted average of f (x, y) over the area A i,j should yield the measurement value ρi,j within the error εi,j .
-Since trace gas distributions are continuous, the reconstructed function f u,v should be continuous and have continuous derivatives.
-The shape of plumes and their maximum values should be reconstructed accurately within the measurement error εi,j .
In this paper, the term "lattice" is used to describe the ground pixel in across-and along-track direction (Level 2), while the term "grid" is used to describe the longitudelatitude grid of the Level 3 product.

Constant value method (CVM)
This section briefly describes a gridding algorithm which assumes constant values of f (x, y) within each pixel's boundaries.This constant value method (CVM) is described by Wenig et al. (2008) for OMI NO 2 column densities.CVM is used for quantitative comparison with the parabolic spline method (PSM) presented in this paper.
CVM uses a latitude-longitude grid with a resolution much smaller than the ground pixel size.Each measurement value ρi,j is stored in all grid cells lying within the pixel's Introduction

Conclusions References
Tables Figures

Back Close
Full boundaries.If pixels overlap, more than one value may be assigned to a given grid cell.In such cases, the weighted average of these values is stored in the cell.Overlap occurs when several orbits are averaged or in consecutive measurements of OMI's FoV75 pixel product.The error in the average is minimized, if the weight ω i,j is proportional to the inverse of the squared measurement uncertainty δ 2 i,j .Therefore, pixel weights are computed by with pixel size A i,j (Wenig et al., 2008).The disadvantage of CVM is that the reconstructed distribution can have discontinuities at the pixel boundaries.Furthermore, maximum values can be underestimated (see Fig. 5), because the mean value of a trace gas distribution is usually lower than the largest value within a ground pixel.This effect is stronger for pixels with a wider instrument function, because the area of the pixel is larger and thus CVM averages pixels with larger overlap.

Parabolic spline method (PSM)
In order to derive more realistic distributions, we developed a new gridding algorithm.
The trace gas distribution is described by a parabolic spline surface on the rectangular lattice formed by OMI's tiled pixel product.Parabolic splines are used, because they are the simplest, continuous functions which can preserve mean values.In this section, we first consider one-dimensional splines in across-and along-track direction, which are used to calculate two-dimensional surfaces afterwards.

Across-track direction
In across-track direction, we have n pixels, whose boundaries (x i , x i+1 ) are an increasing sequence of n + 1 knots

Conclusions References
Tables Figures

Back Close
Full with x 0 and x n being the left and right edge of the swath.A parabolic spline f (x) on the interval [x 0 , x n ] is written as the sum of piecewise-defined, parabolic polynomials f i which are defined by Each polynomial is described by three coefficients a i and a basis of parabolic polynomials Φ k .It is common (e.g.Kobza and Mlcak, 1994) to introduce the dimensionless variable s given by and to formulate f i (s) on the unit interval using the coefficients d i , p i and p i+1 defined by With these definitions, the polynomial f i (s) is with base functions

Conclusions References
Tables Figures

Back Close
Full The advantage of this formulation is that the coefficients have an obvious physical meaning as values of the distribution and thus their consistency can be validated easily.The values at the knots are given by the n + 1 coefficients p i and the mean values over the intervals are given by the n coefficients d i .Furthermore, it follows directly from Eqs. ( 8) and ( 9) that the spline is continuous.
In across-track direction, we assume that the mean values d i are given by the measurement values ρi .Therefore, we only need to determine the knot values p i , which is typically done by requiring continuity of the first derivative given by leading to the following equation On an interval with boundaries x 0 and x n , two additional boundary conditions are required to close the system.Common choices are the natural boundary conditions: To compute the n + 1 coefficients p i , the n − 1 continuity conditions (Eq.14) and the two boundary conditions (Eq.15) are written as matrix equation The parameter vector x ∈ R 2n+1 is defined as and C ∈ R with α i = 1/h i .Since the mean values d i are given by the measurement values ρ i , matrix C can be simplified to a tridiagonal matrix.The parameter vector x is computed by inverting this simplified matrix using a decomposition algorithm for a tridiagonal matrix.

Along-track direction
In along-track direction the same spline, Eq. ( 11), is used on the m +1 knots of the tiled pixel The m measurement values ρj are obtained by the convolution of the true distribution ρ(y ) with the instrument function W j (y ), see Eq.The integer r is chosen such that the integral over W j (y) contains 99 % of the total area under the curve.In this case r = 1 (i.e.integral over three pixels) at nadir and r = 2 (integral over five pixels) at the edge of the swath.Examples of instrument functions are shown in Fig. 3. Substituting Eq. ( 11) in Eq. ( 21) shows that the equation is linear in the spline coefficients p j and d j .The factors of each coefficient are computed by numerical integration of Eq. ( 21) for a given instrument function.
To compute the m coefficients d j , the m measurement equations (Eq.21) are written as matrix equation with parameter vector x ∈ R 2m+1 and data vector y ∈ R m : Equation ( 22) relates the data vector y to the weighted mean values fj using matrix M. Each row of the band matrix M has entries of the factors of the spline coefficients as calculated by Eq. ( 21).If the instrument function W j (y) is computed over three pixels (r = 1) the matrix M is as follows with base function Φ d (y) given by Eq. ( 12).The factor κ p j,k of coefficient p j+k is computed by The factors depend on the knots y j , the distance between satellite and ground pixel and the exposure time, see Eq. ( 2).They can be computed by numerical integration.

Quadratic programming problem
The parameter vector x can be computed from Eqs. ( 16) and ( 22).The inverse problem is ill-posed, if the FWHM of the instrument function is about twice the length of the pixel size.The problematic solutions have mean value coefficients d j which fluctuate around their true values, that is, the values are too small and too large, alternately.
This issue gets worse with increasing measurement uncertainty.To prevent over-and undershooting, the inverse problem of retrieving x is written as least square problem, which minimizes a quadratic objective function ϕ(x) with an additional penalty term to regularise the inversion: S ε being the covariance matrix, which here is assumed to be diagonal.The second term (on the right hand side) is the penalty term, which is weighted by the smoothing parameter γ.The matrix 2m+1) is the second-order difference matrix, which is used for the mean values d j in x only: The penalty term computes the sum of the squared deviations of d j from the local mean of three consecutive coefficients and thus prevents strong local oscillation.The diagonal matrix B is used to cancel the units of the penalty term and to reduce the dependency of the optimal smoothing on the measurement uncertainty δ ε (see Sect. 3.2).
We set the entries of B as where ρ est is an estimate of the maximum value of the distribution.
Adding the continuity constraint given by Eq. ( 16) the computation of the spline can be written as a linearly constrained quadratic programming problem: where the Hessian matrix H γ ∈ R (2m+1)×(2m+1) depends on the smoothing parameter γ.The problem 30 is rewritten using the method of Lagrangian multipliers as where λ is the vector of Lagrangian multipliers (Nocedal and Wright, 2006).The parameter vector is calculated by inverting the sparse matrix of Eq. ( 32) using a standard solver for sparse linear systems.Other strategies to solve quadratic programming problems can be found in Nocedal and Wright (2006).Introduction

Conclusions References
Tables Figures

Back Close
Full

Surface spline algorithm
A two-dimensional spline surface can be described on the nearly-rectangular lattice formed by the tiled pixel product: Following Kobza and Mlcak (1994) each piecewise defined polynomial is written as + Ψ x1 q x i,j+1 + Ψ y0 q y i,j + Ψ y1 q y i+1,j + Ψ d d i,j using nine coefficients (see also Fig. 4) defined by x i y j+1 y j f i,j (x, y)dydx and the basis functions

Conclusions References
Tables Figures

Back Close
Full If the mean values d i,j are known, the surface spline Eq. ( 34) can be computed uniquely using the one-dimensional splines described in Sect.2.4.1 by employing an algorithm described by Kobza and Mlcak (1994).Using three one-dimensional spline computations is advantageous because one-dimensional splines are calculated more efficiently than two-dimensional splines when the lattice is very large.In our application, the mean values are unknown and are computed by solving the quadratic programming problem Eq. (32) for each row (i = 0, . .., n − 1) in along-track direction using measured values ρi,j and estimated errors δ ε i,j .The resulting parameter vectors contain the coefficients q x i,j and d i,j .Thus the mean values are determined and can be used to compute the remaining coefficients.For each swath (j = 0, . .., m − 1), the mean values d i,j are used in Eq. ( 16) to compute the coefficients q y i,j .Finally, the line integrals q x i,j are used in the same way in Eq. ( 16) to compute the coefficients p i,j .Now all coefficients are known and the values of the distribution f u,v can be calculated at each grid cell.The borders of each pixel are calculated to obtain all grid cells which lie within the pixel.Next, the relative coordinates s and t are computed for the centre of each grid cell to calculate the value of the distribution f u,v = f i,j (s, t) by Eq. (34).

Handling of missing values
In observational data, trace gas values may be missing for some pixels due to unavailable measurement values or retrieval errors.In along-track direction, missing values Introduction

Conclusions References
Tables Figures

Back Close
Full could be omitted in the minimization problem Eq. ( 27) and thus the spline coefficients at these pixels would be determined by the penalty term and neighbouring pixels.However, the reconstruction of the surface spline should compute missing values using the measurement in both along-and in across-track direction.Therefore, we estimate all missing values from available measurements using bilinear interpolation.The estimated values are included with low weight in Eq. ( 27) using an (estimated) measurement uncertainty δ ε equal to the estimated maximum value ρ est .By this method the entire spline surface of an orbit can be determined.However, only the parts of the spline surface are projected on the grid for which measurement values are available.Optionally, interpolated values can be kept in the gridded data product to fill small data gaps.

Performance evaluation
In this section, test scenarios are used to study the performance of PSM and to determine the optimal smoothing parameter γ.Another important aspect addressed in this section is an analysis of the impact of measurement noise on the algorithm performances using Monte Carlo simulations.Finally, the performance of PSM is compared with CVM under different test scenarios.

Reconstruction errors
Four types of errors may occur when reconstructing a distribution.First we have to consider that the parametrization of the surface using constant values or parabolic polynomials may not be able to provide an adequate representation of the true distribution (parametrization error).Furthermore, errors may occur due to the inaccurate description of the measurement process and its numeric implementation (forward model error).Another error type is the perturbation error, which results from the propagation of the measurement error through the model.The oscillation of the mean values caused by Introduction

Conclusions References
Tables Figures

Back Close
Full measurement noise is a perturbation error.Finally, smoothing errors can occur in PSM, if the penalty term is given too much weight.We define the best possible representation as the one without (or with negligible) perturbation and smoothing errors.We use the root mean square deviation (RMSD) between true and reconstructed distribution as a measure for the overall error of the gridded result: where ρ u,v and f u,v are the true and reconstructed distributions, respectively.u and v are the grid indices and N is the number of grid cells.Furthermore, the absolute difference between maximum and minimum value is used to identify how well maxima are reconstructed by the algorithms: where u max and v max are the grid coordinates of the maximum in the true distribution.These error measures depend on the grid resolution.The continuous distribution created by PSM will only be of advantage, if the grid resolution is smaller than the ground pixel size.Here we use a grid resolution 10 % smaller than the pixel size, which is 1.3 km for OMI.

Performance for test scenarios
In this section the performances of PSM and CVM are compared for simple test scenarios and the characteristics of the different types of errors are discussed.Since PSM is based on separate one-dimensional reconstructions, we first evaluate the performance of this algorithm for the two dimensions separately and then for two-dimensional reconstructions.Introduction

Conclusions References
Tables Figures

Back Close
Full

Monte Carlo simulation
Our test scenarios consist of single plumes and a constant background distribution.
Trace gas plumes are created by localised emissions over cities or industrial areas.Constant values can be expected in remote regions without major point emissions.
The plumes are modelled by Gaussian functions with standard deviations σ between 1h and 2h (with pixel size h) which are typical sizes of plumes found in OMI's NO 2 product.The maximum value of the plumes is one.
The distributions are discretized on a one-or two-dimensional grid with x and y being across-and along-track direction, respectively.Simulated measurement values are obtained by convolving the test distribution with the instrument function, Eq. ( 1), and selecting discrete values on an equidistant lattice of 11 pixels in across-and 11 pixels in along-track direction.The position of the lattice on the test distribution is shifted between samples to mimic different instrument positions over ground.In along-track direction, instrument functions with different widths (FWHM) are tested for positions between nadir and the edge of the swath.
The measurement errors are treated as random Gaussian noise with standard deviations from 1 to 75 % of the plumes' maximum value.Our main objective is to reconstruct tropospheric NO 2 column densities, whose uncertainties are mainly below 25 % under good viewing conditions (cloud cover ≤ 50 %).This assumption is based on typical errors which are 0.1×10 16 cm −2 for clear sky and about 0.5×10 16 cm −2 for cloud fractions of 50 % (Bucsela et al., 2013).Maximum NO 2 column densities are typically about 2 to 10 × 10 16 cm −2 over polluted regions.
The performance is studied for single reconstructions and for averages of 100 reconstructions with different errors.Averaging measurements is common practise to reduce noise and to study for example annual distribution.For the different test scenarios a range of smoothing parameters γ is analysed to determine optimal smoothing.Introduction

Conclusions References
Tables Figures

Back Close
Full

Across-track direction
In across-track direction, the spline f u (x) is obtained by computing its parameters using Eq. ( 16).If a constant distribution is reconstructed, the l 2 errors are growing linearly with increasing δ ε for CVM and PSM (not shown).The reason is that forward model and parametrization errors are small and the l 2 error is dominated by the perturbation error which originates from the propagation of the measurement errors δ ε .In this scenario PSM is more sensitive to δ ε than CVM and thus l 2 errors of PSM are slightly larger.Since the measurement uncertainty δ ε is modelled as random noise, averaging k measurements will reduce the l 2 error by the inverse square root of k.Consequently, l 2 is reduced by ten for k = 100 for both algorithms.
Figure 5a illustrates the reconstruction of a plume (σ = 1.5h) using PSM and CVM for small errors.It can be seen that the reconstruction with PSM is more accurate than with CVM.The dependency of l 2 and l max on the uncertainty δ ε is shown in Fig. 6a.As for the constant distribution, the error grows with increasing uncertainty.However, in this test scenario the best possible representation is limited by the parametrization error as can be seen for small errors (δ ε < 20 %).If CVM is used, the error of the best possible reconstruction is considerably larger than for PSM.For this reason, PSM performs better than CVM for small and moderate errors.On the other hand, a larger sensitivity of PSM to δ ε causes larger errors for δ ε greater than 25 %, if PSM is used instead of CVM.The error of k averaged measurements decreases with the inverse square root of k, but converges towards a best possible representation (see Fig. 6a for k = 100).The l 2 error of CVM for k = 100 averages is smaller than the best possible representation of a single measurement.The reason is that the movement of the lattice between samples improves the parametrization of the distribution.This error reduction cannot be found for the l max error.Overall, the performance of PSM is better than CVM, if the measurement error is small or is reduced by averaging.

Conclusions References
Tables Figures

Back Close
Full

Along-track direction
In along-track direction, the coefficients of spline f v (y ) are computed by solving Eq. ( 32) for a given smoothing parameter γ.The penalty term is scaled using the diagonal matrix B. Using the scaling defined by Eq. ( 29) with ρ est equal to one, l 2 and l max have minimum values at similar γ independently of the measurement uncertainty δ ε (Fig. 7).
This suggests that the smoothing parameter γ can be selected independently of the measurement uncertainty δ ε .
For one measurement (k = 1) the performance curves shown in Fig. 7 have a distinct minimum.Larger values for γ increase the influence of the penalty term and lead to solutions which are smoothed excessively such that l 2 increases (smoothing error).
Using a smaller γ increases the sensitivity to measurement uncertainty δ ε which also causes a larger l 2 error (perturbation error).At the swath edge, where the FWHM of the instrument function is large, the l 2 error increases strongly for small γ.As discussed in Sect.2.4.2, solutions with γ equal zero at the swath edge are not usable, because they show severe over-and undershooting.Assessing the performance by l max the curve has a minimum at a smaller γ, because l max measures only the reduction of the maximum value, which is affected stronger by the smoothing error.Figure 7 shows that averaging k measurement decreases the perturbation error (for small γ) with the inverse square root of k.The systematic error of the penalty term (for large γ) is not reduced by averaging.As a result, the minimum of the l 2 and l max errors for averaged measurements is at a smaller γ than for a single measurement.Using a wider plume, smoothing errors are smaller and as a result it is possible to choose larger smoothing parameters (not shown).
Since our main interest is the reconstruction of distributions in source regions, γ is chosen such that the plume is reconstructed best for small and moderate errors.Furthermore, we want avoid systematic errors resulting from large γ.Therefore, γ is chosen such that it is slightly smaller than the minimum of l 2 for a single reconstruction.Since the position of the optimal γ depends on the pixel overlap and thus on the

Conclusions References
Tables Figures

Back Close
Full FWHM of the instrument function, this optimal γ will be different for other values of FWHM.Here, two effects are important: (1) a greater overlap increases the sensitivity to measurement uncertainties, which requires more smoothing.
(2) The same width of the plume appears narrower, if the FWHM of the instrument function is larger.A wider plume has a larger optimal γ.Monte Carlo simulations are used for a fixed width of the plume (σ = 1h) to determine the optimal γ for varying FWHM.The optimal value for γ decreases with increasing FWHM (not shown).This suggests that the second effect is dominant.However, it is more meaningful to measure the performance, when the width of the plume is directly proportional to the FWHM, since this determines the reconstructability of a distribution.Using this scenario, the optimal smoothing parameter γ increases with the FWHM of the instrument function.Figure 7 is used to pick the optimal γ at nadir (σ = 1h and FWHM about 1h) and at the swath edge (σ = 2h and FWHM about 2h) to be 1 and 10, respectively.The smoothing parameters for other FWHM are computed by linear interpolation between the optimal γ at nadir and the swath edge.Figure 5b-d illustrate the reconstruction of the plume in along-track direction using optimal smoothing parameters for which the distributions are reconstructed well.In addition, the reconstruction without regularisation (γ = 0) is shown, where over-and undershooting is indiscernible at nadir, weak near the edge and clearly visible at the swath edge (Fig. 5b-d).Hence the regularisation is mainly important for the last three rows on both edges of the swath.However, the penalty term is also suitable to reduce "fitting to noise" for all rows as is shown next.The dependency of l 2 and l max on the measurement uncertainty δ ε at nadir and at the swath edge is shown in Fig. 6b and c, respectively.l 2 and l max errors are similar to the errors in across-track direction (Fig. 6a).
The performance of PSM is improved by the penalty term, which reduces perturbation error especially for large measurement errors (δ ε > 25 %).On the other hand, the additional error due to the penalty term (smoothing error) can reduce the performance of PSM.After k = 100 averages the total error of the reconstructed plume is not a tenth of the error of a single reconstruction, because of the systematic error of the penalty term.However, PSM still performs similarly or better than CVM in all test scenarios.Figures

Back Close
Full inventory has been used in some earlier studies.The results of the simulations agree well with the observations (Huang et al., 2005;Lam et al., 2005).Figure 8a shows the modelled NO 2 column densities on 8 October 2012, 14:00 LT (local time).The distribution is normalized by its maximum value.This model distribution is sampled using instrument functions of OMI at nadir and at the swath edge.The pixel size is assumed to be 13 × 24 km 2 in along-and acrosstrack direction, respectively.At the swath edge, OMI pixels are actually 160 km long in across-track direction.However, with this coarse resolution the reconstruction of the modelled plumes would cause large parametrization errors.For this reason, we choose a smaller pixel size to better investigate the impact of the instrument function.Again, a Monte Carlo simulation is performed to study the performance of CVM and PSM for different measurement uncertainties δ ε , smoothing parameters γ and number of measurements k.The optimal smoothing parameters γ are identified as 1 and 10 at nadir and at the swath edge, respectively, which are the same as in the two-dimensional test scenarios.In Fig. 10 the l 2 and l max errors are shown for CVM and PSM using different measurement uncertainties δ ε .The results are similar to the findings obtained from the test scenarios.Figure 8 shows the reconstruction from a single measurement using a small error (δ ε = 0.05) for CVM and PSM.The structure of the plume, that is outflow towards the south west, is clearly recognisable using PSM.A cross section through the plume is given in Fig. 9.The reconstruction of the maximum value using PSM is better than for CVM.We conclude that PSM can reconstruct complex distribution over source regions with a better performance than CVM.Figures

Back Close
Full the performance without the influence of the row anomaly, our analysis concentrates on measurements before 2007.
The smoothing parameters γ is set to 1 at nadir and 10 at the edge of the swath.In between parameters are obtained by linear interpolation as function of the FWHM.The scaling factor ρ est is set to 10 16 cm −2 , which is a typical size for high NO 2 column densities.The entries of the measurement matrix (κ p and κ d ) are pre-computed for satellite-to-ground distances between 700 and 1700 km and an exposure time of 2 s using an averaged along-track pixel size of 13.5 km.Measurement uncertainties, provided by SP2, scale with NO 2 column densities due to the air mass factor dependency of the error (Wenig et al., 2008;Bucsela et al., 2013).Using these uncertainties gives lower weights to smaller densities and hence biases distributions towards low values.
To avoid this bias we adopt empirically estimated measurement uncertainties δ ε based on radiative cloud fractions (RCF) (Wenig et al., 2008): Where multiple orbits overlap, column densities at a grid cell f u,v are the weighted average of the values in all orbits.All grid cells within the boundaries of an OMI pixel are weighted by the inverse pixel area A i,j .This weighting function is suitable to preserve the continuity of column densities in the averaged products.In addition, the area dependency guarantees that individual grid cell errors are consistent with the error for the total pixel (Wenig et al., 2008).Using this weighting scheme, smaller OMI pixels are weighted higher, which enhances details on the map.Daily, monthly and annual maps of tropospheric NO 2 column densities are plotted for the Pearl River Delta region in southern China (Fig. 11).Locations of major cities and fossil-fuel power stations are marked.The daily map created by PSM shows several distinct plumes, which have larger maximum values than the plumes in the maps produced with CVM.The higher maximum values are likely to be closer to the true maximums.This was demonstrated in the performance evaluation in Sect.may differ substantially.This spatial noisiness of the distribution makes it difficult to use column amounts at a certain location, for instance for validation studies.To obtain a reasonable value, grid cells are often averaged within a certain radius around the location.This additional step is not necessary if the maps are created by PSM, because the gridding algorithm already produces a continuous distribution.As a final example, Figs. 12 and 13 show a daily and monthly distribution of NO 2 column densities over Europe and the eastern USA.

Discussion and conclusions
In this study, the parabolic spline method (PSM) was developed in order to improve the accuracy and spatial resolution of trace gas maps obtained from satellite measurements.The new gridding algorithm reconstructs the trace gas distribution by a continuous parabolic spline surface using measurements from a single orbit.The spline coefficients are computed using one-dimensional splines, which allow fast computation for a large number of pixels.PSM is the first gridding algorithm that explicitly treats the spatial sensitivity distribution of the instrument over ground using an instrument function.If the FWHM of the instrument function is about twice the length of the tiled pixel, the inverse problem is very sensitive to measurement errors.Therefore, a penalty term is used to regularise the inversion by a second-order difference matrix.This regularisation is similar to the penalty term in a smoothing spline and hence is also used to smooth measurement noise (de Boor, 2001).Monte Carlo simulations are conducted to study the performance of the algorithm for different distributions of trace gas column densities.The optimal weight of the penalty term is found to be proportional to the measurement uncertainty and the width of the instrument function.
The algorithm calculates the unknown mean values over the tiled pixel product from the measurements.While we use these values as spline coefficients, it should be noted that they can also be used directly as a better approximation to visualize the mean values over the tiled pixels.This approach is more accurate than the often used as-Introduction

Conclusions References
Tables Figures

Back Close
Full sumption that the mean values over the tiled pixel are equivalent to the measurement values.A related study on the best approximation of mean values over tiled pixels is planned in the future.The performance of PSM is compared with a standard algorithm, described by Wenig et al. (2008), which assumes constant measurement values within pixel boundaries and thus is named constant value method (CVM) here.Due to the smaller parametrization errors, reconstructions using PSM are better than CVM, if the perturbation error is small or is reduced by averaging.Trace gas maps created by PSM are smoother and maximum values are better reconstructed.Better performance is also illustrated by reconstructing the high-resolution distribution calculated with the CMAQ chemistry model.Application to OMI NO 2 column densities suggests that PSM is more suitable than CVM to create high-resolution maps.The use of PSM in satellite validation studies is being planned for in the future.
The general approach described in this study can be easily extended for higherorder splines, for instance quartic splines, to further improve the parametrization of the distribution.The gridding algorithm can also be applied to other OMI-like instruments such as the TROPOspheric Monitoring Instrument (TROPOMI), planned for launch in 2014 (Veefkind et al., 2012).The high spatial resolution of TROPOMI (7×7 km 2 at nadir) allows for the detection of small-scale features and is likely to reduce the measurement uncertainty due to more cloud-free observations.
In this paper, we demonstrate that the newly developed gridding algorithm improves the reconstruction of regional trace gas maps.Its application could be very helpful when satellite-derived trace gas distributions are studied on a regional scale.Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | (n+1)×(2n+1) is a band matrix of the form Introduction Discussion Paper | Discussion Paper | Discussion Paper | (2) which now takes the form ρj = +∞ −∞ ρ(y)W j (y )dy + εj (20) with measurement error εj .Likewise, the weighted mean value fj can be derived from the spline function f (yl (y) W j (y ) dy Discussion Paper | Discussion Paper | Discussion Paper | is the factor of coefficient d j+k calculated by κ d j,k = y j+k +1 y j+k Φ d (y )W j (y )dy (25Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Ψ 11 (s, t) = st(4 − 6s − 6t + 9st) Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3. The monthly map for October 2006 produced with CVM is coarse and values of adjacent grid cells Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 5 .Fig. 6 .
Fig.5.Reconstruction of a one-dimensional plume (σ = 1.5h) using CVM and PSM in (a) across-track direction and in along-track direction for instrument functions (b) at nadir, (c) near the swath edge and (d) at the edge of the swath.The reconstruction using PSM is shown with optimal and without (γ = 0) regularisation.The measurement values ρ j (with δ ε = 5 %) are drawn at the centre of each pixel.