Pyplis — A Python Software Toolbox for the Analysis of SO 2 Camera Images for Emission Rate Retrievals from Point Sources

Ultraviolet (UV) SO2 cameras have become a common tool to measure and monitor SO2 emission rates, mostly from volcanoes but also from anthropogenic sources (e.g., power plants or ships). Over the past decade, the analysis of UV SO2 camera data has seen many improvements. As a result, for many of the required analysis steps, several alternatives exist today (e.g., cell vs. DOAS based camera calibration; optical flow vs. cross-correlation based gas-velocity retrieval). This inspired the development of Pyplis (Python plume imaging software), an open-source software toolbox written in Python 2.7, which unifies the most prevalent methods from literature within a single, cross-platform analysis framework. Pyplis comprises a vast collection of algorithms relevant for the analysis of UV SO2 camera data. These include several routines to retrieve plume background radiances as well as routines for cell and DOAS based camera calibration. The latter includes two independent methods to identify the DOAS field-of-view (FOV) within the camera images (based on (1) Pearson correlation and (2) IFR inversion method). Plume velocities can be retrieved using an optical flow algorithm as well as signal cross-correlation. Furthermore, Pyplis includes a routine to perform a first order correction of the signal dilution effect (also referred to as light dilution). All required geometrical calculations are performed within a 3D model environment allowing for distance retrievals to plume and local terrain features on a pixel basis. SO2 emission rates can be retrieved simultaneously for an arbitrary number of plume intersections. Hence, Pyplis provides a state-of-the-art framework for more efficient and flexible analyses of UV SO2 camera data and, therefore, marks an important step forward towards more transparency, reliability and inter-comparability of the results. Pyplis has been extensively and successfully tested using data from several field campaigns. Here, the main features are introduced using a dataset obtained at Mt. Etna, Italy on 16 September 2015.


Introduction
Sulfur dioxide (SO 2 ) is a toxic gas emitted by anthropogenic and natural sources (e.g., power plants, ships, volcanoes).The pollutant has impacts on the atmosphere both on local and global scales (e.g., particle formation, radiation budget, e.g., [1,2]).Furthermore, the monitoring of SO 2 emissions from active volcanoes can provide insight into the magmatic degassing behaviour and hence plays an important role for the development of new risk assessment approaches (e.g., [3][4][5][6], and references therein).
The gas composition of the emission plumes can, for instance, be studied using ground-based passive remote sensing techniques.The column-densities (CDs) of the gases in the plumes are quantified based on the absorption signature carried by scattered sunlight that has penetrated the plume.SO 2 -CDs, for instance, can be retrieved at ultraviolet (UV) wavelengths (i.e., around 310 nm) where it exhibits distinct absorption bands.Prominent examples for passive remote sensing instrumentation are the correlation spectrometer (COSPEC, [7]), or instruments based on the technique of Differential Optical Absorption Spectroscopy (DOAS, [8], e.g., [9,10]).Over the past decade, the comparatively young technique of UV SO 2 cameras has gained in importance, since it enables the study of volcanic SO 2 emissions at unprecedented spatial and temporal resolution (e.g., [11][12][13][14][15]).This is particularly helpful to study multiple sources independently (e.g., [16]) or to investigate volcanic degassing characteristics by studying periodicities in the SO 2 emission rates (e.g., [17][18][19] and references therein).The technique of UV SO 2 cameras has seen remarkable improvements over the past decade (e.g., [20][21][22][23][24]) and can now be considered one of the standard methods for ground-based remote sensing of SO 2 plumes.A drawback, however, is the low spectral resolution, restricting the technique to a single species and furthermore requiring external calibration.
The retrieval of SO 2 emission rates from plume imagery comprises several analysis steps that are summarised in Table 1 and are illustrated in the flowchart shown in Appendix Figure A1.Thanks to ongoing developments, today, researchers can choose between several methods for nearly all of the required steps (e.g., cell vs. DOAS calibration, velocity retrieval using optical flow vs. cross-correlation method).
Available software solutions include Vulcamera [25], the IDL R source code provided by [22] and Plumetrack [26].The first two programs include routines for cell calibration and cross-correlation based plume velocity retrievals.The IDL R program also includes a routine to perform a first order correction for the signal dilution effect (commonly referred to as light dilution, e.g., [22]).The software Plumetrack provides an interface to calculate gas velocities using an optical flow algorithm and can be applied to pre-calibrated SO 2 -CD images in order to retrieve SO 2 emission rates.Pyplis ( [27,28]) is a cross-platform, open-source software toolbox for the analysis of UV SO 2 camera data.The code is written in Python 2.7 and emerged from the idea of a common software platform incorporating the most relevant analysis methods, including recent developments.The most important features of Pyplis 1.0.0 are (details follow in Section 3): 1. 3D distance retrievals to plume and local terrain features at pixel-level, 2. several methods to retrieve plume background radiances, 3. cell and DOAS based camera calibration including two independent DOAS FOV search routines, 4. cross-correlation and optical flow based plume velocity retrievals, 5. histogram based correction for ill-posed optical flow vectors in low-contrast image regions, 6. image based correction for the signal dilution effect, 7. automated emission rate retrievals along linear plume intersections.Pyplis comes with numerous example scripts providing an easy and comprehensive introduction into the software.The following Section 2 introduces the technique of UV SO 2 cameras and the required analysis steps for SO 2 emission rate retrievals.The implementation of the individual analysis methods is discussed in Section 3.

UV SO 2 Cameras
UV SO 2 cameras analyse scattered sunlight that has penetrated a plume containing SO 2 gas.Plume optical densities (ODs) τ are retrieved in two wavelength windows of approximately 10 nm width, using optical bandpass filters.One filter is centred around 310-315 nm, where SO 2 has strong absorption bands (referred to as SO 2 on-band).A second filter is situated around 330 nm, where SO 2 absorption is weak (SO 2 off-band).The latter is used to correct for additional broadband light extinction, for instance, resulting from aerosols or water droplets in the plume.From the retrieved ODs in both channels, the apparent absorbance (AA) of SO 2 (τ AA , e.g., [20]) can be retrieved as where I and I 0 denote the plume and plume background radiances, respectively, in both channels (on, off).Note that the method requires all additional optical densities in the on and off-band regime to be of broadband nature.

Image Analysis-Retrieval of S SO2 Images
AA images are determined from a set of pre-processed (e.g., dark / offset corrected) plume and background images using Equation (1).Next, the AA images are converted into SO 2 column-density (CD) images, where denotes the SO 2 -CD in the viewing direction C ij of the image pixel i, j. c(x, y, z) is the concentration distribution of SO 2 in real world coordinates x, y, z (cf. Figure 1a) and ds = dx 2 + dy 2 + dz 2 is the integration differential.The camera calibration (i.e., the conversion of AA values into SO 2 -CDs) can be performed either using gas cells (e.g., [11]) or a co-located DOAS spectrometer (e.g., [21]).The latter is more accurate in case aerosols are present in the plume [21,29].The position and shape of the DOAS FOV within the camera images are required in order to perform the DOAS calibration.The DOAS FOV can either be measured experimentally (e.g., in the lab) or can be retrieved directly from the field data [21,30].
The camera calibration curve depends on the pixel-position in the images (e.g., [21]).This is due to shifts in the filter transmission curves for non-perpendicular illumination (i.e., off-axis image pixels), resulting in an increased SO 2 sensitivity towards the image edges (e.g., [31]).A first order correction for this effect can be performed using a normalised sensitivity correction mask retrieved from cell calibration data, as described in [21].The impact of this effect can be dramatically reduced by placing the transmission filters between lens and detector rather than in front of the lens [31].Please note that this effect of an increased SO 2 sensitivity towards the image edges is not to be confused with lens-vignetting effects, as discussed by [31].

Emission Rate Retrieval
SO 2 emission rates Φ are retrieved from SO 2 -CD images by performing a discrete path integration along a suitable plume cross section (PCS) projected into the image plane, for instance a straight line (illustrated in Figure 1b).Then, corresponds to the SO 2 emission rate through , where m denotes one of a total of M sample positions along in the image plane and ∆s is the integration step length, measured in physical distances on the detector.f is the focal length of the camera, d pl is the distance between camera and plume and n is the normal of .v is a 2-vector containing projected plume velocities averaged in the viewing direction.
The plume distances d pl can be estimated from the measurement geometry and require information about the camera position and viewing direction as well as the source position and the meteorological wind direction.The gas velocities are usually retrieved from the images directly either using cross-correlation based methods (e.g., [32]) or optical flow algorithms (e.g., [23]).
Equation ( 3) is equivalent to commonly used retrieval methods (e.g., [11,29]) and is based on the assumption that over or underestimations of the measured quantities S SO2 and v (e.g., due to non-perpendicular plume transects) cancel out in the emission rate retrieval.This is a valid approximation for typical measurement conditions (i.e., plume nearly perpendicular to the optical axis, and plume extent small compared to plume distance).However, care has to be taken for unfavourable geometries (e.g., tilted or overhead plume; retrieval close to the image edges), which may require additional corrections (e.g., [33]).The software Plumetrack [26] includes an alternative "2D" method to compute the emission rates, which considers all image pixels that are crossing the PCS line between consecutive frames.This method is not included in the current version of Pyplis but may be a valuable extension in a future release of the software.
If a locally uniform gas velocity can be assumed (i.e., v(i, j) → v) and a planar PCS is used for the retrieval (i.e., n(i, j) → n), then, Equation (3) can be further approximated as denotes the integrated column amount (ICA) along .

Radiative Transfer Corrections
Radiative transfer effects may introduce systematic errors in the retrieved emission rates (see e.g., [34,35]).The signal dilution effect describes a decrease in the measured CDs due to scattering between instrument and plume.The magnitude of this effect primarily depends on the local visibility (i.e., amount of molecules and particles in the ambient atmosphere) and on the distance between camera and plume.A first order correction can be performed using the atmospheric scattering extinction coefficients in the viewing direction between camera and source (e.g., [12,36]).The latter can be retrieved, for instance, by studying brightness variations of topographic terrain features in the images as a function of their distance to the instrument [22].More complicated radiative transfer issues (e.g., optically thick plumes) require corrections using radiative transfer models (e.g., [29]).Pyplis can correct for the signal dilution effect based on the method suggested by [22].

Implementation
In the following, the main features and modules of Pyplis are presented.The application programming interface (API) is designed to be modular using an object oriented architecture.It can therefore be used in parts or as a whole.Figure 2 illustrates the Pyplis API for SO 2 emission rate retrievals (see also Appendix E).Pyplis includes 19 example scripts, which provide an introduction to the main features of the software (summarised in Table A2).The scripts are based on an example dataset recorded at Mt. Etna in Italy on 16 September 2015 between 6:40 a.m. and 7:30 a.m.UTC.These data are freely available and can be downloaded from the website (for details, see Appendix F.1). include relevant meta information and can be used to create Dataset objects (blue).The latter perform file separation by image type and create ImgList objects (green) for each type (e.g., on, off, dark).Further analysis classes are indicated in yellow.Note that the routine for signal dilution correction is not shown here (cf. Figure A1).

Geometrical Calculations
Geometrical calculations are performed within an instance of the MeasGeometry class (part of the geometry.pymodule).The plume distances d pl (cf.Equation ( 3)) can be retrieved per pixel-column i using intersections of the respective viewing azimuth with the plume azimuth (see Figure 3).This requires specification of the camera position, viewing direction and optics (e.g., detector specifications, focal length), source coordinates and meteorological wind direction.The distances are determined based on the horizontal plume distance and the altitude difference between source and camera (i.e., assuming a horizontally aligned plume).It is pointed out that this approach may be inapplicable for complicated measurement geometries (e.g., overhead plumes), which would require a more detailed knowledge of the plume shape and altitude.The map includes plume orientation (dark green line), camera azimuth retrieved using a compass (dashed light green line), the corrected camera azimuth (light green line) and the corresponding FOV (semi-transparent green), retrieved automatically using the pixel-position i, j of the Etna south-east crater within the images and the corresponding coordinates of the crater (longitude, latitude, altitude).The map was generated using an instance of the MeasGeometry class (see Section 3.1).
Further features of the MeasGeometry class include a routine to retrieve the camera viewing direction based on the position of distinct objects in the images (e.g., summit of volcano, cf. Figure 3), or the calculation of distances to topographic features in the images (cf. Figure 10).The MeasGeometry class makes use of the Geonum library [37], which is briefly introduced in Appendix A, and which provides automatic online access to topographic data from the NASA shuttle radar topography mission (SRTM, [38]) and also supports the ETOPO1 dataset [39].Furthermore, 2D and 3D topographic overview maps of the measurement setup can be created automatically (as shown in Figure 3, see also Figure 10c).

Image Representation and Pre-Processing Routines
Images are represented by the Img class (image.pymodule).The Img class includes reading routines for all image formats supported by the Python Imaging Library (PIL, e.g., png, tiff, jpeg, bmp) as well as the FITS format (Flexible Image Transport System).It further allows to store relevant meta information (e.g., exposure and acquisition time, focal length, f-number) and includes several basic processing methods (e.g., smoothing, cropping or resizing using a Gaussian pyramid approach).The image data itself is loaded and stored as a 2D-NumPy array within an Img object.Automatic dark and offset correction can be performed as described in Appendix C. Img objects can be saved as FITS files at any stage of the analysis (e.g., dark corrected image, τ AA image, calibrated SO 2 -CD image).

Retrieval of Plume Background Radiances
The calculation of the OD images in both wavelength channels requires the retrieval of the sky-background intensities I 0 behind the plume (cf.Equation ( 1)).The plumebackground.pymodule provides several alternatives to retrieve the background intensities, either from the plume images directly or using an additional sky reference image (I 0 -image).For the latter, several methods are available to correct for offsets and non-uniformities in the sky-background between the plume and I 0 -image.The corrections are based on suitable sky-reference-areas in the plume image (rectangles or profile-lines, cf.Table 2) and use polynomials of a suitable order (e.g., linear or quadratic) to model the I 0 -image such that the corresponding OD image satisfies τ = 0 within the specified sky reference areas.
If no I 0 -image is available, the plume background radiances can also be estimated from the plume images directly using a masked 2D polynomial surface fit.The required mask specifies clear-sky image pixels that are considered during the fit.The mask can either be provided by the user or can be retrieved automatically using the function find_sky_background of the plumebackground.pymodule.
The PlumeBackgroundModel class (part of plumebackground.py)provides high level access to eight default methods for the background retrieval (in the code denoted with mode).These different retrieval methods are summarised in Table 2. Figure 4 shows four example plume on-band OD images calculated using the background modelling methods 0, 1, 4 and 6 (cf.Table 2).It is not intended to give a recommendation for a "best" method here, as this strongly depends on the data (e.g., availability of suitable reference areas; acquisition time and relative viewing direction of I 0 -image; solar position).

Camera Calibration
The camera calibration can either be performed using data based on cuvettes (gas cells) filled with a known amount of SO 2 -gas, or using plume SO 2 -CDs retrieved from a co-located DOAS spectrometer (cf.Table 1).

Calibration Using SO 2 Cells
Cell calibration can be performed using the CellCalibEngine class (cellcalib.pymodule) and requires a dataset containing both cell and sky background images.It is assumed that the camera is pointed into a gas and cloud free area of the sky and that the cells (containing different SO 2 -CDs) are consecutively placed in front of the lens, such that they cover the whole FOV of the camera.Figure 5a shows a time-series of the image mean intensities retrieved from such a dataset.The individual cells can be identified by abrupt intensity drops in the time-series while the ambient background only changes gradually.The corresponding time-intervals for cell and background images can be specified by the user or, alternatively, detected automatically within an instance of the CellCalibEngine class (shown in Figure 5a, see also example scripts 0.7 and 5).Cell OD images in both channels (τ on , τ off ) can be determined using a suitable background image.Care has to be taken for measurements performed at large solar zenith angles (early morning, late afternoon) due to rapid changes of the ambient sky intensity (cf. Figure 5a).In this case, the background image needs to be scaled to the sky intensity present at the acquisition time of each cell in order to calculate the OD images.This correction was performed for the data shown in Figure 5a (i.e., dashed line "Fitted BG polynomial").It requires at least two background images, one recorded before, and a second one after the cells were put in front of the lens.
The calibration results (e.g., one AA image for each cell) are stored within CellCalibData objects together with the corresponding cell SO 2 -CDs (which need to be provided by the user).Calibration curves can then be retrieved per image pixel or within a certain pixel neighbourhood.Figure 5b shows example calibration curves for τ on , τ off and τ AA images.

Calibration Using DOAS Data
The DOAS calibration is performed using a set of plume optical density images (usually AA images) and a corresponding time-series of SO 2 -CDs retrieved from a DOAS spectrometer.In a first step, the AA images are stacked into a 3D-NumPy array and merged in time with the DOAS data.The latter can be performed in three different ways (for details see Appendix D.1).The calibration data (i.e., merged time-series of SO 2 -CDs and camera AA values within the DOAS FOV) can be retrieved by convolving the AA stack with a mask specifying position and shape of the DOAS field-of-view (FOV) within the image plane.The calibration data is stored within instances of the DoasCalibData class, which is also used to retrieve the calibration curve.The latter is done by fitting a polynomial of appropriate order to the calibration data.Optionally, the fit can be performed using a weighted regression, to account for statistical uncertainties in the DOAS SO 2 -CDs (e.g., fit-errors, cf. Figure 6).DoasCalibData objects can be stored using the FITS standard (including the FOV mask).
Note that the DOAS calibration curve is only valid within the image pixel area covered by the DOAS FOV.This is due to cross-detector variations in the SO 2 sensitivity (see Section 2.2) and can be corrected for using a mask retrieved from a cell calibration dataset.The mask is determined by fitting a 2D polynomial to a cell AA image (see prev.Section 3.4.1),which is then normalised to the centre-position of the DOAS FOV (illustrated in example script 7).Please also note that Pyplis cannot perform the DOAS analysis.Thus, the SO 2 -CDs need to be retrieved using a suitable 3rd party software (e.g., DOASIS, [40]).7. The fit was performed using a first-order weighted polynomial fit.The weights were calculated using the relative errors ∆S SO2 / S SO2 of the DOAS SO 2 -CDs.The y-axis offset is likely due to uncertainties in the DOAS retrieval (e.g., due to O 3 interference) and is compensated by the calibration.

DOAS FOV Search
Pyplis includes two routines to retrieve the DOAS FOV mask (included in the DoasFOVEngine class) based on a stack of AA images and a DOAS data vector: 1. Pearson routine: this method loops over all image pixels in the AA stack and determines the Pearson correlation coefficient between the corresponding AA time-series (τ AA (t)) and the DOAS SO 2 -CD vector (S SO2 (t)).The method yields a correlation image as shown in Figure 7a, from which the pixel coordinate with highest correlation (i 0 , j 0 ) is extracted (see also [21]).Assuming a circular FOV shape, the pixel extent of the FOV is estimated around i 0 , j 0 , by iteratively searching the disk radius with highest correlation between the AA and the DOAS time-series.2. In-operation field-of-view retrieval (IFR) routine: this method is based on [30] and uses an inversion algorithm to retrieve the FOV.Position and shape of the FOV is parametrised by fitting a 2D Super-Gaussian to the retrieved FOV distribution (shown in Figure 7b, see Appendix D.2 for details).
The retrieved FOV information (position, shape, convolution mask) is represented by the DoasFOV class and can be stored as a FITS file.

Plume Velocity Analysis
Plume velocities can be retrieved either using the ICA cross-correlation method or using an optical flow algorithm.Both methods yield displacement estimates in units of pixels / time.These are converted into plume gas velocities based on the measurement geometry (MeasGeometry, see Section 3.1).The relevant code is implemented in the plumespeed.pymodule.The retrieved FOV positions show good agreement.Note that the FOV was retrieved from downscaled images (Gauss pyramid level 2).

Velocity Retrieval Using the ICA Cross-Correlation Method
For the cross-correlation method, ICA time-series (see Equation ( 4)) are determined using two PCS lines located at different positions downwind the emission source.Ideally, the PCS lines should be parallel and should cover an entire plume cross-section (indicated in Figure 8, left).
In a first step, the two time-series are re-sampled onto a regular grid (default is 1 s resolution).In a second step, a correlation analysis is applied to the re-sampled data vectors in order to find the time lag corresponding to the highest correlation between both signals.The method yields one average velocity vector, which is representative for the used image region and time-series.Note that the method intrinsically assumes that the average plume propagation direction φ in the i, j-plane is aligned with the PCS normal (i.e., φ n) and furthermore that SO 2 is conserved between the two PCS lines used.Figure 8 shows results from an example cross-correlation analysis, resulting in a plume velocity of 3.5 m/s.

Optical Flow Based Velocity Retrievals
Optical flow velocity retrievals are performed using the Farnebäck algorithm [41], which is implemented in the OpenCV library [42].The algorithm can estimate motion at the pixel-level by tracking local contrast features between consecutive frames.Note, however, that optical flow algorithms may fail to detect motion in extended homogeneous image areas.In this case, appropriate corrections may be required (e.g., [43]).
All relevant calculations for optical flow based velocity retrievals are performed within the OptflowFarneback class.The class includes the Farnebäck algorithm itself as well as several pre-and post-analysis routines.The latter includes a routine that can detect and correct for unphysical motion estimates in homogeneous image regions based on the method proposed by [43].The correction identifies the local average velocity vector using peaks in histograms of an optical flow motion field.The Farnebäck algorithm itself requires specification of several input parameters (see e.g., [23], or OpenCV documentation).The Pyplis default settings follow the suggestions of [23].An example flow field is shown in Figure 9 including results from the histogram post analysis within a narrow region-of-interest (ROI) around an example retrieval line.The latter yields an average velocity magnitude of 4.2 (± 0.4) m/s and a predominant movement direction of ϕ = −65 (± 14) • within the image plane.Optical flow plume velocity retrievals, including the histogram post analysis method, are introduced in example script 9 (see Table A2).

Image Based Signal Dilution Correction
A correction for the signal dilution effect (see Section 2.4) can be performed using the DilutionCorr class.The method is based on [22] and uses the model function to retrieve atmospheric extinction coefficients in the on and off-band regime.Here, λ denotes the wavelength, I meas are measured intensities of dark terrain features within the images, d is the distance to these features and I 0 their initial intensity.The ambient intensity I A can be approximated using gas free sky areas in the plume images [22].For the retrieval, a set of measured intensities I meas is extracted from vignetting corrected plume images using suitable terrain features (e.g., pixels along a volcanic flank).These intensities are fitted to Equation ( 5) (as a function of their distance d) in order to retrieve the extinction coefficients in both wavelength regimes.The required distances d to the features can be retrieved automatically, based on intersections of individual pixel viewing directions with the local topography (using the Geonum library; for details, see Appendix A).The plume images can then be corrected for the signal dilution effect using the retrieved extinction coefficients (see Equation ( 4) in [22]).The correction is only applied to plume pixels (e.g., using a τ AA threshold mask).The required plume distances can be retrieved from the measurement geometry (for details, see Section 3.1).
Figure 10 shows results of an example dilution correction using an on and off-band image from the Etna example data, recorded at 6:45 a.m.UTC.Extinction coefficients of on = 0.0743 km (Figure 10a) and off = 0.0654 km −1 (Figure 10b) were retrieved using the two topographic profile lines shown in Figure 10c,d

Emission Rate Retrieval
Emission rates can be determined using the EmissionRateAnalysis class.The analysis is performed based on an image list containing calibrated SO 2 -CD images (see also Appendix E.2) and a specification of one or more retrieval PCS lines (LineOnImage objects, cf. Figure 2).Plume velocities can be provided by the user (e.g., using the result of a cross-correlation analysis, see Section 3.5.1)or can be retrieved automatically during the emission rate analysis using the Farnebäck optical flow routine (see Section 3.5.2).For the latter, three options are available: 1. flow_raw → the raw output of the Farnebäck algorithm is used.This should only be done if it can be assured that the algorithm yields reliable output in the considered plume area (i.e., ROI around the PCS line) and for all images of the time-series. 2. flow_histo → performs the histogram post-analysis proposed by [43] (cf.Section 3.5.2).The retrieved local average velocity vector for each PCS line is then used as a velocity estimate for the corresponding retrieval line.3. flow_hybrid → reliable motion vectors from the flow field are used while unreliable ones are identified and replaced based on the histogram post-analysis (see previous point).
Figure 11 shows results from an emission rate analysis for the Etna test data (upper panel) including plume velocities retrieved using the flow_hybrid method (lower panel).The AA images were calculated from dilution corrected on and off-band plume images using background modelling method 6 (cf.Section 3.3) and were corrected for cross-detector sensitivity variations using a mask retrieved from cell "a53" (cf. Figure 5a).The AA images were calibrated using the Pearson DOAS calibration curve shown in Figure 6.(3) using the flow_histo method (thin dashed brown) and (4) using the flow_hybrid method (bold orange).The retrieved effective velocities are plotted in the lower panel and correspond to the average velocities along the PCS line using the flow_hybrid method.Uncertainties are displayed as shaded areas.

Remark on Performance
The emission rate analysis (cf. Figure 11) was performed using a Gaussian pyramid level of 1 (i.e., reduction of image size by a factor of 2).For these settings, the typical computation time to calculate a calibrated SO 2 -CD image from the raw data amounts to approximately 0.2 s using an Intel(R) Core(TM) i7-6500U CPU (2.50 GHz, 8 GB RAM) (Santa Clara, CA, USA).Approximately the same time of 0.2 s is required to compute the optical flow between two frames at this pyramid level (cf.Appendix Table A1).For the discussed dataset of 209 on-band plume images, this results in a total computation time of about two minutes, including the four different velocity retrieval methods introduced above and shown in Figure 11.Note that the latter does not include the time to compute the DOAS FOV and calibration data, nor the time required for the cross-correlation velocity retrieval, since these analyses were performed beforehand (for details, see Appendix B and example script 12).

Remark on Uncertainties
The uncertainties in the emission rates (cf.shaded areas in Figure 11) were calculated using Gaussian error propagation (see method det_emission_rate of the fluxcalc.pymodule) considering (1) uncertainties in the meteorological wind direction and camera viewing direction (accessible using the method plume_dist_err of the MeasGeometry class); (2) uncertainties in the calibration curve (based on slope error retrieved from covariance matrix of fit result, accessible using the attribute slope_err of either the DoasCalibData or the CellCalibData class); and (3) the uncertainties in the velocity retrieval.For optical flow based velocity analyses, the uncertainties are calculated as discussed in Section 2.4.1 in [43], based on the retrieved effective velocities v eff = v • n along (cf.Equation ( 3)).Note that Pyplis does not provide a method to estimate the uncertainties of cross-correlation based velocity retrievals.These need to be provided by the user (cf.example script 12).The errors in the effective velocities can be accessed via the velo_eff_err attribute of the EmissionRates class.Note that all uncertainties computed by Pyplis are assumed to be of statistical nature.Hence, they do not account for any potential systematic errors, for instance arising from ill-constrained optical flow vectors when using the flow_raw retrieval method (see discussion in Section 2.4.1 in [43]).

Conclusions
In this paper, the software package Pyplis was introduced.Pyplis contains an extensive collection of relevant algorithms for the analysis of UV SO 2 camera data, particularly for the retrieval of emission rates from SO 2 -emitters (e.g., volcanoes, power plants, ships).
Apart from established analysis methods, such as cross-correlation based velocity retrievals (e.g., [32]) or cell and DOAS calibration (e.g., [21]), Pyplis incorporates more recent developments.These include an implementation of the DOAS FOV retrieval algorithm proposed by [30], a routine to correct for the signal dilution effect based on [22], or pixel based gas velocity retrievals using an optical flow algorithm (e.g., [23]).The latter incorporates a method to detect and correct for ill constrained optical flow motion-vectors based on [43].Furthermore, Pyplis includes a framework to perform a detailed 3D-analysis of the observed camera scene including automatic online access to high resolution topography data from the SRTM dataset.This enables, for example, to retrieve the camera viewing direction based on distinct topographic features in the images (e.g., volcano summit), or to calculate distances between the camera and the local topography at a pixel-level.The latter is of particular relevance for the image based correction of the signal dilution effect.
Due to this extensive collection of algorithms, Pyplis provides flexibility with regard to the analysis strategy and is highly adaptable to different data situations.The object oriented architecture of the API gives intuitive access to the main features and makes it easy to compare individual methods (e.g., different plume velocity retrievals, as illustrated in this paper).Pyplis is open-source and can be operated both on Windows and Unix machines.Thus, it is well suited for inter-comparison studies, the exchange of analysis results or for the development and verification of new methods.
The Pyplis installation includes numerous example scripts that were used in this paper to introduce the main features of the software.The examples are based on a 15 min dataset recorded at Mt. Etna, Italy in September 2015, which is freely available and can be downloaded from the website.
Finally, the authors wish to point out that Pyplis may also be used for other applications based on the same measurement principle (e.g., NO 2 cameras) and that parts of it can also be useful for other remote sensing applications (e.g., the engine for geometrical calculations).The Pyplis software is hosted on GitHub ( [27,28]).The code documentation and further information (e.g., installation instructions) can be found on the documentation website (see [27,28]).
Table A1 summarises typical computation times for these steps that are performed for each image during the emission rate analysis.Further settings for this performance analysis (e.g., background modelling method) correspond to the settings provided in Section 3.7.Note that here we omit the discussion of potential additional analysis procedures that may be required and that are typically performed prior to the emission rate analysis (e.g., dilution correction, DOAS and / or cell calibration, cross-correlation velocity retrieval, cf. Figure 2).These highly depend on the chosen settings and often include looping over a number of images and pre-computation of OD images (e.g., DOAS FOV search, cross-correlation velocity analysis) and we encourage the reader to use the corresponding example scripts for individualised performance tests of these additional analysis procedures.

• Second Method: Vice Versa Interpolation of Both Grids
This method uses the unified sampling grid (all time stamps from both datasets) and performs interpolation of the DOAS data vector (at image acquisition time stamps) and vice versa.
The method is slow compared to method 1 since each image pixel of the stack is interpolated.However, it results in more data points, which can be an advantage especially for short time series.This method can be significantly accelerated by reducing the image size or by only performing the analysis within a certain image region (c.f.example script no.6, Table A2, script option: DO_FINE_SEARCH).The time series interpolation is done using the pandas Python library.

•
Third Method: Nearest Data Point This method loops over all spectra and for each spectrum, finds the image which is nearest in time.This method is for instance used, if only the acquisition time stamps are provided and not the start / stop stamps of each exposure (which is required for the first method).

Appendix D.2. FOV Determination Applying the IFR Method
The In-operation Field of view Retrieval (IFR) method is an implementation of the method proposed by [30].IFR applies a linear camera model to invert the FOV of a low-resolution instrument (in this case, a DOAS spectrometer) from imager data without a priori information (e.g., FOV position, size and shape).The inversion problem is typically under-determined for SO 2 camera applications.Therefore, the iterative LSMR method [44] is applied to retrieve the FOV coefficients depending on the regularisation parameter λ.
The choice of λ is somewhat arbitrary but may influence the IFR results depending on the input data quality.The default value is λ = 10 −6 .However, in case only a small sample size is available, λ may need to be increased (e.g., λ = 10 −3 ) in order to produce meaningful results.A side effect of increasing λ is a spatial smoothing of the results potentially leading to unrealistic large FOVs. Figure 7b shows a sample FOV distribution retrieved from the Etna test data (88 images) using λ = 2 × 10 −3 .
In order to reach the final goal of gaining a calibration curve from the image stack containing AA images, individual images need to be convolved with the FOV of the low-resolution instrument.Therefore, a parametrised FOV is fitted to the IFR results, which is more applicable to ground-based instruments than the parametrisation proposed by [30].We propose the following elliptical Super-Gaussian FOV parametrisation g of the IFR result depending on pixel coordinates i, j in horizontal and vertical direction, respectively: where C is a constant offset, A the amplitude, i m and j m define the centre position, σ measures the width in i-direction, the asymmetry parameter a measures the ratio of the semi-axes (e.g., a = 1 yields a circular FOV), and b is the shape parameter of the Super Gaussian (e.g., b = 1 yields a Gaussian FOV, b = 10 approximates a flat-disk FOV).
If an additional tilt of the FOV is required in case of an elliptical FOV, the above fit may be performed in a transformed coordinate system and i defines the transformation into tilted coordinates i and j .Equation (A2) is then replaced by • EmissionRateAnalysis (fluxcalc.py)→ Performs emission rate analysis based on an ImgList containing calibrated images.Emission rates can be retrieved along one (or more) plume cross section lines (LineOnImage objects) and has several options related to the plume velocity retrieval (Section 3.7).

•
EmissionRates (fluxcalc.py)→ Contains results (time series) of an emission rate analysis (i.e., including plume velocity data), specific for one PCS line and one velocity retrieval (e.g., the analysis shown in Figure 11 creates three EmissionRates objects for each of the three different velocity retrievals, Section 3.7).
Further important classes (not shown in Figure 2) are: • PlumeBackgroundModel (plumebackground.py)→ performs τ image modelling using either of the available modelling methods (Section 3.3).

•
VeloCrossCorrEngine (plumespeed.py)→ high level class to calculate the plume velocity using the cross-correlation method (Section 3.5.1).Most of the example and introduction scripts provided with Pyplis are based on a short example dataset recorded at Mt. Etna, Italy on 16 September 2015, using a type Envicam-2 camera.It includes ∼15 min of plume data (between 7:06 a.m. and 7:22 a.m.UTC, see e.g., Figure 11.) as well as cell calibration data including suitable background images (between 6:59 a.m. and 7:04 a.m.UTC, see Figure 5).These data are used for demonstration purposes in the provided example scripts, which are summarised in Table A2.The data is not part of the Pyplis installation and can be downloaded from the website.

Figure 1 .
Figure 1.Measurement geometry-sketched (a) in three dimensions and (b) as the camera sees it.The emission plume is indicated in yellow colours, gas velocities in red.The orientation of the plume cross-section (PCS) A (green colours) is aligned with the camera optical axis k.

Figure 2 .
Figure 2. Pyplis API flowchart illustrating the central analysis steps and the corresponding classes of the Pyplis API for SO 2 emission rate retrievals.Italic denotations correspond to class names in Pyplis.Optional/alternative analysis procedures are indicated by dashed boxes.Setup classes (red)include relevant meta information and can be used to create Dataset objects (blue).The latter perform file separation by image type and create ImgList objects (green) for each type (e.g., on, off, dark).Further analysis classes are indicated in yellow.Note that the routine for signal dilution correction is not shown here (cf.FigureA1).

Figure 3 .
Figure 3. Measurement geometry 2D overview map of the measurement setup at Mt. Etna from the Pyplis example data.The map includes plume orientation (dark green line), camera azimuth retrieved using a compass (dashed light green line), the corrected camera azimuth (light green line) and the corresponding FOV (semi-transparent green), retrieved automatically using the pixel-position i, j of the Etna south-east crater within the images and the corresponding coordinates of the crater (longitude, latitude, altitude).The map was generated using an instance of the MeasGeometry class (see Section 3.1).
fitting polynomial of n-th order using sky reference pixels along horizontal profile line (default: n = 2, i.e., quadratic polynomial) 0 no Masked 2D polynomial surface fit 99 yes None (use I 0 -img as is) (a) Method 0 (b) Method 1 (c) Method 4 (d) Method 6

Figure 4 .
Figure 4. Plume background modelling.Four examples of on-band OD images (τ on ) determined using background modelling methods 0 (polynomial surface fit) and 1, 4 and 6 (based on I 0 -image, cf.Table2).Horizontal and vertical profiles are plotted on the top and on the right, respectively.The top-left plot (Method 0) further includes the mask specifying sky-reference pixels (green area) that was used for the polynomial surface fit.The plume and background images used for the displayed examples were recorded at 7:06 a.m.UTC and at 7:02 a.m.UTC, respectively.

Figure 5 .
Figure 5. Cell calibration -(a) shows the output of the automatic cell search routine (included in the CellCalibEngine class); (b) shows example calibration curves for τ on , τ off and τ AA images for the image center pixel from the dataset shown in (a).The error bars in (b) correspond to uncertainties in the cell SO 2 -CDs from the DOAS analysis (y-axis errors) and the x-errorbars correspond to the min / max-range of the corresponding τ-values over the entire image (arising from effects due to non-perpendicular illumination, see Section 2.2 for details).

Figure 6 .
Figure 6.DOAS calibration curves.Plot showing example DOAS calibration data corresponding to the two different FOV parametrisations shown in Figure7.The fit was performed using a first-order weighted polynomial fit.The weights were calculated using the relative errors ∆S SO2 / S SO2 of the DOAS SO 2 -CDs.The y-axis offset is likely due to uncertainties in the DOAS retrieval (e.g., due to O 3 interference) and is compensated by the calibration.

Figure 7 .
Figure 7. DOAS FOV search.DOAS FOV search results using (a), the Pearson and (b), the IFR method.The Pearson method (a) yields an FOV centered at i = 159, j = 124 and a disk radius of four pixels.For the IFR retrieval (b), a tolerance factor of λ = 2 × 10 −3 was chosen and a Super-Gaussian (without tilt) was fitted to the correlation image yielding an FOV centered at i = 159.3,j = 123.8,σ = 7.1, asymmetry parameter a = 1.9 and a shape parameter of b = 0.3 (for details see Equation (A2)).The retrieved FOV positions show good agreement.Note that the FOV was retrieved from downscaled images (Gauss pyramid level 2).

Figure 8 .
Figure 8. Plume velocity retrieval using the cross-correlation method applied to a time-series of on-band OD (τ On ) images.Left: example plume image including the two used PCS lines.Right: Time-series of the integrated optical densities along both lines (orange dashed and cyan line) in addition to the PCS signal shifted using the retrieved correlation lag of 22.1 s (orange profile).The analysis yields an average gas velocity of 3.5 m/s.

Figure 9 .
Figure 9. Plume velocity optical flow.Left: example optical flow vector field (blue lines with red dots depicting the movement direction) calculated using the Farnebäck algorithm.Middle, right: histograms of the flow vector orientation angles ϕ and magnitudes |f|, respectively, corresponding to the ROI shown in the left image (semi-transparent rectangular area around the displayed PCS line).From the latter, a plume velocity of 4.2 (±0.4) m/s and a predominant movement direction of θ = −65 (±14) • was retrieved using first and second moments of the displayed histogram distributions.
. The correction yields an emission rate of 4.8 (±1.2) kg/s using the PCS line shown in Figure10d.The emission rate of the uncorrected image is 1.6 (±0.8) kg/s, corresponding to an underestimation of approximately 67 % at an average plume distance of 10.4 (±0.1) km.

Figure 10 .
Figure 10.Signal dilution correction.(a) and (b) show the fit results of an example dilution analysis for the on and off-band regime, respectively; (c) shows a 3D map of the camera scene and (d) a dilution corrected SO 2 -CD image.The terrain features used for the dilution analysis are indicated in (c) and (d) (blue and lime coloured lines); (d) further includes the image region used to estimate the ambient intensity I A (blue rectangle) and an example PCS line used to compare emission rates before and after the correction.

Figure 11 .
Figure 11.Emission rate retrieval.Emission rates (top) and plume velocities (bottom) of the Etna example dataset on 16 September 2015 between 7:06 a.m. and 7:22 a.m.UTC using the retrieval line "PCS" shown in Figure 8.The analysis was performed using (1) a global velocity of 3.5 m/s (cyan, from cross-correlation analysis); (2) the raw output of the Farnebäck algorithm (thin orange);(3) using the flow_histo method (thin dashed brown) and (4) using the flow_hybrid method (bold orange).The retrieved effective velocities are plotted in the lower panel and correspond to the average velocities along the PCS line using the flow_hybrid method.Uncertainties are displayed as shaded areas.

Table 1 .
Analysis blocks for emission rate retrievals (for details see Section 2).

Table 2 .
Available plume background modelling methods of the PlumeBackgroundModel class (cf.names of sky reference areas in Figure4).

Table A1 .
Performance analysis of the five image preparation steps listed above, dependent on Gaussian pyramid level (pyrlevel) and Gaussian blurring (blur) including the relative percentage impact of the optical flow computation.
• DilutionCorr (dilutioncorr.py)→engine to perform signal dilution correction (Section 3.6).•ImgStack (processing.py)→ contains a series of images (stored as 3D numpy array) as well as supplementary data (e.g., acq.time stamps, exposure times of all stacked images) and basic processing operations (time merging with other data, up / downscaling), can be saved as FITS file.
Appendix F. Supplementary Information and Test DataAppendix F.1.Example Dataset and Example Scripts