Example scripts
Pyplis example scripts. The scripts require downloading the pyplis test data. A selection of outputs from the example scripts is shown in the Plot gallery.
The examples are build upon each other and are categorized into a set of basic examples to get started with the core features of pyplis and a set of advanced examples that introduce the main analysis steps and available routines to perform SO2 emission rate retrievals, including all aspects of the analysis workflow.
Basic examples for getting started
These scripts give an introduction into basic features and classes of pyplis. More advanced examples tailored to SO2 emission rate analysis are provided below, in the second part of this document.
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 datetime import datetime
import matplotlib.pyplot as plt
import pathlib
import pyplis
# imports from SETTINGS.py
from SETTINGS import ARGPARSER, SAVE_DIR
# file name of test image stored in data folder
IMG_FILE_NAME = "test_201509160708_F01_335.fts"
IMG_DIR = pathlib.Path(pyplis.__dir__) / "data"
def main():
plt.close("all")
img_path = 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(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 = ARGPARSER.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(f"All tests passed in script: {pathlib.Path(__file__).name}")
if __name__ == "__main__":
main()
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
"""
import pathlib
from SETTINGS import ARGPARSER
from numpy.testing import assert_array_equal
import os
import tempfile
import pyplis
# ## SCRIPT FUNCTION DEFINITIONS
def create_ecII_cam_new_filters(cam_id: str) -> pyplis.Camera:
# Start with creating an empty Camera object
cam = pyplis.setupclasses.Camera(cam_id=cam_id, try_load_from_registry=False)
# 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 = cam_id
# 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 = None
# 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
def main():
# Create a test cam_info.txt file for this script in order not to mess with the camera info file shipped with pyplis
with tempfile.TemporaryDirectory() as temp_data_dir:
os.environ["PYPLIS_DATADIR"] = temp_data_dir
cam_info_file_tmp = pathlib.Path(temp_data_dir) / "cam_info.txt"
cam_info_file_tmp.touch()
available_cam_defs = pyplis.utils.get_cam_ids()
new_cam_id = "test_cam"
if new_cam_id in available_cam_defs:
raise ValueError(f"Camera with ID {new_cam_id} already exists, please choose a different ID")
cam = create_ecII_cam_new_filters(new_cam_id)
print(cam)
# 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(cam_info_file=cam_info_file_tmp)
cam_reload = pyplis.Camera(new_cam_id, cam_info_file=cam_info_file_tmp)
print(cam_reload)
# ## IMPORTANT STUFF FINISHED - everything below is of minor importance
# for educational purposes
options = ARGPARSER.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': new_cam_id,
'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(f"All tests passed in script: {pathlib.Path(__file__).name}")
if __name__ == "__main__":
main()
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).
"""
import pathlib
import pyplis
import matplotlib.pyplot as plt
from datetime import datetime
from SETTINGS import IMG_DIR, SAVEFIGS, SAVE_DIR, FORMAT, DPI, ARGPARSER
# ## RELEVANT DIRECTORIES AND PATHS
OFFSET_FILE = IMG_DIR / "EC2_1106307_1R02_2015091607064723_D0L_Etna.fts"
DARK_FILE = IMG_DIR / "EC2_1106307_1R02_2015091607064865_D1L_Etna.fts"
def main():
plt.close("all")
# ## Get all images in the image path which are FITS files (actually all)
all_paths = list(IMG_DIR.glob("*.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 = plt.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(SAVE_DIR / f"ex0_3_out_1.{FORMAT}",
format=FORMAT, dpi=DPI)
# Import script options
options = ARGPARSER.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.current_img().edit_log["darkcorr"],
sum(on_list.current_img().shape),
on_list.gaussian_blurring -
on_list.current_img().edit_log["blurring"]])
npt.assert_allclose([402.66284],
[off_img.mean() - on_img.mean()],
rtol=1e-7, atol=0)
print(f"All tests passed in script: {pathlib.Path(__file__).name}")
try:
if int(options.show) == 1:
plt.show()
except Exception:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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)
"""
import pathlib
import pyplis
from SETTINGS import IMG_DIR, ARGPARSER
from ex0_2_camera_setup import create_ecII_cam_new_filters
def main():
# create the camera object using the function defined in ex0_2_camera_setup.py
cam = create_ecII_cam_new_filters("test_cam")
# 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(f"list_id: {lst.list_id}, list_type: {lst.list_type}, number_of_files: {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(f"ImgLists linked to ImgList on: {on_list.linked_lists.keys()}")
print(f"Current file number on / off list: {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 = ARGPARSER.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.current_img().is_darkcorr +
off_list.current_img().is_darkcorr,
sum(on_list.current_img().shape),
on_list.gaussian_blurring -
on_list.current_img().edit_log["blurring"],
on_list.cfn])
npt.assert_allclose(actual=[on_list.get_dark_image().mean()],
desired=[190.56119],
rtol=1e-7)
print(f"All tests passed in script: {pathlib.Path(__file__).name}")
if __name__ == "__main__":
main()
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)
"""
import signal
import sys
import pyplis
flow = pyplis.plumespeed.OptflowFarneback()
def signal_handler(sig, frame):
print("\nExiting gracefully...")
sys.exit(0) # Exit without error
def main():
signal.signal(signal.SIGINT, signal_handler)
try:
flow.live_example()
except KeyboardInterrupt:
sys.exit(0)
if __name__ == "__main__":
main()
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 SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, ARGPARSER
import pathlib
import matplotlib.pyplot as plt
from pyplis import LineOnImage
from matplotlib.cm import get_cmap
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
def main():
plt.close("all")
fig, ax = plt.subplots(1, 2, figsize=(18, 9))
lines_r, lines_l = create_example_lines()
for k in range(len(lines_r)):
line = lines_r[k]
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]
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:
outfile = SAVE_DIR / f"ex0_6_out_1.{FORMAT}"
fig.savefig(outfile, format=FORMAT, dpi=DPI)
# Import script options
options= ARGPARSER.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(f"All tests passed in script: {pathlib.Path(__file__).name}")
try:
if int(options.show) == 1:
plt.show()
except Exception:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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.
"""
import pathlib
import pyplis
import matplotlib.pyplot as plt
from time import time
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, ARGPARSER, IMG_DIR
# 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"
def main():
plt.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=IMG_DIR / A53_ON,
cell_gas_cd=4.15e17,
cell_id="a53", filter_id="on")
cellcalib.set_cell_images(img_paths=IMG_DIR / A53_OFF,
cell_gas_cd=4.15e17,
cell_id="a53", filter_id="off")
cellcalib.set_cell_images(img_paths=IMG_DIR / A37_ON,
cell_gas_cd=8.59e17,
cell_id="a37", filter_id="on")
cellcalib.set_cell_images(img_paths=IMG_DIR / A37_OFF,
cell_gas_cd=8.59e17,
cell_id="a37", filter_id="off")
cellcalib.set_cell_images(img_paths=IMG_DIR / A57_ON,
cell_gas_cd=1.92e18,
cell_id="a57", filter_id="on")
cellcalib.set_cell_images(img_paths=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 = [
IMG_DIR / BG_BEFORE_ON,
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 = [
IMG_DIR / BG_BEFORE_OFF,
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(f"Time elapsed for preparing calibration data: {stop - start:.4f} s")
### IMPORTANT STUFF FINISHED
if SAVEFIGS:
outfile = SAVE_DIR / f"ex0_7_out_1.{FORMAT}"
ax.figure.savefig(outfile, format=FORMAT, dpi=DPI)
# Import script options
options = ARGPARSER.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]
desired = [1007480.56, 0.24401481, 3.194000052267778e+18, 4.779782684462987e+18, -2.7244424951657216e+16]
npt.assert_allclose(actual=actual, desired=desired, rtol=1e-5)
print(f"All tests passed in script: {pathlib.Path(__file__).name}")
try:
if int(options.show) == 1:
plt.show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
Example 0.8 - Parameterising optical flow histograms - The MultiGaussFit class
This script introduces the class pyplis.optimisation.MultiGaussFit which
is central for robust optical flow velocity retrievals.
Code
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.pyplot import rc_context
from pyplis import print_log as logger
from pyplis.model_functions import multi_gaussian_no_offset, multi_gaussian_same_offset
from pyplis.optimisation import MultiGaussFit
rc_context({'font.size': '12'})
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, ARGPARSER, IMG_DIR
def create_multigauss_testdata(x, params, offset, add_noise, noise_frac = 0.01):
"""Create test data set containing multiple Gaussians.
Parameters
----------
x : ndarray
index values for the signal
params : list
list of parameters for the Gaussians, each Gaussian is defined by
3 parameters: [amplitude, position, width]
offset : float
offset to be added to the signal
add_noise : bool
add noise to test data
noise_frac : float
determines noise amplitude (fraction relative to max amplitude of
Gaussian), defaults to 0.01 (1%). Only used if `add_noise` is True.
Returns
-------
tuple
2-element tuple containing
- the signal (y)
- the index (x)
"""
y = multi_gaussian_no_offset(x, *params) + offset
if add_noise:
y = y + max(y) * noise_frac * np.random.normal(0, 1, size=len(x))
return y, x
def create_noise_dataset():
"""Make pure noise dataset"""
x = np.linspace(0, 400, 401)
y = 5 * np.random.normal(0, 1, size=len(x))
return y, x
def plot_result(f: MultiGaussFit, mu_main_peak_analysis, sigma_main_peak_analysis, add_gaussians):
fig, axes = f.plot_result(True, figsize=(12, 10))
amp_main_peak = f.get_value(mu_main_peak_analysis)
for g in add_gaussians:
amp_tot = g[0] + f.offset
axes[0].annotate(
f"Additional peak:\nμ={int(g[1])}, σ={int(g[2])}",
xy=(g[1], amp_tot), # point to annotate
xytext=(-20, 0), # offset from point in points
textcoords='offset points',
bbox=dict(boxstyle='round,pad=0.2', fc='white', alpha=0.7),
ha='right', # horizontal alignment
va='top' if amp_tot > 300 else "bottom", # vertical alignment
zorder=10, # draw on top of other elements,
fontsize=8
)
axes[0].annotate(
f"Main peak:\nμ={int(mu_main_peak_analysis)}, σ={int(sigma_main_peak_analysis)}",
xy=(mu_main_peak_analysis, amp_main_peak), # point to annotate
xytext=(mu_main_peak_analysis, amp_main_peak*0.1), # offset from point in points
textcoords='data',
bbox=dict(boxstyle='round,pad=0.2', fc='white', alpha=0.7),
ha='center', # horizontal alignment
va='center', # vertical alignment
zorder=10, # draw on top of other elements,
fontsize=8
)
return fig
def main():
"""Example script to demonstrate the use of the MultiGaussFit class."""
figs_to_save = []
# Index (x axis) for the multigauss test data
signal_index = np.linspace(0, 400, 401)
# Case 1: Single Gaussian without offset
testdata = create_multigauss_testdata(
x=signal_index,
params=[100, 200, 10], # Amplitude, mu, sigma
offset=0,
add_noise=False)
# Create MultiGaussFit object and fit the data
# If data are provided and do_fit is True, the fit
# is performed automatically during initialisation
f = MultiGaussFit(data=testdata[0], index=testdata[1], do_fit=True)
# Get the number of Gaussians that were found in the data
n_gaussians = f.num_of_gaussians
logger.info(f"Number of Gaussians found: {n_gaussians}")
# you can also retrieve the parameters (amp, mu, sigma) of each Gaussian
for gaussian in f.gaussians():
logger.info(f"Gaussian parameters: (A={gaussian[0]:.2f}, "
f"μ={gaussian[1]:.2f}, σ={gaussian[2]:.2f})")
# For more complex data comprising a combination of multiple
# overlapping or distinct peaks, you can also run a post analysis
# on the fit result to identify the main peak (and its width) of the
# parameterised distribution, as well as additional peaks that may be present.
# See below for an example.
(mu_main_peak_analysis,
sigma_main_peak_analysis,
_,
add_gaussians) = f.analyse_fit_result(sigma_tol_overlaps=3)
# Visualise the result including the residual
fig = plot_result(f, mu_main_peak_analysis, sigma_main_peak_analysis, add_gaussians)
figs_to_save.append(fig)
# Case 2: 5 Gaussians with offset, some overlapping
gaussians = [
# Amplitude, mu, sigma (5 Gaussians)
150, 30, 8,
200, 110, 3,
300, 150, 20,
75, 370, 40,
300, 250, 1]
testdata = create_multigauss_testdata(
x=signal_index,
params=gaussians,
offset=45,
add_noise=True,
noise_frac=0.03
)
f = MultiGaussFit()
f.set_data(*testdata)
f.run_optimisation()
# For illustration: compute the 1st and 2nd moment of the parameterised
# multi Gaussian distribution, i.e. the most distinctive peak. Note that
# the 1st moment is the mean and the 2nd moment is the variance.
# Note that this considers all detected Gaussians, not just the main peak (see fit result)
# Regardless of this, the result is still a good approximation of the main peak as can
# be seen in the console output at the end of this script.
p_norm = f.normalise_params()
x = f.index
data_norm = multi_gaussian_no_offset(x, *p_norm)
mu_all_peaks = f.det_moment(x, data_norm, 0, 1)
sigma_all_peaks = np.sqrt(f.det_moment(x, data_norm, mu_all_peaks, 2))
(mu_main_peak_analysis,
sigma_main_peak_analysis,
_,
add_gaussians) = f.analyse_fit_result(sigma_tol_overlaps=3)
# Visualise the result including the residual
fig = plot_result(f, mu_main_peak_analysis, sigma_main_peak_analysis, add_gaussians)
figs_to_save.append(fig)
logger.info(f"\nMu, sigma (all Gaussians): {mu_all_peaks:.2f}, {sigma_all_peaks:.2f}")
logger.info(f"\nMu, sigma (Gaussian main peak): {mu_main_peak_analysis:.2f}, {sigma_main_peak_analysis:.2f}")
### IMPORTANT STUFF FINISHED
if SAVEFIGS:
for i, fig in enumerate(figs_to_save):
outfile = SAVE_DIR / f"ex0_8_out_{i+1}.png"
fig.savefig(outfile, format=FORMAT, dpi=DPI)
# Import script options
options = ARGPARSER.parse_args()
try:
if int(options.show) == 1:
plt.show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
plt.close("all")
main()
Advanced examples for emission rate analysis
The following scripts demonstrate a complete workflow for emission rate analysis using pyplis. Each example builds upon the previous ones, culminating in a full emission rate analysis in Example 12 - Emission rate analysis (Etna example data). These scripts cover essential aspects including:
Setting up the analysis environment
Handling measurement geometry
Processing plume background data
Performing calibrations (cell-based and DOAS)
Retrieving plume velocities
Correcting for signal dilution
Computing final emission rates
The examples use data from Mt. Etna (Italy) to demonstrate real-world applications of the pyplis software.
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 SETTINGS import IMG_DIR, ARGPARSER
import pathlib
import pyplis as pyplis
from datetime import datetime
import matplotlib.pyplot as plt
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:
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}
# Create BaseSetup object (which creates the MeasGeometry object)
stp = pyplis.setupclasses.MeasSetup(
base_dir=IMG_DIR,
start=start,
stop=stop,
camera=cam,
source=source,
wind_info=wind_info
)
# Create analysis object (from BaseSetup)
# The dataset takes care of finding all vali
return pyplis.Dataset(stp)
# SCRIPT MAIN FUNCTION
def main():
plt.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
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 = ARGPARSER.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
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.current_img().is_darkcorr +
off_list.current_img().is_darkcorr,
sum(on_list.current_img().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(f"All tests passed in script: {pathlib.Path(__file__).name}")
try:
if int(options.show) == 1:
plt.show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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 geonum import GeoPoint
import pathlib
import matplotlib.pyplot as plt
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, ARGPARSER
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex01_analysis_setup import create_dataset
# 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", auto_topo_access=True)
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: "
f"elev = {new_elev:.1f}, azim = {new_azim:.1f}")
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 = plt.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
def main():
plt.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(dict(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"
f"Previous: {plume_dist_cfov:.3f}m\n"
f"New: {plume_dist_cfov_new:.3f}m")
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
outfile = SAVE_DIR / f"ex02_out_1.{FORMAT}"
map_.ax.figure.savefig(outfile,format=FORMAT,dpi=DPI)
outfile = SAVE_DIR / f"ex02_out_2.{FORMAT}"
fig.savefig(outfile,format=FORMAT,dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
options = ARGPARSER.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
# 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.058292, 37.753504, 0.182902, 0.145056],
rtol=1e-5)
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
]
# check some basic properties / values of the geometry
npt.assert_allclose(actual=actual,desired=desired,rtol=1e-5)
print(f"All tests passed in script: {pathlib.Path(__file__).name}")
try:
if int(options.show) == 1:
plt.show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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 typing import Tuple
from matplotlib.figure import Figure
import numpy as np
import pathlib
import matplotlib.pyplot as plt
import pyplis
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, ARGPARSER
# 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 = IMG_DIR / 'EC2_1106307_1R02_2015091607065477_F01_Etna.fts'
BG_FILE = IMG_DIR /'EC2_1106307_1R02_2015091607022602_F01_Etna.fts'
OFFSET_FILE = IMG_DIR / 'EC2_1106307_1R02_2015091607064723_D0L_Etna.fts'
DARK_FILE = IMG_DIR / 'EC2_1106307_1R02_2015091607064865_D1L_Etna.fts'
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: pyplis.PlumeBackgroundModel,
plume_img: pyplis.Img) -> Tuple[dict, Figure]:
"""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_img)
current_params = bg_model.settings_dict()
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
axes[0].set_title("Manually set parameters")
pyplis.plumebackground.plot_sky_reference_areas(plume_img, current_params,
ax=axes[0])
pyplis.plumebackground.plot_sky_reference_areas(plume_img, 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 = plt.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
def main():
plt.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=bg_model, plume_img=plume)
# Script option
if USE_AUTO_SETTINGS:
bg_model.update(**auto_params)
# Model 4 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:
outfile = SAVE_DIR / f"ex03_out_1.{FORMAT}"
fig0.savefig(outfile, format=FORMAT,
dpi=DPI, transparent=True)
for k in range(len(_tau_figs)):
outfile = SAVE_DIR / f"ex03_out_{k+2}.{FORMAT}"
_tau_figs[k].savefig(outfile, format=FORMAT, dpi=DPI)
outfile = outfile = SAVE_DIR / f"ex03_out_{k+3}.{FORMAT}"
fig6.savefig(outfile, format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
options = ARGPARSER.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
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(f"All tests passed in script: {pathlib.Path(__file__).name}")
try:
if int(options.show) == 1:
plt.show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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
"""
import pyplis
import pathlib
import matplotlib.pyplot as plt
from time import time
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, ARGPARSER, PCS1
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex01_analysis_setup import create_dataset
from ex02_meas_geometry import find_viewing_direction
# 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 = IMG_DIR / 'EC2_1106307_1R02_2015091607022602_F01_Etna.fts'
path_bg_off = IMG_DIR / 'EC2_1106307_1R02_2015091607022820_F02_Etna.fts'
# Get on and off lists and activate dark correction
on_list = dataset.get_list("on")
off_list = dataset.get_list("off")
# Deactivate automatic reload in list while changing some list
# attributes
on_list.auto_reload = False
off_list.auto_reload = False
on_list.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_list.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_list.set_bg_img(bg_on)
off_list.set_bg_img(bg_off)
# automatically set gas free areas
on_list.bg_model.set_missing_ref_areas(on_list.current_img())
# Now update some of the information from the automatically set sky ref
# areas
on_list.bg_model.xgrad_line_startcol = 20
on_list.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_list.bg_model.mode = bg_corr_mode
off_list.bg_model.mode = bg_corr_mode
on_list.aa_mode = True # activate AA mode
off_list.auto_reload = True
on_list.auto_reload = True
on_list.meas_geometry = geom
return on_list
def main():
plt.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.current_img().shape), aa_list.current_img().mean()
ax = aa_list.show_current(zlabel=r"$\tau_{AA}$")
print(f"Elapsed time: {time() - t0} s")
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.current_img()
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:
outfile = SAVE_DIR / f"ex04_out_1.{FORMAT}"
ax.figure.savefig(outfile, format=FORMAT, dpi=DPI)
outfile = SAVE_DIR / f"ex04_out_2.{FORMAT}"
ax1.figure.savefig(outfile, format=FORMAT, dpi=DPI)
outfile = SAVE_DIR / f"ex04_out_3.{FORMAT}"
fig.savefig(outfile, format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
options = ARGPARSER.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
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.current_img().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(f"All tests passed in script: {pathlib.Path(__file__).name}")
try:
if int(options.show) == 1:
plt.show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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).
"""
import pathlib
import numpy as np
import matplotlib.pyplot as plt
import pyplis
from datetime import datetime
from time import time
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, ARGPARSER
# File path for storing cell AA calibration data including sensitivity
# correction mask
AA_CALIB_FILE = 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
def main():
plt.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])
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)
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
ax0.figure.savefig(SAVE_DIR / f"ex05_2_out_1.{FORMAT}", format=FORMAT, dpi=DPI)
ax1.figure.savefig(SAVE_DIR / f"ex05_2_out_2.{FORMAT}", format=FORMAT, dpi=DPI)
ax2.figure.savefig(SAVE_DIR / f"ex05_2_out_3.{FORMAT}", format=FORMAT, dpi=DPI)
# 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(f"Time elapsed for preparing calibration data: {stop - start:.4f} s")
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
options = ARGPARSER.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
calib_reload = pyplis.CellCalibData()
calib_reload.load_from_fits(AA_CALIB_FILE)
calib_reload.fit_calib_data(polyorder=2, through_origin=True)
timestamps = np.array([datetime(2015, 9, 16, 7, 0, 25, 127400), datetime(2015, 9, 16, 7, 0, 58, 637400), datetime(2015, 9, 16, 7, 1, 32, 647400)])
npt.assert_array_equal(timestamps, calib_reload.time_stamps)
# 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)]
desired=[8.59e+17,4.15e+17,1.924e+18,-1.831327e+16,0.0]
npt.assert_allclose(actual=vals_approx,desired=desired,rtol=1e-5)
# 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()]
desired=[1.892681e+18, -3.742539e+19, 2.153654e+18, 2.119768e+18]
npt.assert_allclose(actual=actual,desired=desired,rtol=1e-5)
print(f"All tests passed in script: {pathlib.Path(__file__).name}")
try:
if int(options.show) == 1:
plt.show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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 pathlib import Path
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 pyplis.imagelists import ImgList
from pyplis.processing import ImgStack
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, ARGPARSER
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list
# RELEVANT DIRECTORIES AND PATHS
# Directory containing DOAS result files
DOAS_DATA_DIR = IMG_DIR / ".." / "spectra" / "plume_prep" / "min10Scans" / "ResultFiles"
# Output file path for storing the AA image stack
AA_IMG_STACK_PATH = 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: ImgList,
output_path: Path,
roi_abs=None,
pyrlevel=None) -> ImgStack:
"""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()
stack.save_as_fits(save_dir=output_path.parent, save_name=output_path.name)
return stack
# 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.fov.sigma_x_abs,
calib.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)
def main():
# close all plots
close("all")
# import DOAS SO2 column density timeseries
doas_time_series = load_doas_results()
# get the AA image list
aa_list = prepare_aa_image_list()
# Calculate AA image stack from the list of images in the list
stack = make_aa_stack_from_list(
aa_list=aa_list,
output_path=AA_IMG_STACK_PATH,
pyrlevel=2
)
# Instantiate the engine for retrieving the FOV and calibration data
fov_engine = pyplis.doascalib.DoasFOVEngine(stack, doas_time_series)
calib_pears = fov_engine.perform_fov_search(method="pearson", mergeopt="nearest")
calib_ifr = fov_engine.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()
_, 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]
#Perform FOV search within ROI around result from pearson fov
#search at full resolution (pyrlevel=0)
aa_list = prepare_aa_image_list()
num_merge, h, w = fov_engine.img_stack.shape
s_fine = fov_engine.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("")
outfile = SAVE_DIR / f"ex06_out_{k + 1}.{FORMAT}"
ax.figure.savefig(outfile, 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
# Import script options
options = ARGPARSER.parse_args()
if int(options.test):
# test re-loading of stack from FITS file
stack = pyplis.processing.ImgStack()
stack.load_stack_fits(AA_IMG_STACK_PATH)
num, h, w = stack.shape
num2 = fov_engine.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(f"All tests passed in script: {Path(__file__).name}")
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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.
"""
import pyplis
from pathlib import Path
import numpy as np
from matplotlib.pyplot import close, subplots, show
from matplotlib.patches import Circle
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, ARGPARSER
# 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
CELL_AA_CALIB_FILE = SAVE_DIR / "ex05_cellcalib_aa.fts"
# RELEVANT DIRECTORIES AND PATHS
# fits file containing DOAS calibration information (from ex6)
DOAS_CALIB_FILE = 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(f"Line {pcs1.line_id}")
axes[1].set_title(f"Line {pcs2.line_id}")
phi_p1_init = sum(p10) / num
phi_p2_init = sum(p20) / num
axes[0].plot(p10, "-", label=f"Init φ={phi_p1_init:.3f}")
axes[1].plot(p20, "-", label=f"Init φ={phi_p2_init:.3f}")
# store main results for testing in main()
results_p1 = [(np.nan, float(phi_p1_init))]
results_p2 = [(np.nan, float(phi_p2_init))]
for cd, aa_corr in aa_imgs_corr.items():
p1 = pcs1.get_line_profile(aa_corr.img)
p2 = pcs2.get_line_profile(aa_corr.img)
phi_p1 = sum(p1) / num
phi_p2 = sum(p2) / num
results_p1.append((float(cd), float(phi_p1)))
results_p2.append((float(cd), float(phi_p2)))
axes[0].plot(p1, "-", label=f"Cell CD: {cd:.2e} φ={phi_p1:.3f}")
axes[1].plot(p2, "-", label=f"Cell CD: {cd:.2e} φ={phi_p2:.3f}")
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, results_p1, results_p2
# SCRIPT MAIN FUNCTION
def main():
close("all")
if not DOAS_CALIB_FILE.exists():
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 calib")
ax0 = cellcalib.calib_data["aa"].plot(ax=ax0, c="r", add_label_str="Cell calib")
ax0.set_title("")
ax0.set_xlim([0, 0.5])
# Get current AA image from image list
aa_init_img = 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.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, _, results_p1, results_p2 = plot_pcs_comparison(aa_init_img, aa_imgs_corr, pcs1, pcs2)
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
ax0.figure.savefig(SAVE_DIR / f"ex07_out_1.{FORMAT}", format=FORMAT, dpi=DPI)
ax.figure.savefig(SAVE_DIR / f"ex07_out_2.{FORMAT}", format=FORMAT, dpi=DPI)
fig.savefig(SAVE_DIR / f"ex07_out_3.{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 = ARGPARSER.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
assert len(doascalib.tau_vec) == 88
assert len(doascalib.cd_vec) == 88
npt.assert_allclose(
actual=[
doascalib.tau_vec.mean(),
doascalib.cd_vec.mean()
],
desired=[
0.10984407750550997,
1.2129406193969774e+18
],
rtol=1e-7
)
npt.assert_allclose(
actual=[
cellcalib.calib_data["aa"].cd_vec,
cellcalib.calib_data["aa"].tau_vec,
cellcalib.calib_data["on"].tau_vec,
cellcalib.calib_data["off"].tau_vec
],
desired=[
[4.150e+17, 8.590e+17, 1.924e+18],
[0.11438841, 0.20960312, 0.45484966],
[0.25561193, 0.34218314, 0.58974224],
[0.14122352, 0.13258003, 0.13489263]
],
rtol=1e-7
)
npt.assert_allclose(
actual=results_p1,
desired=[
(np.nan, 0.06704542158954685),
(4.15e+17, 0.06718042466335893),
(8.59e+17, 0.0672710026562296),
(1.924e+18, 0.06713881262426598)],
rtol=1e-7
)
npt.assert_allclose(
actual=results_p2,
desired=[
(np.nan, 0.09656624766757335),
(4.15e+17, 0.08434758808651711),
(8.59e+17, 0.08673337715111279),
(1.924e+18, 0.08756188495164571)],
rtol=1e-7
)
print(f"All tests passed in script: {Path(__file__).name}")
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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 matplotlib.pyplot import close, show, subplots
from time import time
from pathlib import Path
from pyplis.plumespeed import VeloCrossCorrEngine
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, ARGPARSER, LINES
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list
# SCRIPT OPTONS
# distance in pixels between two lines used for cross correlation analysis
OFFSET_PIXNUM = 40
# 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
def 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(pcs=PCS, imglist=aa_list)
cc.create_parallel_pcs_offset(offset_pix=40,
color=COLOR_PCS_OFFS,
linestyle="--")
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)
cc.load_pcs_profile_img(SAVE_DIR / PCS_PROFILES_PIC_NAME, line_id="pcs")
cc.load_pcs_profile_img(SAVE_DIR / OFFSET_PROFILES_PIC_NAME, line_id="pcs_offset")
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()
_, 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(f"Result performance analysis\n"
f"Number of images: {STOP_IDX - START_IDX}\n"
f"Create ICA images: {t1 - t0:.3f} s\n"
f"Cross-corr analysis: {t2 - t1:.3f} s")
print(f"Retrieved plume velocity of v = {velo:.2f} m/s")
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
for k in range(len(axes)):
outfile= SAVE_DIR / f"ex08_out_{k + 1}.{FORMAT}"
axes[k].figure.savefig(outfile, format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
options = ARGPARSER.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
results = cc.results
npt.assert_array_equal([len(results),
len(results["coeffs"]),
len(results["ica_tseries_offset"]),
len(results["ica_tseries_shift"])],
[6, 1578, 7891, 7712])
npt.assert_allclose(actual=[results["velo"], results["lag"]],
desired=[4.434, 17.9],
rtol=1e-2)
print(f"All tests passed in script: {Path(__file__).name}")
try:
if int(options.show) == 1:
show()
except Exception:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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."""
import pyplis
from pathlib import Path
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, ARGPARSER, 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
rcParams["font.size"] = 16
def analyse_and_plot(lst, lines):
fig = figure(figsize=(14, 8))
ax0 = fig.add_axes([0.01, 0.15, 0.59, 0.8])
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(0.05)
fl = lst.optflow
fl.plot(ax=ax0)
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)
_, 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)
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
def 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(pyrlevel=1)
# set the pyramid level in the list
aa_list.pyrlevel = 1
# 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 = 2
s.roi_rad_abs = [0, 0, 1344, 730]
aa_list.optflow_mode = True
plume_mask = pyplis.Img(aa_list.get_thresh_mask(0.05))
plume_mask.show(tit="AA threshold mask")
figs.append(analyse_and_plot(aa_list, LINES))
PCS1, PCS2 = 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)
aa_list.goto_img(0)
stop_idx = aa_list.nof - 1
print("Computing local optical flow properties for 2 PCS lines for each image in the image list")
for k in range(0, stop_idx):
if k%20==0:
print(f"Computation ongoing ({k}/{stop_idx})")
plume_mask = aa_list.get_thresh_mask(0.05)
plume_props_l1.get_and_append_from_farneback(
fl, line=PCS1, pix_mask=plume_mask,
dir_multi_gauss=True)
plume_props_l2.get_and_append_from_farneback(
fl, line=PCS2, pix_mask=plume_mask,
dir_multi_gauss=True)
aa_list.goto_next()
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([])
figs.append(fig)
# Save the time series as txt
plume_props_l1.save_txt(SAVE_DIR / "ex09_plumeprops_young_plume.txt")
plume_props_l2.save_txt(SAVE_DIR / "ex09_plumeprops_aged_plume.txt")
if SAVEFIGS:
for k in range(len(figs)):
outfile = SAVE_DIR / f"ex09_out_{k + 1}.{FORMAT}"
figs[k].savefig(outfile, format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
options = ARGPARSER.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
import numpy as np
res = plume_props_l1
assert res.roi_id == "young_plume"
assert np.isnan(res._len_mu_norm).sum() == 5
assert np.isnan(res._len_sigma_norm).sum() == 5
assert np.isnan(res._dir_mu).sum() == 5
assert (np.array(res._dir_sigma) == 180).sum() == 5
assert (np.array(res._significance) == 0).sum() == 5
npt.assert_allclose(
actual=[
np.nanmean(res._len_mu_norm),
np.nanmean(res._len_sigma_norm),
np.nanmean(res._dir_mu),
np.nanmean(res._dir_sigma),
np.nanmean(res._del_t),
np.nanmean(res._significance),
],
desired=[
0.8446354828438595,
0.288896828513025,
-56.972104109591605,
19.873233016577196,
4.206971153846154,
0.7650539016892137
],
rtol=1e-5)
res = plume_props_l2
assert res.roi_id == "old_plume"
assert np.isnan(res._len_mu_norm).sum() == 5
assert np.isnan(res._len_sigma_norm).sum() == 5
assert np.isnan(res._dir_mu).sum() == 5
assert (np.array(res._dir_sigma) == 180).sum() == 5
assert (np.array(res._significance) == 0).sum() == 5
npt.assert_allclose(
actual=[
np.nanmean(res._len_mu_norm),
np.nanmean(res._len_sigma_norm),
np.nanmean(res._dir_mu),
np.nanmean(res._dir_sigma),
np.nanmean(res._del_t),
np.nanmean(res._significance),
],
desired=[
1.0966500134369748,
0.39550772843047144,
-93.31145142654839,
20.139529500869035,
4.206971153846154,
0.6430045694337098
],
rtol=1e-5)
print(f"All tests passed in script: {Path(__file__).name}")
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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."""
import pyplis
from datetime import datetime
from matplotlib.pyplot import show
from pathlib import Path
# IMPORT GLOBAL SETTINGS
from SETTINGS import IMG_DIR, ARGPARSER
def get_bg_image_lists():
# start time of sky background image acquisition
start = datetime(2015, 9, 16, 7, 2, 5)
# stop time of sky background image acquisition
stop = datetime(2015, 9, 16, 7, 2, 30)
# Define camera (in the example dataset from Etna, the ECII camera is used)
cam_id = "ecII"
# Declare the on and offband camera filters used
filters = [pyplis.utils.Filter(type="on", acronym="F01"),
pyplis.utils.Filter(type="off", acronym="F02")]
# Create Camera instance using the ECII camera type and filters
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_list, off_list = ds.get_list("on"), ds.get_list("off")
on_list.darkcorr_mode = True
off_list.darkcorr_mode = True
return on_list, off_list
def main():
on_list, off_list = get_bg_image_lists()
on_list.show_current()
off_list.show_current()
# Import script options
options = ARGPARSER.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_allclose(actual=[on_list.current_img().mean(), off_list.current_img().mean()],
desired=[2555.4597, 2826.8848],
rtol=1e-3)
print(f"All tests passed in script: {Path(__file__).name}")
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
# Execute main function
main()
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)
"""
import pyplis as pyplis
from geonum import GeoPoint
from matplotlib.pyplot import show, close, subplots, Rectangle, plot
from datetime import datetime
from pathlib import Path
# IMPORT GLOBAL SETTINGS
from SETTINGS import IMG_DIR, SAVEFIGS, SAVE_DIR, FORMAT, DPI, ARGPARSER
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex10_bg_imglists import get_bg_image_lists
from pyplis.dataset import Dataset
from pyplis.geometry import MeasGeometry
from pyplis.imagelists import ImgList
# SCRIPT OPTONS
# lower boundary for I0 value in dilution fit
I0_MIN = 0.0
AA_THRESH = 0.03
# 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")
# Create lines covering topographic terrain features in the images used
# 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="-")
# Lines covering topographic terrain features used for the signal dilution analysis
LOCAL_TOPOGRAPHY_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, from example script 8)
PLUME_VELO = 4.14 # m/s
# FITS file containing DOAS calibration data (from example script 6)
CALIB_FILE = SAVE_DIR / "ex06_doascalib_aa.fts"
def create_dataset_dilution():
"""Create a :class:`pyplis.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, 0)
stop = datetime(2015, 9, 16, 6, 47, 0)
# 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": 0} # 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(
base_dir=IMG_DIR,
start=start,
stop=stop,
camera=cam,
source=source,
wind_info=wind_info
)
return pyplis.dataset.Dataset(stp)
def find_view_dir(geom: MeasGeometry, plot: bool) -> MeasGeometry:
"""Perform a correction of the viewing direction using crater in img.
Args:
geom: measurement geometry
plot: if True, the result is plotted
Returns:
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=plot)
return geom
def prepare_lists(dataset: Dataset) -> tuple[ImgList, ImgList]:
"""Prepare on and off lists for dilution analysis.
Steps:
1. get on and offband list
2. load background image list on and off (from ex10)
3. set image preparation and assign background images to on / off list
4. configure plume background model settings
Args:
- dataset: the dilution dataset (see :func:`create_dataset_dilution`)
Returns:
- onlist: ImgList object for onband filter
- offlist: ImgList object for offband filter
"""
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, dil):
"""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 LOCAL_TOPOGRAPHY_LINES:
line.plot_line_on_grid(ax=ax, marker="", color=line.color,
lw=2, ls=line.linestyle)
for x, y, _ in dil._add_points:
plot(x, y, " or")
return ax
def main():
if not CALIB_FILE.exists():
raise IOError("Calibration file could not be found at specified "
"location:\n %s\nPlease run example 6 first")
close("all")
pcs_line = PCS_LINE
calib = pyplis.doascalib.DoasCalibData()
calib.load_from_fits(CALIB_FILE)
# create dataset and correct viewing direction
ds = create_dataset_dilution()
geom = find_view_dir(ds.meas_geometry, plot=True)
# get plume distance image
pix_dists, _, _ = geom.compute_all_integration_step_lengths()
# Create dilution correction class
dil = pyplis.dilutioncorr.DilutionCorr(LOCAL_TOPOGRAPHY_LINES, geom, skip_pix=SKIP_PIX_LINES)
# you may 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 via input arg `dist`,
# since in some regions, the topographic data used is incomplete or
# has too low resolution). Here, `dist` is not provided, which means that the
# distance will be determined based on the measurement
# geometry (cam position, viewing direction, and optics) and the local topography.
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,_ = 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)
fit_result_on = dil.fit(
img=on_vigncorr,
rad_ambient=ia_on,
i0_min=I0_MIN)
ax0 = dil.plot_fit_result(
dists=fit_result_on.dists,
rads=fit_result_on.rads,
rad_ambient=ia_on,
i0=fit_result_on.i0,
ext=fit_result_on.ext)
ax0.set_ylabel("Terrain radiances (on band)", fontsize=14)
ax0.set_ylim([0, 2500])
# perform dilution anlysis and retrieve extinction coefficients (off-band)
fit_result_off = dil.fit(
img=off_vigncorr,
rad_ambient=ia_off,
i0_min=I0_MIN)
ax1 = dil.plot_fit_result(
dists=fit_result_off.dists,
rads=fit_result_off.rads,
rad_ambient=ia_off,
i0=fit_result_off.i0,
ext=fit_result_off.ext)
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
ext_on = fit_result_on.ext
ext_off = fit_result_off.ext
# assign the just retrieved extinction coefficients together with the
# acq. time of images they were retrieved from
onlist.set_ext_coeffs(values=[ext_on], timestamps=[on_vigncorr.start_acq])
offlist.set_ext_coeffs(values=[ext_off], timestamps=[off_vigncorr.start_acq])
# 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(SAVE_DIR / "ex11_ext_scat_on.csv", header=False)
offlist.ext_coeffs.to_csv(SAVE_DIR / "ex11_ext_scat_off.csv", header=False)
# Activate automatic dilution correction in both lists
# get dilution corrected on and off-band image
onlist.dilcorr_mode = True
offlist.dilcorr_mode = True
# 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.current_img()
tau_off_corr = offlist.current_img()
# determine corrected SO2-CD image from the image lists
# by applying the calibration function to the AA image
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
# Get uncorrected SO2-CD image
so2_img_uncorr = calib(onlist.current_img() - offlist.current_img())
# 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, dil)
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=f"Uncorr: $\\Phi_{{SO2}}=$ {phi_uncorr / 1000.0:.2f} (+/- {phi_uncorr_err / 1000.0:.2f}) kg/s")
ax3.plot(so2_cds_corr, "-g", lw=3,
label=f"Corr: $\\Phi_{{SO2}}=$ {phi_corr / 1000.0:.2f} (+/- {phi_corr_err / 1000.0:.2f}) kg/s")
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(SAVE_DIR / f"ex11_out_{k}.{FORMAT}",
format=FORMAT, dpi=DPI)
basemap.ax.set_axis_off()
basemap.ax.view_init(15, 345)
basemap.ax.figure.savefig(SAVE_DIR / f"ex11_out_5.{FORMAT}",
format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
options = ARGPARSER.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
# check correction of viewing direction in geometry
npt.assert_allclose(
actual=[geom._cam["elev"], geom._cam["azim"], geom._cam["elev_err"], geom._cam["azim_err"]],
desired=[14.212510910775158, 280.3278110428114, 1.0656026217149126, 1.0616956101542883],
rtol=1e-5)
# Check inputs to the dilution correction fit engine
npt.assert_allclose(
actual=[on_vigncorr.mean(), off_vigncorr.mean(), ia_on, ia_off],
desired=[2849.1296, 2576.153, 3963.2626953125, 3525.40478515625],
rtol=1e-7)
npt.assert_allclose(
actual=[
fit_result_on.dists.mean(), fit_result_on.rads.mean(),
fit_result_off.dists.mean(), fit_result_off.rads.mean()
],
desired=[6854.925234914583, 1794.5193, 6854.925234914583, 1378.1238],
rtol=1e-7)
# check the retrieved extinction coefficients and i0 values
# for on and off band images
npt.assert_allclose(
actual=[
fit_result_on.ext, fit_result_on.i0,
fit_result_off.ext, fit_result_off.i0,
phi_corr, phi_corr_err,
phi_uncorr, phi_uncorr_err
],
desired=[
7.41e-5,449.9,
6.55e-5, 236.5,
9.67e3, 1.99e3,
3.35e3, 0.81e3],
rtol=1e-2)
print(f"All tests passed in script: {Path(__file__).name}")
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()
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
instances of `EmissionRates` class.
"""
from pathlib import Path
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, ARGPARSER, LINES
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list
rc_context({'font.size': '18'})
PCS = LINES[0]
# If false, then only the working environment is initalised
DO_EVAL = False
# 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 = SAVE_DIR / "ex06_doascalib_aa.fts"
# Scattering extinction coeffcients from example script 11 (stored as txt)
EXT_ON_CACHE_FILE = SAVE_DIR / "ex11_ext_scat_on.csv"
EXT_OFF_CACHE_FILE = SAVE_DIR / "ex11_ext_scat_off.csv"
# AA sensitivity correction mask retrieved from cell calib in script 7
CORR_MASK_FILE = 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_CACHE_FILE = 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(SAVE_DIR / "ex12_flux_velo_glob.txt")
res1.save_txt(SAVE_DIR / "ex12_flux_flow_raw.txt")
res2.save_txt(SAVE_DIR / "ex12_flux_flow_histo.txt")
res3.save_txt(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
def main():
close("all")
figs = []
if not CALIB_FILE.exists():
raise FileNotFoundError(
f"Calibration file {CALIB_FILE} could not be found at specified "
f"Please run example 6 first")
if not CORR_MASK_FILE.exists():
raise FileNotFoundError(
f"Cannot find AA correction mask at {CORR_MASK_FILE}, please run example"
f"script 7 first")
# convert the retrieval line to the specified pyramid level (script option)
pcs = PCS.convert(to_pyrlevel=PYRLEVEL)
# Load results of optical flow histogram analysis performed for
# PCS 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
p = pyplis.LocalPlumeProperties()
p.load_txt(RESULT_PLUMEPROPS_HISTO_CACHE_FILE)
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
# 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_CACHE_FILE)
aa_list.get_off_list().import_ext_coeffs_csv(EXT_OFF_CACHE_FILE)
# 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 ex07)
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=f"Dilution corr: {DILCORR}")
ax = ana.plot_bg_roi_rect(ax=ax, to_pyrlevel=PYRLEVEL)
figs.append(ax.figure)
ana.run_retrieval(start_index=START_INDEX,
stop_index=STOP_INDEX)
figs.append(plot_and_save_results(ana))
# the EmissionRates 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(SAVE_DIR / f"ex12_out_{k + 1}.{FORMAT}",
format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
options = ARGPARSER.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
assert len(ana.results) == 1
assert "young_plume" in ana.results
assert len(ana.results["young_plume"]) == 4
assert list(ana.results["young_plume"]) == ["glob", "flow_raw", "flow_histo", "flow_hybrid"]
res = ana.results["young_plume"]["glob"]
assert isinstance(res, pyplis.EmissionRates)
npt.assert_allclose(actual=[res.phi.mean(), res.phi_err.mean(), res.velo_eff.mean(), res.velo_eff_err.mean()],
desired=[7202.8308, 2551.5447, PLUME_VELO_GLOB, PLUME_VELO_GLOB_ERR],
rtol=1e-3)
res = ana.results["young_plume"]["flow_raw"]
assert isinstance(res, pyplis.EmissionRates)
npt.assert_allclose(actual=[res.phi.mean(), res.phi_err.mean(), res.velo_eff.mean(), res.velo_eff_err.mean()],
desired=[4684.6984, 749.426, 2.7767, 0.4165],
rtol=1e-3)
res = ana.results["young_plume"]["flow_histo"]
assert isinstance(res, pyplis.EmissionRates)
npt.assert_allclose(actual=[res.phi.mean(), res.phi_err.mean(), res.velo_eff.mean(), res.velo_eff_err.mean()],
desired=[5636.7601, 996.3797, 3.3551, 0.5457],
rtol=1e-3)
res = ana.results["young_plume"]["flow_hybrid"]
assert isinstance(res, pyplis.EmissionRates)
npt.assert_allclose(actual=[res.phi.mean(), res.phi_err.mean(), res.velo_eff.mean(), res.velo_eff_err.mean()],
desired=[5853.8232, 954.3205, 3.4865, 0.5288],
rtol=1e-3)
print(f"All tests passed in script: {Path(__file__).name}")
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
if __name__ == "__main__":
main()