Example scripts

pyplis example scripts. The scripts require downloading the pyplis test data.

Note

The scripts are based on the latest commit in the GitHub repo. If you have installed an older version of pyplis, please use the corresponding scripts which are provided here.

Introductory scripts

These scripts give an introduction into basic features and classes of pyplis.

Example 0.1 - Image representation

Introduction into Img object and basic usage including correction for dark current and detector offset.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script 1 - Image import.

This script loads an exemplary image which is part of the pyplis supplied
data. Image data in pyplis is represented by the ``Img`` class, which also
allows for storage of image meta data and keeps track of changes applied to
the image data (e.g. cropping, blurring, dark correction).

The script also illustrates how to manually work with image meta data
encoded in the image file name. The latter can be performed automatically in
pyplis using file name conventions (which can be specified globally, see next
script).
"""
from __future__ import (absolute_import, division)

# Check script version
from SETTINGS import check_version

# imports from SETTINGS.py
from SETTINGS import OPTPARSE, SAVE_DIR

from os.path import join, basename
from datetime import datetime
from matplotlib.pyplot import close
import pyplis

check_version()

# file name of test image stored in data folder
IMG_FILE_NAME = "test_201509160708_F01_335.fts"

IMG_DIR = join(pyplis._LIBDIR, "data")

if __name__ == "__main__":
    close("all")

    img_path = join(IMG_DIR, IMG_FILE_NAME)

    # Create Img object (Img objects can be initiated both with image file
    # paths but also with data in memory in form of a 2D numpy array)
    img = pyplis.image.Img(img_path)

    # log mean of uncropped image for testing mode
    avg = img.mean()

    # The file name of the image includes some image meta information which can
    # be set manually (this is normally done automatically by defining a file
    # name convention, see next script)

    # split filename using delimiter "_"
    spl = IMG_FILE_NAME.split(".")[0].split("_")

    # extract acquisition time and convert to datetime
    acq_time = datetime.strptime(spl[1], "%Y%m%d%H%M")
    # extract and convert exposure time from filename (convert from ms -> s)
    texp = float(spl[3]) / 1000

    img.meta["start_acq"] = acq_time
    img.meta["texp"] = texp
    img.meta["f_num"] = 2.8
    img.meta["focal_length"] = 25e-3

    # add some blurring to the image
    img.add_gaussian_blurring(sigma_final=3)

    # crop the image edges
    roi_crop = [100, 100, 1244, 924]  # [x0, y0, x1, y1]
    img.crop(roi_abs=roi_crop)

    # apply down scaling (gauss pyramid)
    img.to_pyrlevel(2)

    # ## Show image
    img.show()

    img.save_as_fits(SAVE_DIR, "ex0_1_imgsave_test")

    img_reload = pyplis.Img(join(SAVE_DIR, "ex0_1_imgsave_test.fts"))

    # print image information
    print(img)

    # ## IMPORTANT STUFF FINISHED - everything below is of minor importance
    # for educational purposes

    (options, args) = OPTPARSE.parse_args()
    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt

        vals = [img.mean(), avg, int(spl[1]), texp,
                img_reload.meta["texp"],
                img_reload.meta["f_num"],
                img_reload.meta["focal_length"]]

        npt.assert_almost_equal([2526.4623422672885,
                                 2413.0870433989026,
                                 201509160708,
                                 0.335,
                                 0.335,
                                 2.8,
                                 0.025],
                                vals, 4)

        print("All tests passed in script: %s" % basename(__file__))

Example 0.2 - The camera class

Introduction into the Camera object using the example of the ECII camera standard.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Introduction script 2 - The Camera class.

This script introduces the camera class which is of fundamental importance
for image import (e.g. separating on, off, dark background images, etc.) and
also for the data analysis as it includes all relevant detector specifications,
such as number of pixels, pixel size, focal length, etc.

In this script, a newer version of the camera type "ecII" is created manually
in order to illustrate all relevant parameters. The only difference to the
classic ecII camera is, that the filter setup is different.

See also here for more information:

https://pyplis.readthedocs.io/en/latest/tutorials.html#data-import-specifying-
custom-camera-information
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version, OPTPARSE

from os.path import basename
from numpy.testing import assert_array_equal

import pyplis

# Check script version
check_version()

# ## SCRIPT OPTIONS

# Save the new camera as default in database (cam_info.txt file that can be
# found in the "data" directory of the installation.
SAVE_TO_DATABASE = True


# ## SCRIPT FUNCTION DEFINITIONS
def create_ecII_cam_new_filters():
    # Start with creating an empty Camera object
    cam = pyplis.setupclasses.Camera()

    # Specify the camera filter setup

    # Create on and off band filters. Obligatory parameters are "type" and
    # "acronym", "type" specifies the filter type ("on" or
    # "off"), "acronym" specifies how to identify this filter in the file
    # name. If "id" is unspecified it will be equal to "type". The parameter
    # "meas_type_acro" is only important if a measurement type (e.g. M -> meas,
    # C -> calib ...) is explicitely specified in the file name.
    # This is not the case for ECII camera but for the HD camera,
    # see specifications in file cam_info.txt for more info.

    on_band = pyplis.utils.Filter(id="on", type="on", acronym="F01",
                                  meas_type_acro="F01", center_wavelength=310)
    off_band = pyplis.utils.Filter(type="off", acronym="F02",
                                   center_wavelength=330)

    # put the two filter into a list and assign to the camera
    filters = [on_band, off_band]

    cam.default_filters = filters
    cam.prepare_filter_setup()

    # Similar to the filter setup, access info for dark and offset images needs
    # to be specified. The ECII typically records 4 different dark images, two
    # recorded at shortest exposure time -> offset signal predominant, one at
    # low and one at high read gain. The other two recorded at longest possible
    # exposure time -> dark current predominant, also at low and high read gain

    offset_low_gain = pyplis.utils.DarkOffsetInfo(id="offset0", type="offset",
                                                  acronym="D0L", read_gain=0)
    offset_high_gain = pyplis.utils.DarkOffsetInfo(id="offset1", type="offset",
                                                   acronym="D0H", read_gain=1)
    dark_low_gain = pyplis.utils.DarkOffsetInfo(id="dark0", type="dark",
                                                acronym="D1L", read_gain=0)
    dark_high_gain = pyplis.utils.DarkOffsetInfo(id="dark1", type="dark",
                                                 acronym="D1H", read_gain=1)

    # put the 4 dark info objects into a list and assign to the camera
    dark_info = [offset_low_gain, offset_high_gain,
                 dark_low_gain, dark_high_gain]

    cam.dark_info = dark_info

    # Now specify further information about the camera

    # camera ID (needs to be unique, i.e. not included in data base, call
    # pyplis.inout.get_all_valid_cam_ids() to check existing IDs)
    cam.cam_id = "ecII_new_test"

    # image file type
    cam.file_type = "fts"

    # File name delimiter for information extraction
    cam.delim = "_"

    # position of acquisition time (and date) string in file name after
    # splitting with delimiter
    cam.time_info_pos = 3

    # datetime string conversion of acq. time string in file name
    cam.time_info_str = "%Y%m%d%H%M%S%f"

    # position of image filter type acronym in filename
    cam.filter_id_pos = 4

    # position of meas type info
    cam.meas_type_pos = 4

    # Define which dark correction type to use
    # 1: determine a dark image based on image exposure time using a dark img
    # (with long exposure -> dark current predominant) and a dark image with
    # shortest possible exposure (-> detector offset predominant). For more
    # info see function model_dark_image in processing.py module
    # 2: subtraction of a dark image recorded at same exposure time than the
    # actual image
    cam.darkcorr_opt = 1

    # If the file name also includes the exposure time, this can be specified
    # here:
    cam.texp_pos = ""  # the ECII does not...

    # the unit of the exposure time (choose from "s" or "ms")
    cam.texp_unit = ""

    # define the main filter of the camera (this is only important for cameras
    # which include, e.g. several on band filters.). The ID need to be one of
    # the filter IDs specified above
    cam.main_filter_id = "on"

    # camera focal length can be specified here (but does not need to be, in
    # case of the ECII cam, there is no "default" focal length, so this is left
    # empty)
    cam.focal_length = ""

    # Detector geometry
    cam.pix_height = 4.65e-6  # pixel height in m
    cam.pix_width = 4.65e-6  # pixel width in m
    cam.pixnum_x = 1344
    cam.pixnum_y = 1024

    cam._init_access_substring_info()

    cam.io_opts = dict(USE_ALL_FILES=False,
                       SEPARATE_FILTERS=True,
                       INCLUDE_SUB_DIRS=True,
                       LINK_OFF_TO_ON=True)

    # Set the custom image import method
    cam.image_import_method = pyplis.custom_image_import.load_ecII_fits
    # That's it...
    return cam


# ## SCRIPT MAIN FUNCTION
if __name__ == "__main__":

    cam = create_ecII_cam_new_filters()

    print(cam)

    try:
        # you can add the cam to the database (raises error if ID conflict
        # occurs, e.g. if the camera was already added to the database)
        cam.save_as_default()
    except KeyError:
        print("Camera already exists in database")
    cam_reload = pyplis.Camera("ecII_brandnew")
    # ## IMPORTANT STUFF FINISHED - everything below is of minor importance
    # for educational purposes

    (options, args) = OPTPARSE.parse_args()
    # apply some tests. This is done only if TESTMODE is active: testmode can
    # be activated globally (see SETTINGS.py) or can also be activated from
    # the command line when executing the script using the option --test 1
    if int(options.test):
        # quick and dirty test

        cam_dict_nominal = {'darkcorr_opt': 1,
                            '_fid_subnum_max': 1,
                            '_fname_access_flags': {'filter_id': False,
                                                    'meas_type': False,
                                                    'start_acq': False,
                                                    'texp': False},
                            '_mtype_subnum_max': 1,
                            '_time_info_subnum': 1,
                            'cam_id': 'ecII_new_test',
                            'delim': '_',
                            'file_type': 'fts',
                            'filter_id_pos': 4,
                            'focal_length': '',
                            'image_import_method':
                                pyplis.custom_image_import.load_ecII_fits,
                            'main_filter_id': 'on',
                            'meas_type_pos': 4,
                            'pix_height': 4.65e-06,
                            'pix_width': 4.65e-06,
                            'pixnum_x': 1344,
                            'pixnum_y': 1024,
                            'ser_no': 9999,
                            'texp_pos': '',
                            'texp_unit': '',
                            'time_info_pos': 3,
                            'time_info_str': '%Y%m%d%H%M%S%f'}

        from collections import OrderedDict
        geom_data_nominal = OrderedDict([('lon', None),
                                         ('lat', None),
                                         ('altitude', None),
                                         ('azim', None),
                                         ('azim_err', None),
                                         ('elev', None),
                                         ('elev_err', None),
                                         ('alt_offset', 0.0)])

        arr_nominal = list(geom_data_nominal.items())
        arr_nominal.extend(list(cam_dict_nominal.items()))

        arr_vals = list(cam.geom_data.items())
        for k in cam_dict_nominal:
            arr_vals.append((k, cam.__dict__[k]))

        assert_array_equal(arr_nominal, arr_vals)

        print("All tests passed in script: %s" % basename(__file__))

Example 0.3 - Introduction into ImgList objects

Manual creation of ImgList objects and basic features.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script 3 - manual creation of image lists.

This script gives an introduction into the creation and the handling of
ImgList objects. Image lists normally contain images of a certain type (e.g.
on, off, dark, offset). The separation by these file types is normally done
automatically within a Dataset object (see e.g. following script
ex0_4_imglists_auto.py or example script ex01_analysis_setup.py) using a
certain camera type (see previous script ex0_2_camera_setup.py).

In this example these features are NOT used but rather, an on and off band
image list is created manually without previous definition of the Camera type
(i.e. file naming convention, etc.). The script starts with creating an ImgList
containing all images of type FITS found in the test data image base directory.
On and off band lists are then consecutively extracted from this list using
the list method ``separate_by_substr_filename``.

Furthermore, some basic image preparation features of ImgList objects are
introduced (e.g. linking of lists, dark correction, automatic blurring,
cropping, size reduction).
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version, IMG_DIR, SAVEFIGS, SAVE_DIR, FORMAT, DPI,\
    OPTPARSE

import pyplis
from matplotlib.pyplot import subplots, close, show
from datetime import datetime

from os.path import join, isfile, basename
from os import listdir

# Check script version
check_version()

# ## RELEVANT DIRECTORIES AND PATHS
OFFSET_FILE = join(IMG_DIR, "EC2_1106307_1R02_2015091607064723_D0L_Etna.fts")
DARK_FILE = join(IMG_DIR, "EC2_1106307_1R02_2015091607064865_D1L_Etna.fts")

if __name__ == "__main__":
    close("all")

    # ## Get all images in the image path which are FITS files (actually all)
    all_paths = [join(IMG_DIR, f) for f in listdir(IMG_DIR) if
                 isfile(join(IMG_DIR, f)) and f.endswith("fts")]

    # Let pyplis know that this is the ECII camera standard, so that it is
    # able to import meta information from the FITS header
    cam = pyplis.Camera("ecII")  # loads default info for ECII camera

    # ## Now put them all into an image list

    # Note that the files are not separated by filter type, or dark and offset,
    # etc. so the list simply contains all images of type fts which were found
    # in IMG_DIR
    list_all_imgs = pyplis.imagelists.ImgList(all_paths, list_id="all",
                                              camera=cam)

    # Split the list by on band file type (which is identified by acronym
    # "F01" at 4th position in file name after splitting using delimiter "_")
    # creates two new list objects, one containing matches (i.e. on band images
    # the other containing the rest)
    on_list, rest = list_all_imgs.separate_by_substr_filename(sub_str="F01",
                                                              sub_str_pos=4,
                                                              delim="_")
    # now extract all off band images from the "rest" list
    off_list, rest = rest.separate_by_substr_filename(sub_str="F02",
                                                      sub_str_pos=4,
                                                      delim="_")

    # Link the off band list to on band list (the index in the offband list is
    # changed whenever the index is changed in the on band list based on acq.
    # time stamp). Note, that a Camera class was not defined here (see prev.
    # script). Specify the required information in the camera (manually, not
    # required if camera was defined beforehand)
    on_list.camera.delim = "_"
    on_list.camera.time_info_pos = 3
    on_list.camera.time_info_str = "%Y%m%d%H%M%S%f"
    off_list.camera.delim = "_"
    off_list.camera.time_info_pos = 3
    off_list.camera.time_info_str = "%Y%m%d%H%M%S%f"

    on_list.link_imglist(off_list)

    # ## Load dark and offset images and set them in the on-band image list
    dark_img = pyplis.image.Img(DARK_FILE,
                                import_method=cam.image_import_method)
    offset_img = pyplis.image.Img(OFFSET_FILE,
                                  import_method=cam.image_import_method)

    on_list.add_master_dark_image(dark_img)
    on_list.add_master_offset_image(offset_img)

    # Go to image number 100 in on-band list
    on_list.goto_img(100)

    # ## Change image preparation settings (these are applied on image load)

    # region of interest
    on_list.roi_abs = [100, 100, 1300, 900]

    # activate cropping in image list
    on_list.crop = 1

    # scale down to pyramid level 2
    on_list.pyrlevel = 2

    # activate automatic dark correction in list
    on_list.DARK_CORR_OPT = 1  # see previous script
    on_list.darkcorr_mode = True

    # blur the on band images
    on_list.add_gaussian_blurring(3)

    # get current on band image (no. 100 in list), all previously set
    # preparation settings are applied to this image
    on_img = on_list.current_img()

    # get current off band image (no. 100 in list, index was automatically
    # changed since it is linked to on band list. Note that this image is
    # unedited
    off_img = off_list.current_img()

    fig, ax = subplots(1, 2, figsize=(18, 6))

    on_img.show(ax=ax[0])
    on_time_str = datetime.strftime(on_img.start_acq, "%Y-%m-%d %H:%M:%S")
    ax[0].set_title("Current img (on-band list): %s" % on_time_str)

    # show the current off band image as well (this image is unedited)
    off_img.show(ax=ax[1])
    off_time_str = datetime.strftime(off_img.start_acq, "%Y-%m-%d %H:%M:%S")
    ax[1].set_title("Current img (off-band list): %s" % off_time_str)

    on_list.edit_info()

    on_list.skip_files = 3
    on_list.goto_img(15)

    # ## IMPORTANT STUFF FINISHED
    if SAVEFIGS:
        fig.savefig(join(SAVE_DIR, "ex0_3_out_1.%s" % FORMAT),
                    format=FORMAT, dpi=DPI)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt

        npt.assert_array_equal([501, 1, 500, 0],
                               [on_list.nof + off_list.nof,
                                on_list.this.edit_log["darkcorr"],
                                sum(on_list.this.shape),
                                on_list.gaussian_blurring -
                                on_list.this.edit_log["blurring"]])

        npt.assert_allclose([402.66284],
                            [off_img.mean() - on_img.mean()],
                            rtol=1e-7, atol=0)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except Exception:
        print("Use option --show 1 if you want the plots to be displayed")

Example 0.4 - Introduction into the Dataset class

Automatic image type separation using the Dataset object and the ECII camera standard.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script 4 - Automatic creation of image lists.

The previous script gave an introduction into the manual creation of
``ImgList`` objects and some basic image preparation features.
In this script, a number of ImgList objects (on, off, dark low gain,
dark high gain, offset low gain, offset high gain) is created automatically
using the Camera class created in example script ex0_2_camera_setup.py (ECII
camera)

Based on the information stored in the Camera class, a MeasSetup class is
created. The latter class collects all meta information relevant for an
emission rate analysis. Apart from the camera specs, this may include source
definitions  contains information about the camera specs a the
image base directory (note that in this example, start / stop acq. time stamps
are ignored, i.e. all images available in the specified directory are imported)
"""
from __future__ import (absolute_import, division)

# Imports from SETTINGS.py
from SETTINGS import check_version, IMG_DIR, OPTPARSE

import pyplis
from os.path import basename

# ## IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex0_2_camera_setup import create_ecII_cam_new_filters

# Check script version
check_version()


if __name__ == "__main__":
    # create the camera which was
    cam = create_ecII_cam_new_filters()

    # now throw all this stuff into the BaseSetup objec
    stp = pyplis.setupclasses.MeasSetup(IMG_DIR, camera=cam)

    # Create a Dataset which creates separate ImgLists for all types (dark,
    # offset, etc.)
    ds = pyplis.dataset.Dataset(stp)

    # The image lists can be accessed in different ways for instance using
    # the method "all_lists", which returns a Python list containing all
    # ImgList objects that were created within the Dataset
    all_imglists = ds.all_lists()

    # print some information about each of the lists
    for lst in all_imglists:
        print("list_id: %s, list_type: %s, number_of_files: %s"
              % (lst.list_id, lst.list_type, lst.nof))

    # single lists can be accessed using "get_list(<id>)" using a valid ID,
    # e.g.:
    on_list = ds.get_list("on")
    off_list = ds.get_list("off")

    on_list.goto_img(50)  # this also changes the index in the off band list...

    # ... because it is linked to the on band list (automatically set in
    # Dataset)
    print("\nImgLists linked to ImgList on: %s" % on_list.linked_lists.keys())
    print("Current file number on / off list: %d / %d\n" % (on_list.cfn,
                                                            off_list.cfn))

    # Detected dark and offset image lists are also automatically linked to the
    # on and off band image list, such that dark image correction can be
    # applied
    on_list.darkcorr_mode = True
    off_list.darkcorr_mode = True

    # the current image preparation settings can be accessed via the
    # edit_info method
    on_list.edit_info()

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt

        npt.assert_array_equal([501, 2, 2368, 0, 50],
                               [on_list.nof + off_list.nof,
                                on_list.this.is_darkcorr +
                                off_list.this.is_darkcorr,
                                sum(on_list.this.shape),
                                on_list.gaussian_blurring -
                                on_list.this.edit_log["blurring"],
                                on_list.cfn])

        npt.assert_allclose(actual=[on_list.get_dark_image().mean()],
                            desired=[190.56119],
                            rtol=1e-7)

        print("All tests passed in script: %s" % basename(__file__))

Example 0.5 - Optical flow live view

Live view of optical flow calculation using OpticalFlowFarneback object (requires webcam).

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script 5 - Optical flow Farneback liveview using webcam.

Create an OpticalFlowFarneback object and activate live view (requires webcam)
"""
from __future__ import (absolute_import, division)

# Check script version
from SETTINGS import check_version

import pyplis

check_version()

flow = pyplis.plumespeed.OptflowFarneback()
flow.live_example()

Example 0.6 - Plume cross section lines

Creation and orientation of LineOnImage objects for emission rate retrievals.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script no 6 - LineOnImage objects and their orientation.

This script introduces how to create LineOnImage objects within the image plane
and specify their orientation (i.e. the direction into which the normal vector
of the line points). This is mainly important for emission rate retrievals,
where, for instance velcoity displacement vectors (e.g. from an optical flow
algorithm) have to be multiplied with the normal vector of such a line (using
the dot product).
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version, SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE

from pyplis import LineOnImage
from os.path import join, basename
from matplotlib.pyplot import show, subplots, close
from matplotlib.cm import get_cmap

# Check script version
check_version()


def create_example_lines():
    """Create some exemplary lines."""
    lines_r = []  # lines with orientation "right" are stored within this list
    lines_l = []  # lines with orientation "left" are stored within this list

    cmap = get_cmap("jet")
    # horizontal line, normal orientation to the top (0 deg)
    lines_r.append(LineOnImage(x0=10, y0=10, x1=90, y1=10,
                               normal_orientation="right",
                               color=cmap(0),
                               line_id="Line 1"))

    # horizontal line, normal orientation to the bottom (180 deg)
    lines_l.append(LineOnImage(10, 10, 90, 10, normal_orientation="left",
                               color=cmap(0),
                               line_id="Line 1"))

    # Vertical line normal to the right (90 deg)
    lines_r.append(LineOnImage(10, 10, 10, 90, normal_orientation="right",
                               color=cmap(50),
                               line_id="Line 2"))

    # Vertical line normal to the left (270 deg)
    lines_l.append(LineOnImage(10, 10, 10, 90, normal_orientation="left",
                               color=cmap(50),
                               line_id="Line 2"))

    # Slanted line 45 degrees
    lines_r.append(LineOnImage(20, 30, 50, 60, normal_orientation="right",
                               color=cmap(100),
                               line_id="Line 3"))

    # Slanted line 45 degrees
    lines_l.append(LineOnImage(20, 30, 50, 60, normal_orientation="left",
                               color=cmap(100),
                               line_id="Line 3"))

    # Slanted line 45 degrees
    lines_r.append(LineOnImage(90, 10, 10, 90, normal_orientation="right",
                               color=cmap(150),
                               line_id="Line 4"))

    # Slanted line 45 degrees
    lines_l.append(LineOnImage(90, 10, 10, 90, normal_orientation="left",
                               color=cmap(150),
                               line_id="Line 4"))

    lines_r.append(LineOnImage(60, 20, 80, 90, normal_orientation="right",
                               color=cmap(200),
                               line_id="Line 5"))

    lines_l.append(LineOnImage(60, 20, 80, 90, normal_orientation="left",
                               color=cmap(200),
                               line_id="Line 5"))

    lines_r.append(LineOnImage(40, 20, 30, 90, normal_orientation="right",
                               color=cmap(250),
                               line_id="Line 6"))
    lines_l.append(LineOnImage(40, 20, 30, 90, normal_orientation="left",
                               color=cmap(250),
                               line_id="Line 6"))

    return lines_r, lines_l


if __name__ == "__main__":
    close("all")
    fig, ax = subplots(1, 2, figsize=(18, 9))

    lines_r, lines_l = create_example_lines()

    for k in range(len(lines_r)):
        line = lines_r[k]
        # print "%d: %s" %(k, line.orientation_info)
        normal = line.normal_vector
        lbl = "%s" % line.line_id
        line.plot_line_on_grid(ax=ax[0], include_normal=1,
                               include_roi_rot=True, label=lbl)
    for k in range(len(lines_l)):
        line = lines_l[k]
        normal = line.normal_vector
        lbl = "%s" % line.line_id
        line.plot_line_on_grid(ax=ax[1], include_normal=1,
                               include_roi_rot=True, label=lbl)

    ax[0].set_title("Orientation right")
    ax[0].legend(loc="best", fontsize=8, framealpha=0.5)
    ax[0].set_xlim([0, 100])
    ax[0].set_ylim([100, 0])

    ax[1].set_title("Orientation left")
    ax[1].legend(loc="best", fontsize=8, framealpha=0.5)
    ax[1].set_xlim([0, 100])
    ax[1].set_ylim([100, 0])
    # ## IMPORTANT STUFF FINISHED
    if SAVEFIGS:
        fig.savefig(join(SAVE_DIR, "ex0_6_out_1.%s" % FORMAT),
                    format=FORMAT, dpi=DPI)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt

        npt.assert_array_equal([sum([x.length() for x in lines_l]),
                                sum([x.length() for x in lines_r]),
                                sum([sum(x.coords) for x in lines_r]),
                                int(sum([sum(x.normal_vector + y.normal_vector)
                                    for x, y in zip(lines_l, lines_r)])),
                                int(sum([sum(x.center_pix) for x in lines_l])),
                                ],
                               [459, 459, 1030, 0, 515])

        npt.assert_allclose(
            actual=[sum([sum(x) for x in lines_l[2].rect_roi_rot])],
            desired=[368.79393923],
            rtol=1e-7)

        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except Exception:
        print("Use option --show 1 if you want the plots to be displayed")

Example 0.7 - Manual cell calibration

Manual cell calibration based on a set of background and cell images (on / off).

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script no. 5: manual cell calibration.

Perform manual cell calibration based on 3 cell images and one background
image. The calibration data consists of 3 cells which were put in front of the
lense successively and a background images both before and after the
cell images.

This script creates an empty CellCalibEngine object in which 6 cell images
are assigned (3 on and 3 off band) with their corresponding SO2 column
densities. Further, 2 background images are assigned for each filter, one
before and one after the cells were put in front of the camera. These are used
to determine tau images from the cell images. Variations in the background
intensities are corrected for (for details see manuscript).

Note, that this calibration does not include a dark correction of the images
before the calibration, therefore, the results are slightly different compared
to the results from ex05_cell_calib_auto.py.
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

import pyplis
from os.path import join, basename
from matplotlib.pyplot import close, show
from time import time

# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE, IMG_DIR

# Check script version
check_version()

# SPECIFY IMAGE PATHS FOR EACH CELL AND BACKGROUND IMAGES
BG_BEFORE_ON = "EC2_1106307_1R02_2015091607000845_F01_Etna.fts"
BG_BEFORE_OFF = "EC2_1106307_1R02_2015091607001021_F02_Etna.fts"

A53_ON = "EC2_1106307_1R02_2015091607003032_F01_Etna.fts"
A53_OFF = "EC2_1106307_1R02_2015091607003216_F02_Etna.fts"

A37_ON = "EC2_1106307_1R02_2015091607005847_F01_Etna.fts"
A37_OFF = "EC2_1106307_1R02_2015091607010023_F02_Etna.fts"

A57_ON = "EC2_1106307_1R02_2015091607013835_F01_Etna.fts"
A57_OFF = "EC2_1106307_1R02_2015091607014019_F02_Etna.fts"

BG_AFTER_ON = "EC2_1106307_1R02_2015091607020256_F01_Etna.fts"
BG_AFTER_OFF = "EC2_1106307_1R02_2015091607020440_F02_Etna.fts"

# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    close("all")
    start = time()
    # define camera for ECII custom image import
    cam = pyplis.Camera("ecII")
    cellcalib = pyplis.cellcalib.CellCalibEngine(cam)

    # now add all cell images manually, specifying paths, the SO2 column
    # densities of each cell, and the corresponding cell ID as well as image
    # type (on, off)
    cellcalib.set_cell_images(img_paths=join(IMG_DIR, A53_ON),
                              cell_gas_cd=4.15e17,
                              cell_id="a53", filter_id="on")

    cellcalib.set_cell_images(img_paths=join(IMG_DIR, A53_OFF),
                              cell_gas_cd=4.15e17,
                              cell_id="a53", filter_id="off")

    cellcalib.set_cell_images(img_paths=join(IMG_DIR, A37_ON),
                              cell_gas_cd=8.59e17,
                              cell_id="a37", filter_id="on")

    cellcalib.set_cell_images(img_paths=join(IMG_DIR, A37_OFF),
                              cell_gas_cd=8.59e17,
                              cell_id="a37", filter_id="off")

    cellcalib.set_cell_images(img_paths=join(IMG_DIR, A57_ON),
                              cell_gas_cd=1.92e18,
                              cell_id="a57", filter_id="on")

    cellcalib.set_cell_images(img_paths=join(IMG_DIR, A57_OFF),
                              cell_gas_cd=1.92e18,
                              cell_id="a57", filter_id="off")

    # put the onband background images into a Python list ....
    bg_on_paths = [join(IMG_DIR, BG_BEFORE_ON), join(IMG_DIR, BG_AFTER_ON)]

    # ... and add them to the calibration object ...
    cellcalib.set_bg_images(img_paths=bg_on_paths, filter_id="on")

    # ... same for off band background images
    bg_off_paths = [join(IMG_DIR, BG_BEFORE_OFF), join(IMG_DIR, BG_AFTER_OFF)]
    cellcalib.set_bg_images(img_paths=bg_off_paths, filter_id="off")

    # Prepare calibration data (i.e. CellCalibData objets) for on, off and
    # AA images. This function determines tau images for each cell using
    # background images scaled to the present background intensities at the
    # acq. time stamp of each cell using temporal interpolation of the provided
    # background images. This results in 3 tau images for each filter (on, off)
    # and from that, 3 AA images are determined. Each of these collection of
    # 3 tau images (on, off, AA) is then stored within a CellCalibData object
    # which can be accessed using e.g. cellcalib.calib_data["aa"]
    cellcalib.prepare_calib_data(on_id="on", off_id="off", darkcorr=False)
    stop = time()

    ax = cellcalib.plot_all_calib_curves()
    ax.set_title("Manual cell calibration\n(NO dark correction performed)")

    aa_calib = cellcalib.calib_data["aa"]

    print("Time elapsed for preparing calibration data: %.4f s"
          % (stop - start))
    # ## IMPORTANT STUFF FINISHED
    if SAVEFIGS:
        ax.figure.savefig(join(SAVE_DIR, "ex0_7_out_1.%s" % FORMAT),
                          format=FORMAT, dpi=DPI)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt

        slope, offs = aa_calib.calib_coeffs
        actual = [cellcalib.tau_stacks["aa"].sum(),
                  cellcalib.tau_stacks["aa"].mean(),
                  aa_calib.cd_vec.sum(),
                  slope,
                  offs]

        npt.assert_allclose(actual=actual,
                            desired=[1007480.56,
                                     0.24401481,
                                     3.194000052267778e+18,
                                     4.779782684462987e+18,
                                     -2.7244424951657216e+16],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Examples for emission rate analysis

The following scripts are directly related for emission rate analysis and build on top of each other ending with a full emission rate analysis in Example 12 - Emission rate analysis (Etna example data).

Example 1 - Creation of analysis setup and Dataset

This script introduces the MeasSetup object and how it can be used to specify all relevant information to create a Dataset for emission rate analysis.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 1 - Analysis setup for example data set.

In this script an example data set, recorded on the 16/9/15 7:06-7:22 at
Mt. Etna is setup. Most of the following example scripts will use the
information specified here.

A typical analysis setup for plume image data contains information about the
camera (e.g. optics, file naming convention, see also ex0_2_camera_setup.py),
gas source, the measurement geometry and the wind conditions. These information
generally needs to be specified by the user before the analysis. The
information is stored within a MeasSetup object which can be used as a basis
for further analysis.
If not all neccessary information is provided, a MeasSetup object will be
created nonetheless but analysis options might be limited.

Such a MeasSetup object can be used as input for Dataset objects which
creates the analysis environment. The main purpose of Dataset classes
is to automatically separate images by their type (e.g. on / off-band
images, dark, offset) and create ImgList classes from that. ImgList
classes typically contain images of one type. The Dataset also links
relevant ImgList to each other, e.g. if the camera has an off-band
filter, and off-band images can be found in the specified directory,
then it is linked to the list containing on-band images. This means,
that the current image index in the off-band list is automatically
updated whenever the index is changed in the on-band list. If
acquisition time is available in the images, then the index is updated
based on the closest acq. time of the off-band images, based on the
current on-band acq. time. If not, it is updated by index (e.g. on-band
contains 100 images and off-band list 50. The current image number in
the on-band list is 50, then the off-band index is set to 25). Also,
dark and offset images lists are linked both to on and off-band image
lists, such that the image can be correccted for dark and offset
automatically on image load (the latter needs to be activated in the
lists using ``darkcorr_mode=True``).

This script shows how to setup a MeasSetup object and create a Dataset object
from it. As example, the first image of the on-band image time series is
displayed.

The Dataset object created here is used in script  ex04_prep_aa_imglist.py
which shows how to create an image list displaying AA images.
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version, IMG_DIR, OPTPARSE
import pyplis as pyplis
from datetime import datetime
from matplotlib.pyplot import show, close

# Check script version
check_version()

# SCRIPT FUNCTION DEFINITIONS


def create_dataset():
    """Initialize measurement setup and creates dataset from that."""
    start = datetime(2015, 9, 16, 7, 6, 00)
    stop = datetime(2015, 9, 16, 7, 22, 00)
    # Define camera (here the default ecII type is used)
    cam_id = "ecII"

    # the camera filter setup
    filters = [pyplis.utils.Filter(type="on", acronym="F01"),
               pyplis.utils.Filter(type="off", acronym="F02")]

    # camera location and viewing direction (altitude will be retrieved
    # automatically)
    geom_cam = {"lon": 15.1129,
                "lat": 37.73122,
                "elev": 20.0,
                "elev_err": 5.0,
                "azim": 270.0,
                "azim_err": 10.0,
                "alt_offset": 15.0,
                "focal_length": 25e-3}  # altitude offset (above topography)

    # Create camera setup
    # the camera setup includes information about the filename convention in
    # order to identify different image types (e.g. on band, off band, dark..)
    # it furthermore includes information about the detector specifics (e.g.
    # dimension, pixel size, focal length). Measurement specific parameters
    # (e.g. lon, lat, elev, azim) where defined in the dictinary above and
    # can be passed as additional keyword dictionary using **geom_cam
    # Alternatively, they could also be passed directly, e.g.:

    # cam = pyplis.setup.Camera(cam_id, filter_list=filters, lon=15.1129,
    #                           lat=37.73122)

    cam = pyplis.setupclasses.Camera(cam_id, filter_list=filters,
                                     **geom_cam)

    # Load default information for Etna. This information is stored in
    # the source_info.txt file of the Pyplis information. You may also access
    # information about any volcano via the available online access to the NOAA
    # database using the method pyplis.inout.get_source_info_online(source_id).

    source = pyplis.setupclasses.Source("etna")

    # Provide wind direction
    wind_info = {"dir": 0.0,
                 "dir_err": 1.0}

#                "dir_err"  : 15.0}

    # Create BaseSetup object (which creates the MeasGeometry object)
    stp = pyplis.setupclasses.MeasSetup(IMG_DIR, start, stop, camera=cam,
                                        source=source,
                                        wind_info=wind_info)
    print(stp.LINK_OFF_TO_ON)
    # Create analysis object (from BaseSetup)
    # The dataset takes care of finding all vali
    return pyplis.Dataset(stp)


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    close("all")
    ds = create_dataset()

    # get on-band image list
    on_list = ds.get_list("on")
    on_list.goto_next()
    off_list = ds.get_list("off")

    # activate dark correction in both lists. Dark and offset image lists are
    # automatically assigned to plume on and off-band image lists on initiation
    # of the dataset object
    on_list.darkcorr_mode = True
    off_list.darkcorr_mode = True

    print("On-band list contains %d images, current image index: %d"
          % (on_list.nof, on_list.cfn))

    img = on_list.current_img()

    # plume distance image retrieved from MeasGeometry class...
    plume_dists = on_list.plume_dists

    # ...these may be overwritten or set manually if desired
    on_list.plume_dists = 10000

    # The same applies for the integration step lengths for emission rate
    # retrievals
    step_lengths = on_list.integration_step_length
    on_list.integration_step_length = 1.8  # m

    img_shift = img.duplicate()

    # images can be shifted using the scipy.ndimage.interpolation.shift method
    # this may be required for image registration in dual camera systems.
    # Whether this is supposed to be done automatically can be specified using
    # the REG_SHIFT_OFF option in a MeasSetup class. It may also be specified
    # directly for your cam in the custom camera definition file cam_info.txt
    # using io_opts:REG_SHIFT_OFF=1 (see e.g. defintion of camera with ID
    # "usgs"). Also, a default registration offset can be defined here using
    #
    img_shift.shift(dx_abs=-30, dy_abs=55)
    img_shift.show(tit="Shifted")
    # Set pixel intensities below 2000 to 0 (method of Img class)
    img.set_val_below_thresh(val=0, threshold=2000)
    # show modified image
    img.show()
    print(str(img))  # image object has an informative string representation

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename

        actual = [plume_dists.mean(), plume_dists.std(),
                  on_list.get_dark_image().mean()]
        npt.assert_allclose(actual=actual,
                            desired=[10909.873427010458, 221.48844132471388,
                                     190.56119],
                            rtol=1e-7)

        npt.assert_array_equal([418, 2, 2368, 1, 1, 0,
                                20150916070600,
                                20150916072200],
                               [on_list.nof + off_list.nof,
                                on_list.this.is_darkcorr +
                                off_list.this.is_darkcorr,
                                sum(on_list.this.shape),
                                on_list.cfn,
                                off_list.cfn,
                                sum(img.img[img.img < 2000]),
                                int(ds.setup.start.strftime("%Y%m%d%H%M%S")),
                                int(ds.setup.stop.strftime("%Y%m%d%H%M%S"))])

        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 2 - Measurement Geometry

This script introduces the MeasGeometry object including some basic features.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 2 - Features of the MeasGeometry class.

In this script, some important features of the MeasGeometry class are
introduced. The class itself is automatically created in the MeasSetup
object which was created in example script ex01_analysis_setup.py and was
passed as input for a Dataset object. The relevant MeasGeometry class is stored
within the Dataset object and can be accessed via the ``meas_geometry``
attribute.

As a first feature, the viewing direction of the camera is retrieved from the
image using the position of the south east (SE) crater of Mt.Etna. The result
is then visualized in a 2D map to give an overview of the geometry. The map
further includes the initial viewing direction (see example script
ex01_analysis_setup.py) which was logged in the field using a compass and a
mechanical inclinometer.

Further, the distance to the plume is retrieved on a pixel basis (represented
as image).
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

from geonum import GeoPoint
from matplotlib.pyplot import subplots, show, close
from os.path import join

# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE

# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex01_analysis_setup import create_dataset

# Check script version
check_version()

# SCRIPT FUNCTION DEFINITIONS


def find_viewing_direction(meas_geometry, draw_result=True):
    """Correct viewing direction using location of Etna SE crater.

    Defines location of Etna SE crater within images (is plotted into current
    plume onband image of dataset) and uses its geo location to retrieve the
    camera viewing direction

    :param meas_geometry: :class:`MeasGeometry` object

    """
    # Position of SE crater in the image (x, y)
    se_crater_img_pos = [806, 736]

    # Geographic position of SE crater (extracted from Google Earth)
    # The GeoPoint object (geonum library) automatically retrieves the altitude
    # using SRTM data
    se_crater = GeoPoint(37.747757, 15.002643, name="SE crater")

    print("Retrieved altitude SE crater (SRTM): %s" % se_crater.altitude)

    # The following method finds the camera viewing direction based on the
    # position of the south east crater.
    new_elev, new_azim, _, basemap =\
        meas_geometry.find_viewing_direction(pix_x=se_crater_img_pos[0],
                                             pix_y=se_crater_img_pos[1],
                                             # for uncertainty estimate
                                             pix_pos_err=100,
                                             geo_point=se_crater,
                                             draw_result=draw_result,
                                             update=True)  # overwrite settings

    print("Updated camera azimuth and elevation in MeasGeometry, new values: "
          "elev = %.1f, azim = %.1f" % (new_elev, new_azim))

    return meas_geometry, basemap


def plot_plume_distance_image(meas_geometry):
    """Determine and plot image plume distance and pix-to-pix distance imgs."""
    # This function returns three images, the first corresponding to pix-to-pix
    # distances in horizontal direction and the second (ignored here) to
    # the vertical (in this case they are the same since pixel height and
    # pixel width are the same for this detector). The third image gives
    # camera to plume distances on a pixel basis
    (dist_img, _, plume_dist_img) =\
        meas_geometry.compute_all_integration_step_lengths()

    fig, ax = subplots(2, 1, figsize=(7, 8))

    # Show pix-to-pix distance image
    dist_img.show(cmap="gray", ax=ax[0], zlabel="Pix-to-pix distance [m]")
    ax[0].set_title("Parameterised pix-to-pix dists")

    # Show plume distance image (convert pixel values to from m -> km)
    (plume_dist_img / 1000.0).show(cmap="gray",
                                   ax=ax[1], zlabel="Plume distance [km]")
    ax[1].set_title("Plume dists")
    return fig


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    close("all")

    # Create the Dataset object (see ex01)
    ds = create_dataset()

    # execute function defined above (see above for definition and information)
    geom_corr, map_ = find_viewing_direction(ds.meas_geometry)

    # execute 2. script function (see above for definition and information)
    fig = plot_plume_distance_image(ds.meas_geometry)

    # You can compute the plume distance for the camera CFOV pixel column just
    # by calling the method plume_dist() without specifying the input azimuth
    # angle...
    plume_dist_cfov = geom_corr.plume_dist()[0][0]

    # ... and the corresponding uncertainty
    plume_dist_err_cfov = geom_corr.plume_dist_err()

    # You can also retrieve an array containing the camera azimuth angles for
    # each pixel column...
    all_azimuths = geom_corr.all_azimuths_camfov()

    # ... and use this to compute plume distances on a pixel column basis
    plume_dists_all_cols = geom_corr.plume_dist(all_azimuths)

    # If you want, you can update information about camera, source or
    # meteorolgy using either of the following methods (below we apply a
    # change in the wind-direction such that the plume propagation direction
    # is changed from S to SE )
    # geom_corr.update_cam_specs()
    # geom_corr.update_source_specs()

    geom_corr.update_wind_specs(dir=315)
    geom_corr.draw_map_2d()  # this figure is only displayed and not saved

    # recompute plume distance of CFOV pixel
    plume_dist_cfov_new = geom_corr.plume_dist()[0][0]

    print("Comparison of plume distances after change of wind direction:\n"
          "Previous: %.3f m\n"
          "New: %.3f m" % (plume_dist_cfov, plume_dist_cfov_new))
    # Using this a
    # IMPORTANT STUFF FINISHED
    if SAVEFIGS:
        map_.ax.figure.savefig(join(SAVE_DIR, "ex02_out_1.%s" % FORMAT),
                               format=FORMAT, dpi=DPI)
        fig.savefig(join(SAVE_DIR, "ex02_out_2.%s" % FORMAT), format=FORMAT,
                    dpi=DPI)

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename

        npt.assert_array_equal([],
                               [])

        # check some propoerties of the basemap (displayed in figure)

        # map basemap coordinates to lon / lat values
        lon, lat = map_(8335, 9392, inverse=True)
        npt.assert_allclose(actual=[lon, lat, map_.delta_lon, map_.delta_lat],
                            desired=[15.075131135, 37.76678834,
                                     0.149263852982,
                                     0.118462659944],
                            rtol=1e-7)

        # check some basic properties / values of the geometry
        npt.assert_allclose(actual=[geom_corr.cam_elev,
                                    geom_corr.cam_elev_err,
                                    geom_corr.cam_azim,
                                    geom_corr.cam_azim_err,
                                    plume_dist_cfov,
                                    plume_dist_err_cfov,
                                    plume_dists_all_cols.mean(),
                                    plume_dist_cfov_new],
                            desired=[1.547754e+01,
                                     1.064556e+00,
                                     2.793013e+02,
                                     1.065411e+00,
                                     1.073102e+04,
                                     1.645586e+02,
                                     1.076047e+04,
                                     9.961593e+03],
                            rtol=1e-6)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 3 - Plume background analysis

Introduction into PlumeBackgroundModel object the default modes for retrieval of plume background intensities and plume optical density images.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 3 - Plume background analysis.

This example script introduces features related to plume background modelling
and tau image calculations.
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

import numpy as np
from os.path import join
import pyplis
from matplotlib.pyplot import show, subplots, close

# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, OPTPARSE

# Check script version
check_version()

# SCRIPT OPTIONS

# If this is True, then sky reference areas are set in auto mode (note that
# in this case, the tests at the end of the script will fail!)
USE_AUTO_SETTINGS = False

# intensity threshold to init mask for bg surface fit
POLYFIT_2D_MASK_THRESH = 2600

# Choose the background correction modes you want to use

BG_CORR_MODES = [0,  # 2D poly surface fit (without sky radiance image)
                 1,  # Scaling of sky radiance image
                 4,  # Scaling + linear gradient correction in x & y direction
                 6]  # Scaling + quadr. gradient correction in x & y direction


# Image file paths relevant for this script
PLUME_FILE = join(IMG_DIR, 'EC2_1106307_1R02_2015091607065477_F01_Etna.fts')
BG_FILE = join(IMG_DIR, 'EC2_1106307_1R02_2015091607022602_F01_Etna.fts')
OFFSET_FILE = join(IMG_DIR, 'EC2_1106307_1R02_2015091607064723_D0L_Etna.fts')
DARK_FILE = join(IMG_DIR, 'EC2_1106307_1R02_2015091607064865_D1L_Etna.fts')

# SCRIPT FUNCTION DEFINITIONS


def init_background_model():
    """Create background model and define relevant sky reference areas."""
    # Create background modelling object
    m = pyplis.plumebackground.PlumeBackgroundModel()

    # Define default gas free areas in plume image
    w, h = 40, 40  # width/height of rectangles

    m.scale_rect = [1280, 20, 1280 + w, 20 + h]
    m.xgrad_rect = [20, 20, 20 + w, 20 + h]
    m.ygrad_rect = [1280, 660, 1280 + w, 660 + h]

    # Define coordinates of horizontal and vertical profile lines

    # row number of profile line for horizontal corretions in the sky
    # gradient...
    m.xgrad_line_rownum = 40
    # ... and start / stop columns for the corrections
    m.xgrad_line_startcol = 20
    m.xgrad_line_stopcol = 1323

    # col number of profile line for vertical corretions in the sky gradient...
    m.ygrad_line_colnum = 1300
    # ... and start / stop rows for the corrections
    m.ygrad_line_startrow = 10
    m.ygrad_line_stoprow = 700
    # Order of polyonmial fit applied for the gradient correction
    m.ygrad_line_polyorder = 2

    return m


def load_and_prepare_images():
    """Load images defined above and prepare them for the background analysis.

    Returns
    -------
        - Img, plume image
        - Img, plume image vignetting corrected
        - Img, sky radiance image

    """
    # get custom load method for ECII
    fun = pyplis.custom_image_import.load_ecII_fits
    # Load the image objects and peform dark correction
    plume, bg = pyplis.Img(PLUME_FILE, fun), pyplis.Img(BG_FILE, fun)
    dark, offset = pyplis.Img(DARK_FILE, fun), pyplis.Img(OFFSET_FILE, fun)

    # Model dark image for tExp of plume image
    dark_plume = pyplis.image.model_dark_image(plume.meta["texp"],
                                               dark, offset)
    # Model dark image for tExp of background image
    dark_bg = pyplis.image.model_dark_image(bg.meta["texp"],
                                            dark, offset)

    plume.subtract_dark_image(dark_plume)
    bg.subtract_dark_image(dark_bg)
    # Blur the images (sigma = 1)
    plume.add_gaussian_blurring(1)
    bg.add_gaussian_blurring(1)

    # Create vignetting correction mask from background image
    vign = bg.img / bg.img.max()  # NOTE: potentially includes y & x gradients
    plume_vigncorr = pyplis.Img(plume.img / vign)
    return plume, plume_vigncorr, bg


def autosettings_vs_manual_settings(bg_model):
    """Perform automatic retrieval of sky reference areas.

    If you are lazy... (i.e. you dont want to define all these reference areas)
    then you could also use the auto search function, a comparison is plotted
    here.
    """
    auto_params = pyplis.plumebackground.find_sky_reference_areas(plume)
    current_params = bg_model.settings_dict()

    fig, axes = subplots(1, 2, figsize=(16, 6))
    axes[0].set_title("Manually set parameters")
    pyplis.plumebackground.plot_sky_reference_areas(plume, current_params,
                                                    ax=axes[0])
    pyplis.plumebackground.plot_sky_reference_areas(plume, auto_params,
                                                    ax=axes[1])
    axes[1].set_title("Automatically set parameters")

    return auto_params, fig


def plot_pcs_profiles_4_tau_images(tau0, tau1, tau2, tau3, pcs_line):
    """Plot PCS profiles for all 4 methods."""
    fig, ax = subplots(1, 1)
    tau_imgs = [tau0, tau1, tau2, tau3]

    for k in range(4):
        img = tau_imgs[k]
        profile = pcs_line.get_line_profile(img)
        ax.plot(profile, "-", label=r"Mode %d: $\phi=%.3f$"
                % (BG_CORR_MODES[k], np.mean(profile)))

    ax.grid()
    ax.set_ylabel(r"$\tau_{on}$", fontsize=20)
    ax.set_xlim([0, pcs_line.length()])
    ax.set_xticklabels([])
    ax.set_xlabel("PCS", fontsize=16)
    ax.legend(loc="best", fancybox=True, framealpha=0.5, fontsize=12)
    return fig


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    close("all")

    # Create a background model with relevant sky reference areas
    bg_model = init_background_model()

    # Define exemplary plume cross section line
    pcs_line = pyplis.LineOnImage(x0=530,
                                  y0=730,
                                  x1=890,
                                  y1=300,
                                  line_id="example PCS",
                                  color="lime")

    plume, plume_vigncorr, bg = load_and_prepare_images()

    auto_params, fig0 = autosettings_vs_manual_settings(bg_model)

    # Script option
    if USE_AUTO_SETTINGS:
        bg_model.update(**auto_params)

    # Model 4 exemplary tau images

    # list to store figures of tau plotted tau images
    _tau_figs = []

    # mask for corr mode 0 (i.e. 2D polyfit)
    mask = np.ones(plume_vigncorr.img.shape, dtype=np.float32)
    mask[plume_vigncorr.img < POLYFIT_2D_MASK_THRESH] = 0

    # First method: retrieve tau image using poly surface fit
    tau0 = bg_model.get_tau_image(plume_vigncorr,
                                  mode=BG_CORR_MODES[0],
                                  surface_fit_mask=mask,
                                  surface_fit_polyorder=1)

    # Plot the result and append the figure to _tau_figs
    _tau_figs.append(bg_model.plot_tau_result(tau0, PCS=pcs_line))

    # Second method: scale background image to plume image in "scale" rect
    tau1 = bg_model.get_tau_image(plume, bg, mode=BG_CORR_MODES[1])
    _tau_figs.append(bg_model.plot_tau_result(tau1, PCS=pcs_line))

    # Third method: Linear correction for radiance differences based on two
    # rectangles (scale, ygrad)
    tau2 = bg_model.get_tau_image(plume, bg, mode=BG_CORR_MODES[2])
    _tau_figs.append(bg_model.plot_tau_result(tau2, PCS=pcs_line))

    # 4th method: 2nd order polynomial fit along vertical profile line
    # For this method, determine tau on tau off and AA image
    tau3 = bg_model.get_tau_image(plume, bg, mode=BG_CORR_MODES[3])
    _tau_figs.append(bg_model.plot_tau_result(tau3, PCS=pcs_line))

    fig6 = plot_pcs_profiles_4_tau_images(tau0, tau1, tau2, tau3, pcs_line)

    if SAVEFIGS:
        fig0.savefig(join(SAVE_DIR, "ex03_out_1.%s" % FORMAT), format=FORMAT,
                     dpi=DPI, transparent=True)
        for k in range(len(_tau_figs)):
            # _tau_figs[k].suptitle("")
            _tau_figs[k].savefig(join(SAVE_DIR, "ex03_out_%d.%s"
                                      % ((k + 2), FORMAT)),
                                 format=FORMAT, dpi=DPI)

        fig6.savefig(join(SAVE_DIR, "ex03_out_6.%s" % FORMAT), format=FORMAT,
                     dpi=DPI)
    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename
        m = bg_model

        # test settings for clear sky reference areas
        npt.assert_array_equal([2680, 3960, 160, 6, 1300, 10, 700, 40, 20,
                                1323, 567584],
                               [sum(m.scale_rect),
                                sum(m.ygrad_rect),
                                sum(m.xgrad_rect),
                                m.mode,
                                m.ygrad_line_colnum,
                                m.ygrad_line_startrow,
                                m.ygrad_line_stoprow,
                                m.xgrad_line_rownum,
                                m.xgrad_line_startcol,
                                m.xgrad_line_stopcol,
                                int(m.surface_fit_mask.sum())])

        m.update(**auto_params)
        # test settings for clear sky reference areas
        npt.assert_array_equal([2682, 4142, 1380, 6, 1337, 1, 790, 6, 672,
                                1343, 567584],
                               [sum(m.scale_rect),
                                sum(m.ygrad_rect),
                                sum(m.xgrad_rect),
                                m.mode,
                                m.ygrad_line_colnum,
                                m.ygrad_line_startrow,
                                m.ygrad_line_stoprow,
                                m.xgrad_line_rownum,
                                m.xgrad_line_startcol,
                                m.xgrad_line_stopcol,
                                int(m.surface_fit_mask.sum())])

        # test all tau-modelling results
        actual = [tau0.mean(), tau1.mean(), tau2.mean(), tau3.mean()]
        npt.assert_allclose(actual=actual,
                            desired=[0.11395558008662043,
                                     0.25279653,
                                     0.13842879832119934,
                                     0.13943940574220634],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 4 - Preparation of AA image list

Script showing how to prepare an ImgList containing on-band plume images, such that aa_mode can be activated (i.e. images are loaded as AA images).

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 4 - Prepare AA image list (from onband list).

Script showing how to work in AA mode using ImgList object
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

import pyplis
from matplotlib.pyplot import close
from os.path import join

# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, OPTPARSE, PCS1

# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex01_analysis_setup import create_dataset
from ex02_meas_geometry import find_viewing_direction

# Check script version
check_version()


# SCRIPT FUNCTION DEFINITIONS
def prepare_aa_image_list(bg_corr_mode=6):
    """Get and prepare onband list for aa image mode.

    The relevant gas free areas for background image modelling are set
    automatically (see also ex. 3 for details)

    :return: - on list in AA mode
    """
    dataset = create_dataset()
    geom, _ = find_viewing_direction(dataset.meas_geometry, False)

    # Set plume background images for on and off
    # this is the same image which is also used for example script
    # ex03_plume_background.py demonstrating the plume background routines
    path_bg_on = join(IMG_DIR,
                      'EC2_1106307_1R02_2015091607022602_F01_Etna.fts')
    path_bg_off = join(IMG_DIR,
                       'EC2_1106307_1R02_2015091607022820_F02_Etna.fts')

    # Get on and off lists and activate dark correction
    on_lst = dataset.get_list("on")
    off_list = dataset.get_list("off")

    # Deactivate automatic reload in list while changing some list
    # attributes
    on_lst.auto_reload = False
    off_list.auto_reload = False

    on_lst.darkcorr_mode = True
    off_list.darkcorr_mode = True

    # Prepare on and offband background images
    bg_on = pyplis.Img(path_bg_on)
    bg_on.subtract_dark_image(on_lst.get_dark_image())

    bg_off = pyplis.Img(path_bg_off)
    bg_off.subtract_dark_image(off_list.get_dark_image())

    # set the background images within the lists
    on_lst.set_bg_img(bg_on)
    off_list.set_bg_img(bg_off)

    # automatically set gas free areas
    on_lst.bg_model.set_missing_ref_areas(on_lst.current_img())
    # Now update some of the information from the automatically set sky ref
    # areas
    on_lst.bg_model.xgrad_line_startcol = 20
    on_lst.bg_model.xgrad_line_rownum = 25
    off_list.bg_model.xgrad_line_startcol = 20
    off_list.bg_model.xgrad_line_rownum = 25

    # set background modelling mode
    on_lst.bg_model.mode = bg_corr_mode
    off_list.bg_model.mode = bg_corr_mode

    on_lst.aa_mode = True  # activate AA mode

    off_list.auto_reload = True
    on_lst.auto_reload = True
    print("INITIATED AA LIST")
    on_lst.meas_geometry = geom
    return on_lst


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    from matplotlib.pyplot import show
    from time import time

    close("all")
    aa_list = prepare_aa_image_list()

    t0 = time()
    # Deactivate auto reload while changing some settings (else, after each
    # of the following operations the images are reloaded and edited, which)
    aa_list.auto_reload = False
    aa_list.goto_img(50)

    aa_list.add_gaussian_blurring(1)
    # aa_list.pyrlevel = 2
    aa_list.roi_abs = [300, 300, 1120, 1000]
    aa_list.crop = True
    # now reactive image reload in list (loads image no. 50 with all changes
    # that where set in the previous lines)
    aa_list.auto_reload = True

    # store some results for tests below
    shape_log, mean_log = sum(aa_list.this.shape), aa_list.this.mean()

    ax = aa_list.show_current(zlabel=r"$\tau_{AA}$")
    print("Elapsed time: %s s" % (time() - t0))

    aa_list.crop = False
    ax1 = aa_list.bg_model.plot_sky_reference_areas(aa_list.current_img())
    fig = aa_list.bg_model.plot_tau_result(aa_list.current_img())

    # import plume cross section and computed integrated optical density
    # for current image
    img = aa_list.this
    ax2 = img.show(vmin=-0.1, vmax=0.3)
    pcs = PCS1.convert(to_pyrlevel=0)
    pcs.plot_line_on_grid(ax=ax2)
    pix_steps = aa_list.meas_geometry.compute_all_integration_step_lengths()[0]
    integrated_aa = pcs.integrate_profile(img, pix_steps)

    # IMPORTANT STUFF FINISHED
    if SAVEFIGS:
        ax.figure.savefig(join(SAVE_DIR, "ex04_out_1.%s" % FORMAT),
                          format=FORMAT, dpi=DPI)

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename

        m = aa_list.bg_model
        npt.assert_array_equal([2682, 4144, 1380, 6, 1337, 1, 791, 25, 20,
                                1343],
                               [sum(m.scale_rect),
                                sum(m.ygrad_rect),
                                sum(m.xgrad_rect),
                                m.mode,
                                m.ygrad_line_colnum,
                                m.ygrad_line_startrow,
                                m.ygrad_line_stoprow,
                                m.xgrad_line_rownum,
                                m.xgrad_line_startcol,
                                m.xgrad_line_stopcol])

        actual = [aa_list.meas_geometry.cam_elev,
                  aa_list.meas_geometry.cam_azim,
                  aa_list.meas_geometry.plume_dist()[0, 0],
                  aa_list.this.mean(),
                  shape_log, mean_log]

        npt.assert_allclose(actual=actual,
                            desired=[15.477542212645357,
                                     279.30130009369515,
                                     10731.024327931776,
                                     0.009083584068527644,
                                     1520,
                                     0.014380159209694215],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 5 - Automatic cell calibration

This scripts shows how to perform automatic cell calibration based on a time series of on and off band images containing both, suitable background images and images from different SO2 cells (in both wavelength channels, cf. Example 0.7 - Manual cell calibration).

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 5 - Automatic cell calibration.

Script showing how to use the automatic cell calibration engine which only
requires to specify start / stop time stamps of a calibration window. Based on
that sub time windows for each cell as well as suitable background images are
detected and separated into individual image lists (for all filters, i.e. here
on / off).

"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

import pyplis
from datetime import datetime
from time import time
from os.path import join
from matplotlib.pyplot import show, close

# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, OPTPARSE

# Check script version
check_version()

# File path for storing cell AA calibration data including sensitivity
# correction mask
AA_CALIB_FILE = join(SAVE_DIR, "ex05_cellcalib_aa.fts")


# SCRIPT FUNCTION DEFINITIONS
def perform_auto_cell_calib():
    # Calibration time stamps
    start = datetime(2015, 9, 16, 6, 59, 00)
    stop = datetime(2015, 9, 16, 7, 3, 00)

    # Gas CDs in cells and cell ids
    # See supplementary package data about DOAS fit retrieval
    calib_cells = {'a37': [8.59e17, 2.00e17],
                   'a53': [4.15e17, 1.00e17],
                   'a57': [19.24e17, 3.00e17]}

    # the camera used
    cam_id = "ecII"

    # The camera filter setup is different from the ECII default setup and is
    # therefore defined explicitely
    filters = [pyplis.utils.Filter(type="on", acronym="F01"),
               pyplis.utils.Filter(type="off", acronym="F02")]

    # create camera setup, this includes the filename convention for
    # image separation
    cam = pyplis.setupclasses.Camera(cam_id=cam_id, filter_list=filters)

    # Create CellCalibSetup class for initiation of CellCalib object
    setup = pyplis.setupclasses.MeasSetup(IMG_DIR, start, stop,
                                          camera=cam,
                                          cell_info_dict=calib_cells)

    # Create CellCalibEngine object, read on...
    # This is a DataSet object and performs file separation and creation of
    # on / off, dark / offset lists for all images in the specified time window
    c = pyplis.cellcalib.CellCalibEngine(setup)

    # the following high level method calls several functions in the
    # CellCalibEngine class, most importantly the method find_cells for on and
    # off band image time series, which detects sub time windows for each cell
    # and background images. After the individual time windows are detected for
    # each cell and filter, the method _assign_calib_specs is called, which
    # assigns the SO2 CD amount (specified above in dictionary calib_cells)
    # to the detected sub time windows (both for on and off) based on the depth
    # of the intensity dip (in the onband) for each sub time window (should
    # become clear from the plot produced in this script). Then it creates
    # CellImgList objects for each of the cells and for the detected background
    # images (i.e. resulting in (M + 1) x 2 lists, with M being the number of
    # detected intensity dips, the + 1 is the corresponding background list and
    # times 2 for on / off)
    c.find_and_assign_cells_all_filter_lists()

    return c


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    close("all")
    start = time()
    c = perform_auto_cell_calib()

    # prepare CellCalibData objects for on, off and aa images. These can
    # be accessed via c.calib_data[key] where key is "on", "off", "aa"
    c.prepare_calib_data(pos_x_abs=None,  # change for a specific pix
                         pos_y_abs=None,  # change for a specific pix
                         radius_abs=1,  # radius of retrieval disk
                         on_id="on",  # ImgList ID of onband filter
                         off_id="off")  # ImgList ID of offband filter
    stop = time()
    # Plot search result of on
    ax0 = c.plot_cell_search_result("on", include_tit=False)
    ax1 = c.plot_cell_search_result("off", include_tit=False)
    # Plot all calibration curves for center pixel and in a radial
    # neighbourhood of 20 pixels
    ax2 = c.plot_all_calib_curves()
    ax2.set_xlim([0, 0.7])
    ax2.set_ylim([0, 2.25e18])
    # IMPORTANT STUFF FINISHED
    if SAVEFIGS:
        ax0.figure.savefig(join(SAVE_DIR, "ex05_2_out_1.%s" % FORMAT),
                           format=FORMAT, dpi=DPI)
        ax1.figure.savefig(join(SAVE_DIR, "ex05_2_out_2.%s" % FORMAT),
                           format=FORMAT, dpi=DPI)
        ax2.figure.savefig(join(SAVE_DIR, "ex05_2_out_3.%s" % FORMAT),
                           format=FORMAT, dpi=DPI)

    ax0.set_title("Cell search result on band", fontsize=18)
    ax1.set_title("Cell search result off band", fontsize=18)
    ax2.set_title("Calibration polynomials", fontsize=18)

    # You can explicitely access the individual CellCalibData objects
    aa_calib = c.calib_data["aa"]

    aa_calib.fit_calib_data()
    c.plot_calib_curve("on")
    mask = c.get_sensitivity_corr_mask("aa")

    aa_calib.save_as_fits(AA_CALIB_FILE)
    print("Time elapsed for preparing calibration data: %.4f s"
          % (stop - start))

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename
        calib_reload = pyplis.CellCalibData()
        calib_reload.load_from_fits(AA_CALIB_FILE)
        calib_reload.fit_calib_data(polyorder=2,
                                    through_origin=True)
        # test some basic features of calibraiton dataset (e.g. different
        # ImgList classes for on and off and the different cells)
        npt.assert_array_equal([c.cell_search_performed,
                                c.cell_lists["on"]["a37"].nof,
                                c.cell_lists["on"]["a53"].nof,
                                c.cell_lists["on"]["a57"].nof,
                                c.cell_lists["off"]["a37"].nof,
                                c.cell_lists["off"]["a53"].nof,
                                c.cell_lists["off"]["a57"].nof,
                                calib_reload.calib_id,
                                calib_reload.pos_x_abs,
                                calib_reload.pos_y_abs],
                               [True, 2, 3, 3, 2, 3, 3, "aa", 672, 512])
        d = c._cell_info_auto_search

        vals_approx = [d["a37"][0],
                       d["a53"][0],
                       d["a57"][0],
                       aa_calib.calib_fun(0, *aa_calib.calib_coeffs),
                       calib_reload.calib_fun(0, *calib_reload.calib_coeffs)]
        npt.assert_allclose(actual=vals_approx,
                            desired=[8.59e+17,
                                     4.15e+17,
                                     1.924e+18,
                                     -1.8313164653590516e+16,
                                     0.0],
                            rtol=1e-7)

        # explicitely check calibration data for on, off and aa (plotted in
        # this script)
        actual = [c.calib_data["on"].calib_coeffs.mean(),
                  c.calib_data["off"].calib_coeffs.mean(),
                  aa_calib.calib_coeffs.mean(),
                  calib_reload.calib_coeffs.mean()]
        npt.assert_allclose(actual=actual,
                            desired=[1.8926803829028974e+18,
                                     -3.742539080945949e+19,
                                     2.153653759747737e+18,
                                     2.1197681750384312e+18],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 6 - DOAS calibration

Introduction into DOAS calibration including FOV search using both, the Pearson and the IFR method.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 6 - DOAS calibration and FOV search.

Script showing how to work with DOAS calibration data

In this script, a stack of plume AA images from the Etna test data is imported
as well as a time series of DOAS SO2-CDs in the plume retrieved using the
software DOASIS (see directory "spectra" in test data folder for corresponding
analysis details, the folder also contains the RAW data and the jscript code
for analysing the spectra). The DOAS result import is performed using the
Python package ``pydoas``.

Based on these data, position and shape of the DOAS FOV within the camera
uimages is identified using both FOV search methods (IFR and Pearson). The
results of the FOV search are plotted as well as the corresponding calibration
curves retrieved for both FOV parametrisations.

Note
------

In case a MemoryError occurs while determining the AA image stack, then the
stack (3D numpy array) is too large for the RAM. In this case, try
increasing script option PYRLEVEL_ROUGH_SEARCH.

"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

import pyplis
import pydoas
import numpy.testing as npt
import numpy as np
from datetime import timedelta
from matplotlib.pyplot import close, show, subplots
from os.path import join, exists
from os import remove

# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, OPTPARSE

# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list

# Check script version
check_version()

# SCRIPT OPTONS

# reload and save stack in folder SAVE_DIR, results in increased
# running time due to stack calculation (is automatically switched on if
# the stack is not found at this location)
RELOAD_STACK = False

PYRLEVEL_ROUGH_SEARCH = 2
# Default search settings are at pyramid level 2, the FOV results are upscaled
# to original resolution, if the following option is set 1, then, based on
# the result from pyrlevel=2, another stack is determined at pyrlevel = 0
# (i.e. in full resolution) within a ROI around the center position from the
# search performed at pyrlevel=PYRLEVEL_ROUGH_SEARCH (defined 2 lines below)
DO_FINE_SEARCH = False


# RELEVANT DIRECTORIES AND PATHS

# Directory containing DOAS result files
DOAS_DATA_DIR = join(IMG_DIR, "..", "spectra", "plume_prep", "min10Scans",
                     "ResultFiles")

STACK_PATH = join(SAVE_DIR, "ex06_aa_imgstack.fts")

# SCRIPT FUNCTION DEFINITIONS


def load_doas_results(lt_to_utc_shift=timedelta(-1. / 12)):
    """Specify DOAS data import from DOASIS fit result files.

    In order to perform the DOAS FOV search, as many spectrum datapoints
    as possible are needed. Therefore, only 10 spectra were added (to reduce
    noise) per plume spectrum. The DOAS fit was performed in a wavelength
    range between 314 - 326 nm (fit ID: "f01").
    """
    # This dictionary specifies which information is supposed to be imported
    # from the DOAS fit result files stored in DOAS_DATA_DIR. In the example
    # shown here, only the SO2 fit results are imported from fit scenario
    # with ID "f01" (key of dict). The corresponding value of each key is
    # a list of format ["header_id", ["fit_id1", "fit_id2", ..., "fit_idN"]]
    # specifying the identification string of the species in the result file
    # headers and the second entry is a list specifying all fit scenario IDs
    # from which this species is supposed to be imported (here only f01)
    fit_import_info = {"so2": ["SO2_Hermans_298_air_conv_satCorr1e18",
                               ["f01"]
                               ]}

    # Create a result import setup for the DOAS data based on the import
    # dictionary and the image base directory of the result files ...
    doas_import_setup =\
        pydoas.dataimport.ResultImportSetup(DOAS_DATA_DIR,
                                            result_import_dict=fit_import_info)

    # ... and create a result dataset from that
    doas_dataset = pydoas.analysis.DatasetDoasResults(doas_import_setup)

    # get the SO2 fit results from the dataset. Individual results of certain
    # species can be accessed using the species ID (key in ``fit_import_info``
    # dict) and its fit ID (one of the fit IDs specified for this species, here
    # f01).
    # Note, that the DOAS data was stored using local time, thus they need to
    # be shifted (2h back) to match the camera data time stamps (which are in
    # UTC), otherwise the temporal merging of the two datasets (for the DOAS
    # calibration) does not work
    results_utc = doas_dataset.get_results("so2", "f01").shift(lt_to_utc_shift)
    return results_utc


def make_aa_stack_from_list(aa_list, roi_abs=None, pyrlevel=None,
                            save=True, stack_path=STACK_PATH,
                            save_dir=SAVE_DIR):
    """Get and prepare onband list for aa image mode."""
    # Deactivate auto reload to change some settings (if auto_reload is active
    # list images are reloaded whenever a setting is changed in the list. This
    # can slow down things, thus, if you intend to change a couple of settings
    # you might deactivate auto_reload, adapt the settings and then re-activate
    # auto_reload
    aa_list.auto_reload = False
    if roi_abs is not None:
        aa_list.roi_abs = roi_abs
        aa_list.crop = True
    aa_list.pyrlevel = pyrlevel
    aa_list.auto_reload = True

    # Stack all images in image list at pyrlevel 2 and cropped using specified
    # roi (uncropped if roi_abs=None).
    stack = aa_list.make_stack()
    if save:
        try:
            remove(stack_path)
        except BaseException:
            pass
        stack.save_as_fits(save_dir=save_dir,
                           save_name="ex06_aa_imgstack.fts")
    return stack


def get_stack(reload_stack=RELOAD_STACK, stack_path=STACK_PATH,
              pyrlevel=PYRLEVEL_ROUGH_SEARCH):
    """Load stack data based on current settings."""
    if not exists(stack_path):
        reload_stack = 1

    if not reload_stack:
        stack = pyplis.processing.ImgStack()
        stack.load_stack_fits(stack_path)
        if stack.pyrlevel != pyrlevel:
            reload_stack = True
    aa_list = None
    if reload_stack:
        # import AA image list
        aa_list = prepare_aa_image_list()
        # Try creating stack
        stack = make_aa_stack_from_list(aa_list, pyrlevel=pyrlevel)

    return stack, aa_list


# Test functions used at the end of the script
def test_calib_pears_init(calib):
    calib.fit_calib_data(polyorder=1)
    cc = pyplis.helpers.get_img_maximum(calib.fov.corr_img.img)
    assert cc == (124, 159), cc

    pyrl = calib.fov.pyrlevel
    assert pyrl == 2, pyrl
    res_dict = calib.fov.result_pearson
    npt.assert_allclose([res_dict['rad_rel'],
                         np.max(100 * res_dict['corr_curve'].values)],
                        [3, 95], atol=1)

    fov_ext = calib.fov.pixel_extend(abs_coords=True)
    (fov_x, fov_y) = calib.fov.pixel_position_center(abs_coords=True)

    npt.assert_allclose([fov_ext, fov_x, fov_y],
                        [res_dict['rad_rel'] * 2**pyrl, 636, 496], atol=1)

    npt.assert_allclose(calib.calib_coeffs,
                        [8.58e+18, 2.71e+17], rtol=1e-1)


def test_calib_pears_fine(calib):
    cc = pyplis.helpers.get_img_maximum(calib.fov.corr_img.img)
    npt.assert_allclose((186, 180), cc, atol=1)

    pyrl = calib.fov.pyrlevel
    assert pyrl == 0, pyrl
    res_dict = calib.fov.result_pearson
    npt.assert_allclose([res_dict['rad_rel'],
                         np.max(100 * res_dict['corr_curve'].values)],
                        [6, 95], atol=1)

    fov_ext = calib.fov.pixel_extend(abs_coords=True)
    (fov_x, fov_y) = calib.fov.pixel_position_center(abs_coords=True)

    npt.assert_allclose([fov_ext, fov_x, fov_y],
                        [6, 630, 493], atol=1)

    npt.assert_allclose(calib.calib_coeffs,
                        [8.38e+18, 2.92e+17], rtol=1e-1)


def test_calib_ifr(calib):
    cc = pyplis.helpers.get_img_maximum(calib.fov.corr_img.img)
    npt.assert_allclose((123, 157), cc, atol=1)

    pyrl = calib.fov.pyrlevel
    assert pyrl == 2, pyrl
    npt.assert_allclose(calib.fov.result_ifr['popt'][1:5],
                        [158.6, 122.9, 15.4, 1.5], rtol=1e-1)

    (fov_x, fov_y) = calib.fov.pixel_position_center(abs_coords=True)

    npt.assert_allclose([fov_x, fov_y, calib_ifr.fov.sigma_x_abs,
                         calib_ifr.fov.sigma_y_abs],
                        [635, 492, 61.5, 41.6], atol=2)

    npt.assert_allclose(calib.calib_coeffs,
                        [9.38e+18, 1.75e+17], rtol=1e-1)


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    # close all plots
    close("all")

    # import DOAS results
    doas_time_series = load_doas_results()

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    if options.test:
        # if test mode is active, the image stack is always recomputed from
        # scratch and the option DO_FINE_SEARCH is activated, since the tests
        # are based on this. This will lead to an increased computation time
        # Test mode can be activated / deactivated in SETTINGS.py or via
        # the script option --test 1 (on) or --test 0 (off)
        RELOAD_STACK = True
        PYRLEVEL_ROUGH_SEARCH = 2
        DO_FINE_SEARCH = True

    # reload or create the AA image stack based on current script settings
    stack, aa_list = get_stack()

    s = pyplis.doascalib.DoasFOVEngine(stack, doas_time_series)
    calib_pears = s.perform_fov_search(method="pearson")
    calib_ifr = s.perform_fov_search(method="ifr", ifrlbda=4e-3)

    # plot the FOV search results
    ax0 = calib_pears.fov.plot()
    ax1 = calib_ifr.fov.plot()

    calib_pears.fit_calib_data()
    calib_ifr.fit_calib_data()

    fig, ax2 = subplots(1, 1)
    calib_pears.plot(add_label_str="Pearson", color="b", ax=ax2)

    calib_ifr.plot(add_label_str="IFR", color="g", ax=ax2)
    ax2.set_title("Calibration curves Pearson vs. IFR method")
    ax2.grid()
    ax2.set_ylim([0, 1.8e18])
    ax2.set_xlim([0, 0.20])
    ax2.legend(loc=4, fancybox=True, framealpha=0.7, fontsize=11)
    axes = [ax0, ax1, ax2]

    if DO_FINE_SEARCH:
        """Perform FOV search within ROI around result from pearson fov
        search at full resolution (pyrlevel=0)
        """
        if aa_list is None:
            aa_list = prepare_aa_image_list()
        # remember some properties of the current image stack that is stored in
        # the DoasFOVEngine object (i.e. the merged one in low pixel res.)

        num_merge, h, w = s.img_stack.shape
        s_fine = s.run_fov_fine_search(aa_list, doas_time_series,
                                       method="pearson")

        calib_pears_fine = s_fine.calib_data
        calib_pears_fine.plot()
        calib_pears_fine.fov.plot()

    calib_pears.save_as_fits(save_dir=SAVE_DIR,
                             save_name="ex06_doascalib_aa.fts")
    calib_ifr.save_as_fits(save_dir=SAVE_DIR,
                           save_name="ex06_doascalib_aa_ifr_method.fts")

    # you can also change the order of the calibration polynomial and
    # force it to go through the origin
    calib_pears.fit_calib_data(polyorder=2, through_origin=True)
    calib_pears.plot_calib_fun(add_label_str="Pearson (WRONG,\n"
                               "2nd order, through origin)",
                               color="r", ax=ax2)
    ax2.legend(loc=4, fancybox=True, framealpha=0.7, fontsize=11)
    # IMPORTANT STUFF FINISHED
    if SAVEFIGS:
        for k in range(len(axes)):
            ax = axes[k]
            ax.set_title("")
            ax.figure.savefig(join(SAVE_DIR, "ex06_out_%d.%s"
                                   % ((k + 1), FORMAT)),
                              format=FORMAT, dpi=DPI)

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        from os.path import basename

        num, h, w = stack.shape
        num2 = s.img_stack.shape[0]  # stack after fine FOV search
        prep = stack.img_prep

        # check some basic properties of the data used for the different FOV
        # searches
        npt.assert_array_equal(
            [len(doas_time_series), num, num_merge, h, w, stack.pyrlevel,
             prep["darkcorr"] * prep["is_tau"] * prep["is_aa"], num2],
            [120, 209, 88, 256, 336, 2, 1, 209])

        test_calib_pears_init(calib_pears)
        test_calib_ifr(calib_ifr)
        test_calib_pears_fine(calib_pears_fine)

        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 7 - AA sensitivity correction masks

Combine the results from Example 5 - Automatic cell calibration and Example 6 - DOAS calibration in order to retrieve AA sensitivity correction masks normalised to the position of the DOAS FOV.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 7 - AA sensitivity correction masks.

In this script, cell and DOAS calibration (see previous 2 scripts) of the Etna
test dataset are opposed. Furthermore, it is illustrated, how to create
correction masks for pixel variations in the SO2 sensitivity due to shifts in
the filter transmission windows.

The cell calibration is re-performed (using method ``perform_auto_cell_calib``)
from example script 5. The results from the DOAS calibration
(see prev. example) were stored as a FITS file (including FOV information)
and the results are imported here.
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

import pyplis
from os.path import join, exists
import numpy as np
from matplotlib.pyplot import close, subplots, show
from matplotlib.patches import Circle
import six

# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE

# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex05_cell_calib_auto import perform_auto_cell_calib
from ex04_prep_aa_imglist import prepare_aa_image_list

# Check script version
check_version()

CELL_AA_CALIB_FILE = join(SAVE_DIR, "ex05_cellcalib_aa.fts")
# RELEVANT DIRECTORIES AND PATHS

# fits file containing DOAS calibration information (from ex6)
DOAS_CALIB_FILE = join(SAVE_DIR, "ex06_doascalib_aa.fts")


# SCRIPT FUNCTION DEFINITIONS
def draw_doas_fov(fov_x, fov_y, fov_extend, ax):
    # add FOV position to plot of examplary AA image
    c = Circle((fov_x, fov_y), fov_extend, ec="k", fc="lime", alpha=.5)
    ax.add_artist(c)
    ax.text(fov_x, (fov_y - fov_extend * 1.3), "DOAS FOV")
    ax.set_xlim([0, 1343]), ax.set_ylim([1023, 0])
    return ax


def plot_pcs_comparison(aa_init, aa_imgs_corr, pcs1, pcs2):
    fig, axes = subplots(1, 2, figsize=(18, 6))
    p10 = pcs1.get_line_profile(aa_init.img)
    p20 = pcs2.get_line_profile(aa_init.img)

    num = len(p10)

    axes[0].set_title("Line %s" % pcs1.line_id)
    axes[1].set_title("Line %s" % pcs2.line_id)

    axes[0].plot(p10, "-", label=r"Init $\phi=%.3f$" % (sum(p10) / num))
    axes[1].plot(p20, "-", label=r"Init $\phi=%.3f$" % (sum(p20) / num))

    for cd, aa_corr in six.iteritems(aa_imgs_corr):
        p1 = pcs1.get_line_profile(aa_corr.img)
        p2 = pcs2.get_line_profile(aa_corr.img)

        axes[0].plot(p1, "-", label=r"Cell CD: %.2e $\phi=%.3f$"
                     % (cd, sum(p1) / num))
        axes[1].plot(p2, "-", label=r"Cell CD: %.2e $\phi=%.3f$"
                     % (cd, sum(p2) / num))

    axes[0].legend(loc='best', fancybox=True, framealpha=0.5, fontsize=10)
    axes[1].legend(loc='best', fancybox=True, framealpha=0.5, fontsize=10)
    axes[0].grid()
    axes[1].grid()
    return fig, axes


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    close("all")

    if not exists(DOAS_CALIB_FILE):
        raise IOError("Calibration file could not be found at specified "
                      "location:\n %s\nYou might need to run example 6 first")

    # Load AA list
    aa_list = prepare_aa_image_list()
    aa_list.add_gaussian_blurring(2)

    # Load DOAS calbration data and FOV information (see example 6)
    doascalib = pyplis.doascalib.DoasCalibData()
    doascalib.load_from_fits(file_path=DOAS_CALIB_FILE)
    doascalib.fit_calib_data()

    # Get DOAS FOV parameters in absolute coordinates
    fov_x, fov_y = doascalib.fov.pixel_position_center(abs_coords=True)
    fov_extend = doascalib.fov.pixel_extend(abs_coords=True)

    # Load cell calibration (see example 5)
    cellcalib = perform_auto_cell_calib()

    # get cell calibration
    cellcalib.prepare_calib_data(
        pos_x_abs=fov_x,  # change if you want it for a specific pix
        pos_y_abs=fov_y,  # change if you want it for a specific pix
        radius_abs=fov_extend,  # radius of retrieval disk
        on_id="on",  # ImgList ID of onband filter
        off_id="off")  # ImgList ID of offband filter

    cell_aa_calib = cellcalib.calib_data["aa"]

    # Define lines on image for plume profiles
    pcs1 = pyplis.LineOnImage(620, 700, 940, 280,
                              line_id="center")
    pcs2 = pyplis.LineOnImage(40, 40, 40, 600,
                              line_id="edge")

    # Plot DOAS calibration polynomial
    ax0 = doascalib.plot(add_label_str="DOAS")
    ax0 = cellcalib.calib_data["aa"].plot(ax=ax0, c="r")
    ax0.set_title("")
    ax0.set_xlim([0, 0.5])

    # Get current AA image from image list
    aa_init = aa_list.current_img()

    # now determine sensitivity correction masks from the different cells
    masks = {}
    aa_imgs_corr = {}
    for cd in cell_aa_calib.cd_vec:
        mask = cellcalib.get_sensitivity_corr_mask("aa",
                                                   pos_x_abs=fov_x,
                                                   pos_y_abs=fov_y,
                                                   radius_abs=fov_extend,
                                                   cell_cd_closest=cd)
        masks[cd] = mask
        aa_imgs_corr[cd] = pyplis.Img(aa_init.img / mask.img)

    # get mask corresponding to minimum cell CD
    mask = list(masks.values())[np.argmin(list(masks.keys()))]

    # assing mask to aa_list
    aa_list.senscorr_mask = mask

    # activate AA sensitivity correction in list
    aa_list.sensitivity_corr_mode = True

    # set DOAS calibration data in list ...
    aa_list.calib_data = doascalib

    # ... and activate calibration mode
    aa_list.calib_mode = True
    ax = aa_list.current_img().show(zlabel=r"$S_{SO2}$ [cm$^{-2}$]")

    # plot the two lines into the exemplary AA image
    pcs1.plot_line_on_grid(ax=ax, color="r")
    pcs2.plot_line_on_grid(ax=ax, color="g")
    ax.legend(loc='best', fancybox=True, framealpha=0.5, fontsize=10)
    ax = draw_doas_fov(fov_x, fov_y, fov_extend, ax=ax)

    fig, _ = plot_pcs_comparison(aa_init, aa_imgs_corr, pcs1, pcs2)

    # IMPORTANT STUFF FINISHED

    if SAVEFIGS:
        ax0.figure.savefig(join(SAVE_DIR, "ex07_out_1.%s" % FORMAT),
                           format=FORMAT, dpi=DPI)
        ax.figure.savefig(join(SAVE_DIR, "ex07_out_2.%s" % FORMAT),
                          format=FORMAT, dpi=DPI)
        fig.savefig(join(SAVE_DIR, "ex07_out_3.%s" % FORMAT), format=FORMAT,
                    dpi=DPI)

    # Save the sensitivity correction mask from the cell with the lowest SO2 CD
    so2min = np.min(list(masks.keys()))
    mask = masks[so2min]
    mask.save_as_fits(SAVE_DIR, "ex07_aa_corr_mask")

    # assign mask in doascalib and resave
    doascalib.senscorr_mask = mask
    doascalib.save_as_fits(DOAS_CALIB_FILE)

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename

        npt.assert_array_equal([],
                               [])

        npt.assert_allclose(actual=[],
                            desired=[],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 8 - Plume velocity retrieval (Cross correlation)

In this script an exemplary plume velocity retrieval is performed using the signal cross correlation algorithm. The velocity is retrieved based on two plume cross sections and a time series of plume AA images (using the AA ImgList created in Example 4 - Preparation of AA image list).

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 8 - Plume velocity retrieval by cross correlation.

"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

from matplotlib.pyplot import close, show, subplots
from os.path import join
from time import time

from pyplis.plumespeed import VeloCrossCorrEngine

# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE, LINES

# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list

# Check script version
check_version()

# SCRIPT OPTONS

# distance in pixels between two lines used for cross correlation analysis
OFFSET_PIXNUM = 40
RELOAD = 0  # reload AA profile images for PCS lines

# start / stop indices of considered images in image list (only relevant if PCS
# profiles are reloaded, i.e. Opt RELOAD=True)
START_IDX = 10
STOP_IDX = 200

# PCS line for which velocity is supposed to retrieved
PCS = LINES[0]  # orange "young_plume" line

# Color of PCS offset line used to perform cross correlation analysis
# (relevant for illustration)
COLOR_PCS_OFFS = "c"

# RELEVANT DIRECTORIES AND PATHS
# the time series of PCS profiles for both lines are stored as
# ProfileTimeSeriesImg objects using the following names. The images
# contain the PCS profiles (y-axis) for each image in the list (y-axis)
# and have thus dimension MXN where M denotes the pixel number of the lines
# and N denotes the total number of images from which the profiles are
# extracted. The images will be stored in SAVE_DIR after this script is run
# once. After that, re-running the script and applying the cross-correlation
# analysis will be much faster, since the profiles are imported from the
# two precomupted images and do not need to be extracted by looping over
# the image list.
PCS_PROFILES_PIC_NAME = "ex08_ica_tseries_pcs.fts"
OFFSET_PROFILES_PIC_NAME = "ex08_ica_tseries_offset.fts"


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    close("all")
    axes = []
    # prepare the AA image list (see ex4)
    aa_list = prepare_aa_image_list()
    aa_list.pyrlevel = 1

    t0 = time()
    cc = VeloCrossCorrEngine(aa_list, PCS)
    cc.create_parallel_pcs_offset(offset_pix=40,
                                  color=COLOR_PCS_OFFS,
                                  linestyle="--")
    reloaded = False  # just a flag for output below
    try:
        if RELOAD:
            raise Exception
        cc.load_pcs_profile_img(join(SAVE_DIR, PCS_PROFILES_PIC_NAME),
                                line_id="pcs")
        cc.load_pcs_profile_img(join(SAVE_DIR, OFFSET_PROFILES_PIC_NAME),
                                line_id="pcs_offset")
    except BaseException:
        cc.get_pcs_tseries_from_list(start_idx=START_IDX,
                                     stop_idx=STOP_IDX)
        cc.save_pcs_profile_images(save_dir=SAVE_DIR,
                                   fname1=PCS_PROFILES_PIC_NAME,
                                   fname2=OFFSET_PROFILES_PIC_NAME)
        reloaded = True
    t1 = time()
    # the run method of the high level VeloCrossCorrEngine class is
    # basically a wrapper method for the low-level find_signal_correlation
    # function which is part of the plumespeed.py module. Before calling
    # the latter, the ICA time-series are extracted from the two
    # ProfileTimeSeriesImg objects which were computed above from the
    # ImgList class containing AA images, and which are stored as FITS
    # files for fast re-computing of this script. The following run
    # command passes valid input parameters to the find_signal_correlation
    # method.
    velo = cc.run(cut_border_idx=10,
                  reg_grid_tres=100,
                  freq_unit="L",
                  sigma_smooth=2,
                  plot=0)
    t2 = time()
    fig, ax = subplots(1, 2, figsize=(20, 6))
    axes.append(cc.plot_pcs_lines(ax=ax[0]))
    cc.plot_ica_tseries_overlay(ax=ax[1])
    axes.append(cc.plot_corrcoeff_tseries())

    print("Result performance analysis\n"
          "Images reloaded from list: %s\n"
          "Number of images: %d\n"
          "Create ICA images: %.3f s\n"
          "Cross-corr analysis: %.3f s"
          % (reloaded, (STOP_IDX - START_IDX), (t1 - t0), (t2 - t1)))

    print("Retrieved plume velocity of v = %.2f m/s" % velo)

    # IMPORTANT STUFF FINISHED
    if SAVEFIGS:
        for k in range(len(axes)):
            axes[k].figure.savefig(join(SAVE_DIR, "ex08_out_%d.%s"
                                        % ((k + 1), FORMAT)),
                                   format=FORMAT, dpi=DPI)

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename

        npt.assert_array_equal([],
                               [])

        npt.assert_allclose(actual=[],
                            desired=[],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 9 - Plume velocity retrieval (Optical flow Farneback)

This script gives an introduction into plume velocity retrievals using the Farneback optical flow algorithm (OpticalFlowFarneback) and a histogram based post analysis.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 9 - Optical flow Plume velocity retrieval."""
from __future__ import (absolute_import, division)

from os.path import join
import pyplis
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE, LINES
from matplotlib.pyplot import (close, show, subplots, figure, xticks, yticks,
                               sca, rcParams)

# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list
from SETTINGS import check_version

rcParams["font.size"] = 16
PCS1, PCS2 = LINES

# Check script version
check_version()

# SCRIPT OPTIONS

PEAK_SIGMA_TOL = 2

# perform histogram analysis for all images in time series
HISTO_ANALYSIS_ALL = 1
# applies multi gauss fit to retrieve local predominant displacement
# direction, if False, then the latter is calculated from 1. and 2. moment
# of histogram (Faster but more sensitive to additional peaks in histogram)
HISTO_ANALYSIS_MULTIGAUSS = True
HISTO_ANALYSIS_START_IDX = 0
HISTO_ANALYSIS_STOP_IDX = None  # 207

# Gauss pyramid level
PYRLEVEL = 1
BLUR = 0
ROI_CONTRAST = [0, 0, 1344, 730]
MIN_AA = 0.05


def analyse_and_plot(lst, lines):
    fig = figure(figsize=(14, 8))

    ax0 = fig.add_axes([0.01, 0.15, 0.59, 0.8])
    # ax0.set_axis_off()
    ax1 = fig.add_axes([0.61, 0.15, 0.16, 0.8])
    ax2 = fig.add_axes([0.78, 0.15, 0.16, 0.8])
    mask = lst.get_thresh_mask(MIN_AA)
    fl = lst.optflow
    fl.plot(ax=ax0)  # , in_roi=True)
    for line in lines:
        m = mask * line.get_rotated_roi_mask(fl.flow.shape[:2])
        line.plot_line_on_grid(ax=ax0, include_normal=1,
                               include_roi_rot=1)
        try:
            _, mu, sigma = fl.plot_orientation_histo(pix_mask=m,
                                                     apply_fit=True, ax=ax1,
                                                     color=line.color)
            ax1.legend_.remove()
            low, high = mu - sigma, mu + sigma
            fl.plot_length_histo(pix_mask=m, ax=ax2, dir_low=low,
                                 dir_high=high, color=line.color)
        except BaseException:
            pass
    # pyplis.helpers.set_ax_lim_roi(roi_disp, ax0)
    ax0.get_xaxis().set_ticks([])
    ax0.get_yaxis().set_ticks([])
    ax0.set_title("")

    ymax = max([ax1.get_ylim()[1], ax2.get_ylim()[1]])
    ax1.set_title("")
    ax1.set_xlabel(r"$\varphi\,[^\circ]$", fontsize=20)
    ax1.set_ylim([0, ymax])
    ax1.get_yaxis().set_ticks([])
    ax2.yaxis.tick_right()
    ax2.yaxis.set_label_position("right")
    ax2.set_xlabel(r"$|\mathbf{f}|$ [pix]", fontsize=20)
    ax2.set_ylabel("Counts / bin", fontsize=20)
    ax2.set_ylim([0, ymax])

    ax2.set_title("")
    ax2.legend_.remove()

    sca(ax1)
    xticks(rotation=40, ha="right")
    sca(ax2)
    yticks(rotation=90, va="center")

    return fig


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    close("all")
    figs = []
    # Prepare aa image list (see example 4)
    aa_list = prepare_aa_image_list()

    # the aa image list includes the measurement geometry, get pixel
    # distance image where pixel values correspond to step widths in the plume,
    # obviously, the distance values depend on the downscaling factor, which
    # is calculated from the analysis pyramid level (PYRLEVEL)
    dist_img, _, _ = aa_list.meas_geometry.compute_all_integration_step_lengths(  # noqa: E501
        pyrlevel=PYRLEVEL)
    # set the pyramid level in the list
    aa_list.pyrlevel = PYRLEVEL
    # add some blurring.. or not (if BLUR = 0)
    aa_list.add_gaussian_blurring(BLUR)

    # Access to the optical flow module in the image list. If optflow_mode is
    # active in the list, then, whenever the list index changes (e.g. using
    # list.goto_next(), or list.goto_img(100)), the optical flow field is
    # calculated between the current list image and the next one
    fl = aa_list.optflow
    # (! note: fl is only a pointer, i.e. the "=" is not making a copy of the
    # object, meaning, that whenever something changes in "fl", it also does
    # in "aa_list.optflow")

    # Now activate optical flow calculation in list (this slows down the
    # speed of the analysis, since the optical flow calculation is
    # comparatively slow
    s = aa_list.optflow.settings
    s.hist_dir_gnum_max = 10
    s.hist_dir_binres = 10
    s.hist_sigma_tol = PEAK_SIGMA_TOL

    s.roi_rad = ROI_CONTRAST

    aa_list.optflow_mode = True

    plume_mask = pyplis.Img(aa_list.get_thresh_mask(MIN_AA))
    plume_mask.show(tit="AA threshold mask")

    figs.append(analyse_and_plot(aa_list, LINES))

    figs.append(fl.plot_flow_histograms(PCS1, plume_mask.img))
    figs.append(fl.plot_flow_histograms(PCS2, plume_mask.img))

    # Show an image containing plume speed magnitudes (ignoring direction)
    velo_img = pyplis.Img(fl.to_plume_speed(dist_img))
    velo_img.show(vmin=0, vmax=10, cmap="Greens",
                  tit="Optical flow plume velocities",
                  zlabel="Plume velo [m/s]")

    # Create two objects used to store time series information about the
    # retrieved plume properties
    plume_props_l1 = pyplis.plumespeed.LocalPlumeProperties(PCS1.line_id)
    plume_props_l2 = pyplis.plumespeed.LocalPlumeProperties(PCS2.line_id)

    if HISTO_ANALYSIS_ALL:
        aa_list.goto_img(HISTO_ANALYSIS_START_IDX)
        if HISTO_ANALYSIS_STOP_IDX is None:
            HISTO_ANALYSIS_STOP_IDX = aa_list.nof - 1
        for k in range(HISTO_ANALYSIS_START_IDX, HISTO_ANALYSIS_STOP_IDX):
            plume_mask = aa_list.get_thresh_mask(MIN_AA)
            plume_props_l1.get_and_append_from_farneback(
                fl, line=PCS1, pix_mask=plume_mask,
                dir_multi_gauss=HISTO_ANALYSIS_MULTIGAUSS)
            plume_props_l2.get_and_append_from_farneback(
                fl, line=PCS2, pix_mask=plume_mask,
                dir_multi_gauss=HISTO_ANALYSIS_MULTIGAUSS)
            aa_list.goto_next()

# ==============================================================================
#         plume_props_l1 = plume_props_l1.interpolate()
#         plume_props_l2 = plume_props_l2.interpolate()
# ==============================================================================

        fig, ax = subplots(2, 1, figsize=(10, 9))

        plume_props_l1.plot_directions(ax=ax[0],
                                       color=PCS1.color,
                                       label="PCS1")
        plume_props_l2.plot_directions(ax=ax[0], color=PCS2.color,
                                       label="PCS2")

        plume_props_l1.plot_magnitudes(normalised=True, ax=ax[1],
                                       date_fmt="%H:%M:%S", color=PCS1.color,
                                       label="PCS1")
        plume_props_l2.plot_magnitudes(normalised=True, ax=ax[1],
                                       date_fmt="%H:%M:%S", color=PCS2.color,
                                       label="PCS2")
        ax[0].set_xticklabels([])
        # ax[0].legend(loc='best', fancybox=True, framealpha=0.5, fontsize=14)
        # ax[0].set_title("Movement direction")
        # ax[1].set_title("Displacement length")
        figs.append(fig)
        # Save the time series as txt
        plume_props_l1.save_txt(join(SAVE_DIR,
                                     "ex09_plumeprops_young_plume.txt"))
        plume_props_l2.save_txt(join(SAVE_DIR,
                                     "ex09_plumeprops_aged_plume.txt"))

    if SAVEFIGS:
        for k in range(len(figs)):
            figs[k].savefig(join(SAVE_DIR, "ex09_out_%d.%s"
                                 % ((k + 1), FORMAT)),
                            format=FORMAT, dpi=DPI)

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename

        npt.assert_array_equal([],
                               [])

        npt.assert_allclose(actual=[],
                            desired=[],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 10 - Import plume background images

Create a Dataset object for a time interval containing only plume background images (on / off).

Note

Stand alone script that is not required for any of the following scripts

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 10 - Create background image dataset."""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

import pyplis
from datetime import datetime
from matplotlib.pyplot import show

# IMPORT GLOBAL SETTINGS
from SETTINGS import IMG_DIR, OPTPARSE

# Check script version
check_version()

# SCRIPT FUNCTION DEFINITIONS


def get_bg_image_lists():
    """Initialize measurement setup and creates dataset from that."""
    start = datetime(2015, 9, 16, 7, 2, 0o5)
    stop = datetime(2015, 9, 16, 7, 2, 30)
    # Define camera (here the default ecII type is used)
    cam_id = "ecII"

    # the camera filter setup
    filters = [pyplis.utils.Filter(type="on", acronym="F01"),
               pyplis.utils.Filter(type="off", acronym="F02")]

    # create camera setup
    cam = pyplis.setupclasses.Camera(cam_id=cam_id, filter_list=filters)

    # Create BaseSetup object (which creates the MeasGeometry object)
    stp = pyplis.setupclasses.MeasSetup(IMG_DIR, start, stop, camera=cam)

    ds = pyplis.dataset.Dataset(stp)
    on, off = ds.get_list("on"), ds.get_list("off")
    on.darkcorr_mode = True
    off.darkcorr_mode = True
    return on, off


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    on, off = get_bg_image_lists()
    on.show_current()
    off.show_current()

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename

        npt.assert_array_equal([],
                               [])

        npt.assert_allclose(actual=[],
                            desired=[],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 11 - Image based signal dilution correction

This script introduces the image based signal dilution correction including automatic retrieval of terrain distances on a pixel basis.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 11 - Image based signal dilution correction.

This script illustrates how extinction coefficients can be retrieved from
image data using the DilutionCorr class and by specifying suitable terrain
features in the images.

The extinction coefficients are retrieved from one on an from one off band
image recorded ~15 mins before the dataset used in the other examples. The
latter data is less suited for illustrating the feature since it contains
less terrain (therfore more sky background).

The two example images are then corrected for dilution and the results are
plotted (as comparison of the retrieved emission rate along an exemplary
plume cross section)
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

import pyplis as pyplis
from geonum import GeoPoint
from matplotlib.pyplot import show, close, subplots, Rectangle, plot
from datetime import datetime
from os.path import join, exists

from pyplis.dilutioncorr import DilutionCorr
from pyplis.doascalib import DoasCalibData

# IMPORT GLOBAL SETTINGS
from SETTINGS import IMG_DIR, SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex10_bg_imglists import get_bg_image_lists

# Check script version
check_version()

# SCRIPT OPTONS
# lower boundary for I0 value in dilution fit
I0_MIN = 0.0
AA_THRESH = 0.03

# exemplary plume cross section line for emission rate retrieval (is also used
# for full analysis in ex12)
PCS_LINE = pyplis.LineOnImage(x0=530, y0=586, x1=910, y1=200, line_id="pcs")
# Retrieval lines for dilution correction (along these lines, topographic
# distances and image radiances are determined for fitting the atmospheric
# extinction coefficients)
TOPO_LINE1 = pyplis.LineOnImage(1100, 650, 1000, 900, line_id="flank far",
                                color="lime",
                                linestyle="-")

TOPO_LINE2 = pyplis.LineOnImage(1000, 990, 1100, 990, line_id="flank close",
                                color="#ff33e3",
                                linestyle="-")

# all lines in this array are used for the analysis
USE_LINES = [TOPO_LINE1, TOPO_LINE2]

# specify pixel resolution of topographic distance retrieval (every nth pixel
# is used)
SKIP_PIX_LINES = 10

# Specify region of interest used to extract the ambient intensity (required
# for dilution correction)
AMBIENT_ROI = [1240, 10, 1300, 70]

# Specify plume velocity (for emission rate estimate)
PLUME_VELO = 4.14  # m/s (result from ex8)

# RELEVANT DIRECTORIES AND PATHS
CALIB_FILE = join(SAVE_DIR, "ex06_doascalib_aa.fts")

# SCRIPT FUNCTION DEFINITIONS


def create_dataset_dilution():
    """Create a :class:`pyplis.dataset.Dataset` object for dilution analysis.

    The test dataset includes one on and one offband image which are recorded
    around 6:45 UTC at lower camera elevation angle than the time series shown
    in the other examples (7:06 - 7:22 UTC). Since these two images contain
    more topographic features they are used to illustrate the image based
    signal dilution correction.

    This function sets up the measurement (geometry, camera, time stamps) for
    these two images and creates a Dataset object.
    """
    start = datetime(2015, 9, 16, 6, 43, 00)
    stop = datetime(2015, 9, 16, 6, 47, 00)
    # the camera filter setup
    cam_id = "ecII"
    filters = [pyplis.utils.Filter(type="on", acronym="F01"),
               pyplis.utils.Filter(type="off", acronym="F02")]

    geom_cam = {"lon": 15.1129,
                "lat": 37.73122,
                "elev": 15.0,  # from field notes, will be corrected
                "elev_err": 5.0,
                "azim": 274.0,  # from field notes, will be corrected
                "azim_err": 10.0,
                "focal_length": 25e-3,
                "alt_offset": 7}  # meters above topography

    # create camera setup
    cam = pyplis.setupclasses.Camera(cam_id=cam_id, filter_list=filters,
                                     **geom_cam)

    # Load default information for Etna
    source = pyplis.setupclasses.Source("etna")

    # Provide wind direction
    wind_info = {"dir": 0.0,
                 "dir_err": 15.0}

    # Create BaseSetup object (which creates the MeasGeometry object)
    stp = pyplis.setupclasses.MeasSetup(IMG_DIR, start, stop, camera=cam,
                                        source=source,
                                        wind_info=wind_info)
    return pyplis.dataset.Dataset(stp)


def find_view_dir(geom):
    """Perform a correction of the viewing direction using crater in img.

    :param MeasGeometry geom: measurement geometry
    :param str which_crater: use either "ne" (northeast) or "se" (south east)
    :return: - MeasGeometry, corrected geometry
    """
    # Use position of NE crater in image
    posx, posy = 1051, 605  # pixel position of NE crate in image
    # Geo location of NE crater (info from Google Earth)
    ne_crater = GeoPoint(37.754788, 14.996673, 3287, name="NE crater")

    geom.find_viewing_direction(pix_x=posx, pix_y=posy, pix_pos_err=100,
                                geo_point=ne_crater, draw_result=True)
    return geom


def prepare_lists(dataset):
    """Prepare on and off lists for dilution analysis.

    Steps:

        1. get on and offband list
        #. load background image list on and off (from ex10)
        #. set image preparation and assign background images to on / off list
        #. configure plume background model settings

    :param Dataset dataset: the dilution dataset (see
        :func:`create_dataset_dilution`)
    :return:
        - ImgList, onlist
        - ImgList, offlist

    """
    onlist = dataset.get_list("on")
    offlist = dataset.get_list("off")
    # dark_corr_mode already active
    bg_onlist, bg_offlist = get_bg_image_lists()

    # prepare img pre-edit
    onlist.darkcorr_mode = True
    onlist.gaussian_blurring = 2
    offlist.darkcorr_mode = True
    offlist.gaussian_blurring = 2

    # prepare background images in lists (when assigning a background image
    # to a ImgList, then the blurring amount of the BG image is automatically
    # set to the current blurring level of the list, and if the latter is
    # zero, then the BG image is blurred using filter width=1)
    onlist.bg_img = bg_onlist.current_img()
    offlist.bg_img = bg_offlist.current_img()

    # prepare plume background modelling setup in both lists
    onlist.bg_model.mode = 6
    onlist.bg_model.set_missing_ref_areas(onlist.current_img())
    onlist.bg_model.xgrad_line_startcol = 10
    offlist.bg_model.update(**onlist.bg_model.settings_dict())

    return onlist, offlist


def plot_retrieval_points_into_image(img):
    """Plot terrain distance retrieval lines into an image."""
    ax = img.show(vmin=-5e16, vmax=5e18, zlabel=r"$S_{SO2}$ [cm$^{-2}$]")
    ax.set_title("Retrieval lines")
    for line in USE_LINES:
        line.plot_line_on_grid(ax=ax, marker="", color=line.color,
                               lw=2, ls=line.linestyle)
    for x, y, dist in dil._add_points:
        plot(x, y, " or")
    return ax


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    if not exists(CALIB_FILE):
        raise IOError("Calibration file could not be found at specified "
                      "location:\n %s\nYou might need to run example 6 first")

    close("all")
    pcs_line = PCS_LINE
    calib = DoasCalibData()
    calib.load_from_fits(CALIB_FILE)

    # create dataset and correct viewing direction
    ds = create_dataset_dilution()
    geom = find_view_dir(ds.meas_geometry)

    # get plume distance image
    pix_dists, _, plume_dists = geom.compute_all_integration_step_lengths()

    # Create dilution correction class
    dil = DilutionCorr(USE_LINES, geom, skip_pix=SKIP_PIX_LINES)

    # you can also manually add  a single pixel position in the image that is
    # used to retrieve topographic distances (this function allows you also to
    # set the distance to the feature manually, since sometimes the SRTM data-
    # set is incomplete or has too low resolution)
    dil.add_retrieval_point(700, 930)
    dil.add_retrieval_point(730, 607)

    # Determine distances to the two lines defined above (every 6th pixel)
    for line_id in dil.line_ids:
        dil.det_topo_dists_line(line_id)

    # Plot the results in a 3D map
    basemap = dil.plot_distances_3d(alt_offset_m=10, axis_off=False)

    # retrieve pixel distances for pixels on the line
    # (for emission rate estimate)
    pix_dists_line = pcs_line.get_line_profile(pix_dists)

    # get pixel coordinates of PCS center position ...
    col, row = pcs_line.center_pix

    # ... and get uncertainty in plume distance estimate for the column
    pix_dist_err = geom.pix_dist_err(col)

    # Prepare on and off-band list for retrieval of extinction coefficients
    onlist, offlist = prepare_lists(ds)

    # Activate vignetting correction in both lists (required to extract
    # measured intensities along the terrain features)
    onlist.vigncorr_mode = True
    offlist.vigncorr_mode = True

    # get the vignetting corrected images
    on_vigncorr = onlist.current_img()
    off_vigncorr = offlist.current_img()

    # estimate ambient intensity for both filters
    ia_on = on_vigncorr.crop(AMBIENT_ROI, True).mean()
    ia_off = off_vigncorr.crop(AMBIENT_ROI, True).mean()

    # perform dilution anlysis and retrieve extinction coefficients (on-band)
    ext_on, _, _, ax0 = dil.apply_dilution_fit(img=on_vigncorr,
                                               rad_ambient=ia_on,
                                               i0_min=I0_MIN,
                                               plot=True)

    ax0.set_ylabel("Terrain radiances (on band)", fontsize=14)
    ax0.set_ylim([0, 2500])

    # perform dilution anlysis and retrieve extinction coefficients (off-band)
    ext_off, i0_off, _, ax1 = dil.apply_dilution_fit(img=off_vigncorr,
                                                     rad_ambient=ia_off,
                                                     i0_min=I0_MIN,
                                                     plot=True)
    ax1.set_ylabel("Terrain radiances (off band)", fontsize=14)
    ax1.set_ylim([0, 2500])

    # determine plume pixel mask from AA image
    onlist.aa_mode = True

    plume_pix_mask = onlist.get_thresh_mask(AA_THRESH)
    plume_pix_mask[840:, :] = 0  # remove tree in lower part of the image
    onlist.aa_mode = False

    # assign the just retrieved extinction coefficients to the respective
    # image lists
    onlist.ext_coeffs = ext_on
    offlist.ext_coeffs = ext_off

    # save the extinction coefficients into a txt file (re-used in example
    # script 12). They are stored as pandas.Series object in the ImgList
    onlist.ext_coeffs.to_csv(join(SAVE_DIR, "ex11_ext_scat_on.txt"))
    offlist.ext_coeffs.to_csv(join(SAVE_DIR, "ex11_ext_scat_off.txt"))

    # now activate automatic dilution correction in both lists
    # get dilution corrected on and off-band image
    onlist.dilcorr_mode = True
    offlist.dilcorr_mode = True

    # get current dilution corrected raw images (i.e. in intensity space)
    on_corr = onlist.this
    off_corr = offlist.this

    # now activate tau mode (note that dilution correction mode is still
    # active)
    onlist.tau_mode = True
    offlist.tau_mode = True

    # extract tau images
    tau_on_corr = onlist.this
    tau_off_corr = offlist.this

    # determine corrected SO2-CD image from the image lists
    so2_img_corr = calib(tau_on_corr - tau_off_corr)
    so2_img_corr.edit_log["is_tau"] = True  # for plotting
    so2_cds_corr = pcs_line.get_line_profile(so2_img_corr)

    (phi_corr,
     phi_corr_err) = pyplis.fluxcalc.det_emission_rate(
        cds=so2_cds_corr,
        velo=PLUME_VELO,
        pix_dists=pix_dists_line,
        cds_err=calib.err(),
        pix_dists_err=pix_dist_err)

    # determine uncorrected so2-CD image from the image lists
    offlist.dilcorr_mode = False
    onlist.dilcorr_mode = False

    # the "this" attribute returns the current list image (same as
    # method "current_img()")
    so2_img_uncorr = calib(onlist.this - offlist.this)

    # Retrieve column density profile along PCS in uncorrected image
    so2_cds_uncorr = pcs_line.get_line_profile(so2_img_uncorr)
    # Calculate flux and uncertainty
    (phi_uncorr,
     phi_uncorr_err) = pyplis.fluxcalc.det_emission_rate(
        cds=so2_cds_uncorr,
        velo=PLUME_VELO,
        pix_dists=pix_dists_line,
        cds_err=calib.err(),
        pix_dists_err=pix_dist_err)

    # IMPORTANT STUFF FINISHED (below follow some plots)
    ax2 = plot_retrieval_points_into_image(so2_img_corr)
    pcs_line.plot_line_on_grid(ax=ax2, ls="-", color="g")
    ax2.legend(loc="best", framealpha=0.5, fancybox=True, fontsize=20)
    ax2.set_title("Dilution corrected AA image", fontsize=12)
    ax2.get_xaxis().set_ticks([])
    ax2.get_yaxis().set_ticks([])

    x0, y0, w, h = pyplis.helpers.roi2rect(AMBIENT_ROI)
    ax2.add_patch(Rectangle((x0, y0), w, h, fc="none", ec="c"))

    fig, ax3 = subplots(1, 1)
    ax3.plot(so2_cds_uncorr, ls="-", color="#ff33e3",
             label=r"Uncorr: $\Phi_{SO2}=$%.2f (+/- %.2f) kg/s"
                   % (phi_uncorr / 1000.0, phi_uncorr_err / 1000.0))
    ax3.plot(so2_cds_corr, "-g", lw=3,
             label=r"Corr: $\Phi_{SO2}=$%.2f (+/- %.2f) kg/s"
                   % (phi_corr / 1000.0, phi_corr_err / 1000.0))

    ax3.set_title("Cross section profile", fontsize=12)
    ax3.legend(loc="best", framealpha=0.5, fancybox=True, fontsize=12)
    ax3.set_xlim([0, len(pix_dists_line)])
    ax3.set_ylim([0, 5e18])
    ax3.set_ylabel(r"$S_{SO2}$ [cm$^{-2}$]", fontsize=14)
    ax3.set_xlabel("PCS", fontsize=14)
    ax3.grid()

    # also plot plume pixel mask
    ax4 = pyplis.Img(plume_pix_mask).show(cmap="gray", tit="Plume pixel mask")

    if SAVEFIGS:
        ax = [ax0, ax1, ax2, ax3, ax4]
        for k in range(len(ax)):
            ax[k].set_title("")  # remove titles for saving
            ax[k].figure.savefig(join(SAVE_DIR, "ex11_out_%d.%s"
                                 % (k, FORMAT)),
                                 format=FORMAT, dpi=DPI)
        basemap.ax.set_axis_off()
        basemap.ax.view_init(15, 345)
        basemap.ax.figure.savefig(join(SAVE_DIR, "ex11_out_5.%s" % FORMAT),
                                  format=FORMAT, dpi=DPI)

# IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename

        npt.assert_array_equal([],
                               [])

        npt.assert_allclose(actual=[],
                            desired=[],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")

Example 12 - Emission rate analysis (Etna example data)

Perform emission rate analysis for the example data. The analysis is performed along one plume cross section (in the image center) and using three different plume velocity retrievals.

Code

# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script 12 - Etna emission rate retrieval.

This example import results from the previous examples, for instance the AA
image list including measurement geometry (ex 4), the DOAS calibration
information (which was stored as FITS file, see ex. 6) and the AA sensitivity
correction mask retrieved from the cell calibration and normalised to the
position of the DOAS FOV (ex 7). The emission rates are retrieved for three
different plume velocity retrievals: 1. using the global velocity vector
retrieved from the cross correlation algorithm (ex8), 2. using the raw output
of the optical flow Farneback algorithm (``flow_raw``) and 3. using the
histogram based post analysis of the optical flow field (``flow_histo``).
The analysis is performed using the EmissionRateAnalysis class which basically
checks the AA list and activates ``calib_mode`` (-> images are loaded as
calibrated gas CD images) and loops over all images to retrieve the emission
rates for the 3 velocity modes. Here, emission rates are retrieved along 1
exemplary plume cross section. This can be easily extended by adding additional
PCS lines in the EmissionRateAnalysis class using ``add_pcs_line``.
The results for each velocity mode and for each PCS line are stored within
EmissionRateResults classes.
"""
from __future__ import (absolute_import, division)

from SETTINGS import check_version

from os.path import join, exists
from matplotlib.pyplot import close, show, GridSpec, figure, rc_context
from matplotlib.cm import get_cmap

import pyplis
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE, LINES

# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list

rc_context({'font.size': '18'})

# Check script version
check_version()

PCS = LINES[0]

# If false, then only the working environment is initalised
DO_EVAL = True

# Dilution correction
DILCORR = True

# You can specify here if you only want a certain number of images analysed
START_INDEX = 0
STOP_INDEX = None  # 20


# SCRIPT OPTONS
PYRLEVEL = 1
PLUME_VELO_GLOB = 4.29  # m/s
PLUME_VELO_GLOB_ERR = 1.5
# applies multi gauss fit to retrieve local predominant displacement
# direction, if False, then the latter is calculated from 1. and 2. moment
# of histogram (Faster but more sensitive to additional peaks in histogram)
HISTO_ANALYSIS_MULTIGAUSS = True
# molar mass of SO2
MMOL = 64.0638  # g/mol
# minimum required SO2-CD for emission-rate retrieval
CD_MIN = 5e16

# activate background check mode, if True, emission rates are only
# retrieved  for images showing SO2-CDs within specified interval around
# zero in BG reference rectangle LOG_ROI_SKY (see above). This can be
# used to ensure that significant systematic errors are induced in case
# the plume background retrieval failed. The latter could, for instance
# happen, if, for instance a cloud moves through one of the background
# reference areas used to model the background (cf. example script 3)
REF_CHECK_LOWER = -5e16
REF_CHECK_UPPER = 5e16
REF_CHECK_MODE = True

# the following ROI is in the upper right image corner, where no gas occurs in
# the time series. It is used to log mean, min and max for each analysed image
# this information can be used to check, whether the plume background retrieval
# worked well
LOG_ROI_SKY = [530, 30, 600, 100]  # correspond to pyrlevel 1

# RELEVANT DIRECTORIES AND PATHS

# DOAS calibration results from example script 6
CALIB_FILE = join(SAVE_DIR, "ex06_doascalib_aa.fts")

# Scattering extinction coeffcients from example script 11 (stored as txt)
EXT_ON = join(SAVE_DIR, "ex11_ext_scat_on.txt")
EXT_OFF = join(SAVE_DIR, "ex11_ext_scat_off.txt")

# AA sensitivity correction mask retrieved from cell calib in script 7
CORR_MASK_FILE = join(SAVE_DIR, "ex07_aa_corr_mask.fts")

# time series of predominant displacement vector from histogram analysis of
# optical flow field in ROI around the PCS line "young_plume" which is used
# here for the emission rate retrieval. These information is optional, and is
# calculated during the evaluation if not provided
RESULT_PLUMEPROPS_HISTO = join(SAVE_DIR, "ex09_plumeprops_young_plume.txt")

# SCRIPT FUNCTION DEFINITIONS


def plot_and_save_results(ana, line_id="young_plume", date_fmt="%H:%M"):

    # plot colors for different optical flow retrievals
    cmap = get_cmap("Oranges")

    c_optflow_hybrid = cmap(255)
    c_optflow_histo = cmap(175)
    c_optflow_raw = cmap(100)

    fig = figure(figsize=(16, 12))
    gs = GridSpec(4, 1, height_ratios=[.6, .2, .2, .2], hspace=0.05)
    ax3 = fig.add_subplot(gs[3])
    ax0 = fig.add_subplot(gs[0], sharex=ax3)
    ax1 = fig.add_subplot(gs[1], sharex=ax3)
    ax2 = fig.add_subplot(gs[2], sharex=ax3)
    ax1.yaxis.tick_right()
    ax1.yaxis.set_label_position("right")
    ax3.yaxis.tick_right()
    ax3.yaxis.set_label_position("right")

    # Get emission rate results for the PCS line
    res0 = ana.get_results(line_id=line_id, velo_mode="glob")
    res1 = ana.get_results(line_id=line_id, velo_mode="flow_raw")
    res2 = ana.get_results(line_id=line_id, velo_mode="flow_histo")
    res3 = ana.get_results(line_id=line_id, velo_mode="flow_hybrid")

    res0.save_txt(join(SAVE_DIR, "ex12_flux_velo_glob.txt"))
    res1.save_txt(join(SAVE_DIR, "ex12_flux_flow_raw.txt"))
    res2.save_txt(join(SAVE_DIR, "ex12_flux_flow_histo.txt"))
    res3.save_txt(join(SAVE_DIR, "ex12_flux_flow_hybrid.txt"))

    # Plot emission rates for the different plume speed retrievals
    res0.plot(yerr=True, date_fmt=date_fmt, ls="-", ax=ax0,
              color="c", ymin=0, alpha_err=0.08)
    res1.plot(yerr=False, ax=ax0, ls="-", color=c_optflow_raw, ymin=0)
    res2.plot(yerr=False, ax=ax0, ls="--", color=c_optflow_histo, ymin=0)
    res3.plot(yerr=True, ax=ax0, lw=3, ls="-", color=c_optflow_hybrid, ymin=0)

    # ax[0].set_title("Retrieved emission rates")
    ax0.legend(loc='best', fancybox=True, framealpha=0.5, fontsize=12)
    ax0.grid()

    # Plot effective velocity retrieved from optical flow histogram analysis
    res3.plot_velo_eff(ax=ax1, date_fmt=date_fmt, color=c_optflow_hybrid)
    # ax[1].set_title("Effective plume speed
    #                  (from optflow histogram analysis)")
    ax1.set_ylim([0, ax1.get_ylim()[1]])

    # Plot time series of predominant plume direction (retrieved from optical
    # flow histogram analysis and stored in object of type LocalPlumeProperties
    # which is part of plumespeed.py module
    ana.pcs_lines[line_id].plume_props.plot_directions(ax=ax2,
                                                       date_fmt=date_fmt,
                                                       color=c_optflow_hybrid)

    ax2.set_ylim([-180, 180])
    pyplis.helpers.rotate_xtick_labels(ax=ax2)
    ax0.set_xticklabels([])
    ax1.set_xticklabels([])
    ax2.set_xticklabels([])
    # tight_layout()

    ax3 = ana.plot_bg_roi_vals(ax=ax3, date_fmt="%H:%M")
    # gs.tight_layout(fig, h_pad=0)#0.03)
    gs.update(hspace=0.05, top=0.97, bottom=0.07)
    return fig


# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
    close("all")
    figs = []
    if not exists(CALIB_FILE):
        raise IOError("Calibration file could not be found at specified "
                      "location:\n%s\nPlease run example 6 first")
    if not exists(CORR_MASK_FILE):
        raise IOError("Cannot find AA correction mask, please run example"
                      "script 7 first")

    # convert the retrieval line to the specified pyramid level (script option)
    pcs = PCS.convert(to_pyrlevel=PYRLEVEL)

    # now try to load results of optical flow histogram analysis performed for
    # this line in script no. 9. and assign them to the pcs line. This has the
    # advantage, that  missing velocity vectors (i.e. from images where optical
    # flow analysis failed) can be interpolated. It is, however, not
    # necessarily required to do this in advance. In the latter case the
    # emission rates show gaps at all images, where the optical flow was
    # considered not reliable
    try:
        p = pyplis.LocalPlumeProperties()
        p.load_txt(RESULT_PLUMEPROPS_HISTO)
        p = p.to_pyrlevel(PYRLEVEL)
        fig = p.plot(color="r")
        # p.interpolate()
        # p = p.apply_significance_thresh(0.2).interpolate()
        # p = p.apply_median_filter(3).apply_gauss_filter(2)
        fig = p.plot(date_fmt="%H:%M", fig=fig)

        pcs.plume_props = p
    except BaseException:
        print("Local plume properties could not be loaded and will be "
              "calculated during the emission rate analysis")

    # Load AA list
    # includes viewing direction corrected geometry
    aa_list = prepare_aa_image_list()

    aa_list.pyrlevel = PYRLEVEL

    if DILCORR:
        aa_list.import_ext_coeffs_csv(EXT_ON)
        aa_list.get_off_list().import_ext_coeffs_csv(EXT_OFF)

    # Load DOAS calbration data and FOV information (see example 6)
    doascalib = pyplis.doascalib.DoasCalibData()
    doascalib.load_from_fits(file_path=CALIB_FILE)
    doascalib.fit_calib_data()

    # Load AA corr mask and set in image list(is normalised to DOAS FOV see
    # ex7)
    aa_corr_mask = pyplis.Img(CORR_MASK_FILE)

    aa_list.senscorr_mask = aa_corr_mask

    # set DOAS calibration data in image list
    aa_list.calib_data = doascalib

    ana = pyplis.EmissionRateAnalysis(
        imglist=aa_list,
        bg_roi=LOG_ROI_SKY,
        pcs_lines=pcs,
        velo_glob=PLUME_VELO_GLOB,
        velo_glob_err=PLUME_VELO_GLOB_ERR,
        ref_check_lower_lim=REF_CHECK_LOWER,
        ref_check_upper_lim=REF_CHECK_UPPER,
        velo_dir_multigauss=HISTO_ANALYSIS_MULTIGAUSS,
        senscorr=True,
        dilcorr=DILCORR)

    ana.settings.ref_check_mode = REF_CHECK_MODE

    ana.settings.velo_modes["flow_raw"] = 1
    ana.settings.velo_modes["flow_histo"] = True
    ana.settings.velo_modes["flow_hybrid"] = 1
    ana.settings.min_cd = CD_MIN

    # plot all current PCS lines into current list image (feel free to define
    # and add more PCS lines above)
    ax = ana.plot_pcs_lines(
        vmin=-
        5e18,
        vmax=6e18,
        tit="Dilution corr: %s" %
        DILCORR)
    ax = ana.plot_bg_roi_rect(ax=ax, to_pyrlevel=PYRLEVEL)
    figs.append(ax.figure)

    if not DO_EVAL:
        aa_list.dilcorr_mode = not DILCORR
        aa_list.show_current(
            vmin=-
            5e18,
            vmax=6e18,
            tit="Dilution corr: %s" %
            (not DILCORR))
        # you can check the settings first
        print(ana.settings)
        # check if optical flow works
        ana.imglist.optflow_mode = True
        aa_mask = ana.imglist.get_thresh_mask(CD_MIN)
        ana.imglist.optflow.plot_flow_histograms(line=pcs, pix_mask=aa_mask)

    else:
        ana.run_retrieval(start_index=START_INDEX,
                          stop_index=STOP_INDEX)

        figs.append(plot_and_save_results(ana))
        # the EmissionRateResults class has an informative string
        # representation
        print(ana.get_results("young_plume", "flow_histo"))

    if SAVEFIGS:
        for k in range(len(figs)):
            figs[k].savefig(join(SAVE_DIR, "ex12_out_%d.%s" % (k + 1, FORMAT)),
                            format=FORMAT, dpi=DPI)

    # IMPORTANT STUFF FINISHED (Below follow tests and display options)

    # Import script options
    (options, args) = OPTPARSE.parse_args()

    # If applicable, do some tests. This is done only if TESTMODE is active:
    # testmode can be activated globally (see SETTINGS.py) or can also be
    # activated from the command line when executing the script using the
    # option --test 1
    if int(options.test):
        import numpy.testing as npt
        from os.path import basename

        npt.assert_array_equal([],
                               [])

        npt.assert_allclose(actual=[],
                            desired=[],
                            rtol=1e-7)
        print("All tests passed in script: %s" % basename(__file__))
    try:
        if int(options.show) == 1:
            show()
    except BaseException:
        print("Use option --show 1 if you want the plots to be displayed")