Pyplis-A Python Software Toolbox for the Analysis of SO 2 Camera Data

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). In the past years, 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. This inspired the development of Pyplis, an open-source software toolbox written in Python 2.7, which aims to unify the most prevalent 5 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 within the camera images. Plume velocities can be retrieved using an optical flow algorithm as well as 10 signal cross-correlation. Furthermore, Pyplis includes a routine to perform a first order correction of the signal dilution effect. 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. Pyplis has been extensively and successfully tested using data from several field campaigns. Here, 15 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 magma 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, inter alia, 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 DOAS technique (differential optical absorption spectroscopy, [8], e.g.[9], [10]).During the past years, the comparatively young technique of UV SO 2 cameras has gained in importance, since it enables to study volcanic SO 2 emissions at unprecedented spatial and temporal resolution (e.g.[11], [12], [13], [14]).This is particularly helpful to study multiple sources independently (e.g.[15]) or to investigate volcanic degassing characteristics by studying periodicities in the SO 2 -emission-rates (e.g.[16]).The technique of UV SO 2 cameras has seen remarkable improvements in the past years (e.g.[17], [18], [19], [20], [21]) and can nowadays 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 the plume imagery comprises several analysis steps which are summarised in Table 1 and are illustrated in the flowchart shown in the 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 ( [22]), the IDL source code provided by [19] and Plumetrack ( [23]).The first two programs include routines for cell calibration and cross-correlation based plume speed retrievals.The IDL program also includes a routine to perform a first order correction for the signal dilution effect.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 (Python plume imaging software) 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 0.13.4 are (details follow in Sect.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 Sect. 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 Sect.3.

UV SO 2 cameras
UV SO 2 cameras analyse scattered sun light which 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.[17]) 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 Eq. 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.[18]).The latter is more accurate in case aerosols are present in the plume ( [18], [24]).The position and shape of the DOAS FOV within the camera images is required in order to perform the DOAS calibration.The FOV can either be measured experimentally (e.g. in the lab) or can be retrieved directly from the field data ( [18], [25]).The filter transmission curves shift towards lower wavelengths for non-perpendicular illumination.This leads to an increased SO 2 sensitivity towards the image edges.A correction for this effect can be performed during the image calibration using a sensitivity mask retrieved from cell calibration data (e.g.[18]).

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.[26]) or optical flow algorithms (e.g.[20]).Eq. 3 is equivalent to commonly used retrieval methods (e.g.[11], [24]) 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 extend 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.[27]).
If a locally uniform gas velocity can be assumed (i.e.v(i, j) → v) and further, a planar PCS is used for the retrieval (i.e.n(i, j) → n) then, Eq. 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.[28], [29]).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], [30]).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 ( [19]).More complicated radiative transfer issues (e.g.optically thick plumes) require corrections using radiative transfer models (e.g.[24]).Pyplis can correct for the signal dilution effect based on the method suggested by [19].

Implementation
In the following, the main features and modules of Pyplis are presented.The API is designed modular using 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 D).Pyplis includes 19 example scripts, which provide an introduction into the main features of the software (summarised in Table A1).The scripts are based on an example dataset recorded at Mt. Etna in Italy on 16  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.Eq. 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.
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 (shown in Figure 10, for details see Appendix A).Furthermore, 2D and 3D topographic overview maps of the measurement setup can be created automatically (as shown in Figure 3, see also Figure 10c).The MeasGeometry class uses the Python Geonum library ( [31]) which is briefly introduced in Appendix A.

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 B. 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.Eq. 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 1 st or 2 nd order polynomials 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 which are considered during the fit.The mask can either be provided by the user or can be retrieved automatically using the method find_sky_background.The PlumeBackgroundModel class (part of plumebackground.py)provides high level access to eight default methods for the background retrieval.The methods are summarised in Table 2. Figure 4 shows four example plume on-band OD images calculated using the modelling modes 0, 1, 4 and 6.It is not intended to give a recommendation for a "best" method here, as this strongly depends on the data 165 (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 filled with a known amount of SO 2 -gas, or using plume SO 2 -CDs retrieved from a co-located DOAS spectrometer (cf.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 ex.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 1 .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 C.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).7).The fit was performed using a first-order weighted polynomial fit.The weights where 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.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 Sect. 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.Sect.3.4.1)which is then normalised to the centre-position of the DOAS FOV (illustrated in ex.script 7).

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 [18]).Assuming a circular FOV shape, the pixel extend 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. IFR routine: this method is based on [25] 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 C.2 for details).The retrieved FOV information (position, shape, convolution mask) is represented by the DoasFOV class and can be stored as 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 Sect.Peer-reviewed version available at Geosciences 2017, 7, 134; doi:10.3390/geosciences70401343.5.1.Velocity retrieval using the ICA cross-correlation method For the cross-correlation method, ICA time-series (see Eq. 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 2 .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).

Optical flow based velocity retrievals
Optical flow velocity retrievals are performed using the Farnebäck algorithm ( [33]) which is implemented in the OpenCV library ( [34]).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.[35]).
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 include a routine that can detect and correct for unphysical motion estimates in homogeneous image regions based on the method proposed by [35].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.OpenCV documentation).The Pyplis default settings follow the suggestions of [20].An example flow field is shown in Figure 9 including results from the histogram post analysis within a narrow ROI around an example retrieval line.The latter yield 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 ex.script 9 (see Table A1).

Image based signal dilution correction
A correction for the signal dilution effect (see Sect. 2.4) can be performed using the DilutionCorr class.The method is based on [19] and uses the model function to retrieve atmospheric extinction coefficients in the on and off-band regime.Here, 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 ( [19]).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 Eq. 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 Eq. 4 in [19]).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 Sect.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 06:45 UTC.Extinction coefficients of on = 0.0743 km −1 (Figure 10a) and off = 0.0654 km −1 (Figure 10b) were retrieved using the two topographic profile lines shown in Figure 10c and 10d.The correction yields an emission-rate of 4.8 (± 1.2) kg/s using the PCS line shown in Figure 10d.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.

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 D.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 Sect.3.5.1)or can be retrieved automatically during the emission-rate analysis using the Farnebäck optical flow routine (see Sect. 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.Peer-reviewed version available at Geosciences 2017, 7, 134; doi:10.3390/geosciences7040134 2. flow_histo → performs the histogram post-analysis proposed by [35] (cf.Sect.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 prev.point).
Figure 11.Emission-rate retrieval -Emission-rates (top) and plume velocities (bottom) of the Etna example dataset on 16 September 2015 between 07:06 -07:22 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, velo: glob, from the cross-correlation analysis), 2.) the raw output of the Farnebäck algorithm (thin orange) and 3.) 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.
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 mode 6 (cf.Sect.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.

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.[26]) or cell and DOAS calibration (e.g.[18]), Pyplis incorporates more recent developments.These include an implementation of the DOAS FOV retrieval algorithm proposed by [25], a routine to correct for the signal dilution effect based on [19], or pixel based gas velocity retrievals using an optical flow algorithm (e.g.[20]).The latter incorporates a method to detect and correct for ill-constraint optical flow motion-vectors based on [35].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.The Pyplis installation includes numerous example scripts which 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: https://github.com/jgliss/pyplis.The code documentation and further information (e.g.installation instructions) can be found on the documentation website: http://pyplis.readthedocs.io.including suitable background images (between 06:59 -07:04 UTC, see Figure 5).These data is used for demonstration purposes in the provided example scripts, which are summarised in Table A1.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 main 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 analysis steps 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).

PreprintsFigure 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 (red dashed line), camera azimuth retrieved using a compass (purple dashed line) and the corrected camera azimuth (blue dashed line) and corresponding FOV (semi-transparent green), retrieved automatically using the 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 a MeasGeometry object (see Sect.3.1).

Figure 4 .
Figure 4. Plume background modelling -Four example on-band OD images (τ on ) determined using background modelling modes 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 (Mode 0) further includes the mask specifying sky-reference pixels (green area) which was used for the polynomial surface fit.The plume and background images used for the displayed examples were recorded at 07:06 UTC and at 07:02 UTC, respectively.

( a )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).

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 where 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 a FOV centered at i = 159, j = 124 and a disk radius of 4 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 a 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 Eq. A2).The retrieved FOV positions show good agreement.Note, that the FOV was retrieved from downscaled images (Gauss pyramid level 2).

Figure 8
shows results from an example cross-correlation analysis, resulting in a plume velocity of 3.5 m/s.

Figure 8 .
Figure 8. Plume velocity retrieval using cross-correlation method.Left: example on-band OD image and the two used PCS lines.Right: ICA time-series for both lines (orange dashed and cyan line) and 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 -Example output of the Farnebäck optical flow algorithm including histograms of the flow vector orientation angles ϕ (middle) and magnitudes |f| (right) within an example ROI around the indicated PCS line (orange).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.

( a )Figure 10 .
Figure 10.Signal dilution correction -(a) and (b) show the fit result 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.

Table 1 .
Analysis blocks for emission-rate retrievals

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

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 13 October 2017 doi:10.20944/preprints201710.0085.v1
Peer-reviewed version available at Geosciences 2017, 7, 134; doi:10.3390/geosciences7040134Due 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.