Example scripts¶
pyplis example scripts. The scripts require downloading the pyplis test data.
Note
The scripts are based on the latest commit in the GitHub repo. If you have installed an older version of pyplis, please use the corresponding scripts which are provided here.
Introductory scripts¶
These scripts give an introduction into basic features and classes of pyplis.
Example 0.1 - Image representation¶
Introduction into Img
object and basic usage including correction for dark current and detector offset.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script 1 - Image import.
This script loads an exemplary image which is part of the pyplis supplied
data. Image data in pyplis is represented by the ``Img`` class, which also
allows for storage of image meta data and keeps track of changes applied to
the image data (e.g. cropping, blurring, dark correction).
The script also illustrates how to manually work with image meta data
encoded in the image file name. The latter can be performed automatically in
pyplis using file name conventions (which can be specified globally, see next
script).
"""
from __future__ import (absolute_import, division)
# Check script version
from SETTINGS import check_version
# imports from SETTINGS.py
from SETTINGS import OPTPARSE, SAVE_DIR
from os.path import join, basename
from datetime import datetime
from matplotlib.pyplot import close
import pyplis
check_version()
# file name of test image stored in data folder
IMG_FILE_NAME = "test_201509160708_F01_335.fts"
IMG_DIR = join(pyplis._LIBDIR, "data")
if __name__ == "__main__":
close("all")
img_path = join(IMG_DIR, IMG_FILE_NAME)
# Create Img object (Img objects can be initiated both with image file
# paths but also with data in memory in form of a 2D numpy array)
img = pyplis.image.Img(img_path)
# log mean of uncropped image for testing mode
avg = img.mean()
# The file name of the image includes some image meta information which can
# be set manually (this is normally done automatically by defining a file
# name convention, see next script)
# split filename using delimiter "_"
spl = IMG_FILE_NAME.split(".")[0].split("_")
# extract acquisition time and convert to datetime
acq_time = datetime.strptime(spl[1], "%Y%m%d%H%M")
# extract and convert exposure time from filename (convert from ms -> s)
texp = float(spl[3]) / 1000
img.meta["start_acq"] = acq_time
img.meta["texp"] = texp
img.meta["f_num"] = 2.8
img.meta["focal_length"] = 25e-3
# add some blurring to the image
img.add_gaussian_blurring(sigma_final=3)
# crop the image edges
roi_crop = [100, 100, 1244, 924] # [x0, y0, x1, y1]
img.crop(roi_abs=roi_crop)
# apply down scaling (gauss pyramid)
img.to_pyrlevel(2)
# ## Show image
img.show()
img.save_as_fits(SAVE_DIR, "ex0_1_imgsave_test")
img_reload = pyplis.Img(join(SAVE_DIR, "ex0_1_imgsave_test.fts"))
# print image information
print(img)
# ## IMPORTANT STUFF FINISHED - everything below is of minor importance
# for educational purposes
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
vals = [img.mean(), avg, int(spl[1]), texp,
img_reload.meta["texp"],
img_reload.meta["f_num"],
img_reload.meta["focal_length"]]
npt.assert_almost_equal([2526.4623422672885,
2413.0870433989026,
201509160708,
0.335,
0.335,
2.8,
0.025],
vals, 4)
print("All tests passed in script: %s" % basename(__file__))
Example 0.2 - The camera class¶
Introduction into the Camera
object using the example of the ECII camera standard.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Introduction script 2 - The Camera class.
This script introduces the camera class which is of fundamental importance
for image import (e.g. separating on, off, dark background images, etc.) and
also for the data analysis as it includes all relevant detector specifications,
such as number of pixels, pixel size, focal length, etc.
In this script, a newer version of the camera type "ecII" is created manually
in order to illustrate all relevant parameters. The only difference to the
classic ecII camera is, that the filter setup is different.
See also here for more information:
https://pyplis.readthedocs.io/en/latest/tutorials.html#data-import-specifying-
custom-camera-information
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version, OPTPARSE
from os.path import basename
from numpy.testing import assert_array_equal
import pyplis
# Check script version
check_version()
# ## SCRIPT OPTIONS
# Save the new camera as default in database (cam_info.txt file that can be
# found in the "data" directory of the installation.
SAVE_TO_DATABASE = True
# ## SCRIPT FUNCTION DEFINITIONS
def create_ecII_cam_new_filters():
# Start with creating an empty Camera object
cam = pyplis.setupclasses.Camera()
# Specify the camera filter setup
# Create on and off band filters. Obligatory parameters are "type" and
# "acronym", "type" specifies the filter type ("on" or
# "off"), "acronym" specifies how to identify this filter in the file
# name. If "id" is unspecified it will be equal to "type". The parameter
# "meas_type_acro" is only important if a measurement type (e.g. M -> meas,
# C -> calib ...) is explicitely specified in the file name.
# This is not the case for ECII camera but for the HD camera,
# see specifications in file cam_info.txt for more info.
on_band = pyplis.utils.Filter(id="on", type="on", acronym="F01",
meas_type_acro="F01", center_wavelength=310)
off_band = pyplis.utils.Filter(type="off", acronym="F02",
center_wavelength=330)
# put the two filter into a list and assign to the camera
filters = [on_band, off_band]
cam.default_filters = filters
cam.prepare_filter_setup()
# Similar to the filter setup, access info for dark and offset images needs
# to be specified. The ECII typically records 4 different dark images, two
# recorded at shortest exposure time -> offset signal predominant, one at
# low and one at high read gain. The other two recorded at longest possible
# exposure time -> dark current predominant, also at low and high read gain
offset_low_gain = pyplis.utils.DarkOffsetInfo(id="offset0", type="offset",
acronym="D0L", read_gain=0)
offset_high_gain = pyplis.utils.DarkOffsetInfo(id="offset1", type="offset",
acronym="D0H", read_gain=1)
dark_low_gain = pyplis.utils.DarkOffsetInfo(id="dark0", type="dark",
acronym="D1L", read_gain=0)
dark_high_gain = pyplis.utils.DarkOffsetInfo(id="dark1", type="dark",
acronym="D1H", read_gain=1)
# put the 4 dark info objects into a list and assign to the camera
dark_info = [offset_low_gain, offset_high_gain,
dark_low_gain, dark_high_gain]
cam.dark_info = dark_info
# Now specify further information about the camera
# camera ID (needs to be unique, i.e. not included in data base, call
# pyplis.inout.get_all_valid_cam_ids() to check existing IDs)
cam.cam_id = "ecII_new_test"
# image file type
cam.file_type = "fts"
# File name delimiter for information extraction
cam.delim = "_"
# position of acquisition time (and date) string in file name after
# splitting with delimiter
cam.time_info_pos = 3
# datetime string conversion of acq. time string in file name
cam.time_info_str = "%Y%m%d%H%M%S%f"
# position of image filter type acronym in filename
cam.filter_id_pos = 4
# position of meas type info
cam.meas_type_pos = 4
# Define which dark correction type to use
# 1: determine a dark image based on image exposure time using a dark img
# (with long exposure -> dark current predominant) and a dark image with
# shortest possible exposure (-> detector offset predominant). For more
# info see function model_dark_image in processing.py module
# 2: subtraction of a dark image recorded at same exposure time than the
# actual image
cam.darkcorr_opt = 1
# If the file name also includes the exposure time, this can be specified
# here:
cam.texp_pos = "" # the ECII does not...
# the unit of the exposure time (choose from "s" or "ms")
cam.texp_unit = ""
# define the main filter of the camera (this is only important for cameras
# which include, e.g. several on band filters.). The ID need to be one of
# the filter IDs specified above
cam.main_filter_id = "on"
# camera focal length can be specified here (but does not need to be, in
# case of the ECII cam, there is no "default" focal length, so this is left
# empty)
cam.focal_length = ""
# Detector geometry
cam.pix_height = 4.65e-6 # pixel height in m
cam.pix_width = 4.65e-6 # pixel width in m
cam.pixnum_x = 1344
cam.pixnum_y = 1024
cam._init_access_substring_info()
cam.io_opts = dict(USE_ALL_FILES=False,
SEPARATE_FILTERS=True,
INCLUDE_SUB_DIRS=True,
LINK_OFF_TO_ON=True)
# Set the custom image import method
cam.image_import_method = pyplis.custom_image_import.load_ecII_fits
# That's it...
return cam
# ## SCRIPT MAIN FUNCTION
if __name__ == "__main__":
cam = create_ecII_cam_new_filters()
print(cam)
try:
# you can add the cam to the database (raises error if ID conflict
# occurs, e.g. if the camera was already added to the database)
cam.save_as_default()
except KeyError:
print("Camera already exists in database")
cam_reload = pyplis.Camera("ecII_brandnew")
# ## IMPORTANT STUFF FINISHED - everything below is of minor importance
# for educational purposes
(options, args) = OPTPARSE.parse_args()
# apply some tests. This is done only if TESTMODE is active: testmode can
# be activated globally (see SETTINGS.py) or can also be activated from
# the command line when executing the script using the option --test 1
if int(options.test):
# quick and dirty test
cam_dict_nominal = {'darkcorr_opt': 1,
'_fid_subnum_max': 1,
'_fname_access_flags': {'filter_id': False,
'meas_type': False,
'start_acq': False,
'texp': False},
'_mtype_subnum_max': 1,
'_time_info_subnum': 1,
'cam_id': 'ecII_new_test',
'delim': '_',
'file_type': 'fts',
'filter_id_pos': 4,
'focal_length': '',
'image_import_method':
pyplis.custom_image_import.load_ecII_fits,
'main_filter_id': 'on',
'meas_type_pos': 4,
'pix_height': 4.65e-06,
'pix_width': 4.65e-06,
'pixnum_x': 1344,
'pixnum_y': 1024,
'ser_no': 9999,
'texp_pos': '',
'texp_unit': '',
'time_info_pos': 3,
'time_info_str': '%Y%m%d%H%M%S%f'}
from collections import OrderedDict
geom_data_nominal = OrderedDict([('lon', None),
('lat', None),
('altitude', None),
('azim', None),
('azim_err', None),
('elev', None),
('elev_err', None),
('alt_offset', 0.0)])
arr_nominal = list(geom_data_nominal.items())
arr_nominal.extend(list(cam_dict_nominal.items()))
arr_vals = list(cam.geom_data.items())
for k in cam_dict_nominal:
arr_vals.append((k, cam.__dict__[k]))
assert_array_equal(arr_nominal, arr_vals)
print("All tests passed in script: %s" % basename(__file__))
Example 0.3 - Introduction into ImgList objects¶
Manual creation of ImgList
objects and basic features.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script 3 - manual creation of image lists.
This script gives an introduction into the creation and the handling of
ImgList objects. Image lists normally contain images of a certain type (e.g.
on, off, dark, offset). The separation by these file types is normally done
automatically within a Dataset object (see e.g. following script
ex0_4_imglists_auto.py or example script ex01_analysis_setup.py) using a
certain camera type (see previous script ex0_2_camera_setup.py).
In this example these features are NOT used but rather, an on and off band
image list is created manually without previous definition of the Camera type
(i.e. file naming convention, etc.). The script starts with creating an ImgList
containing all images of type FITS found in the test data image base directory.
On and off band lists are then consecutively extracted from this list using
the list method ``separate_by_substr_filename``.
Furthermore, some basic image preparation features of ImgList objects are
introduced (e.g. linking of lists, dark correction, automatic blurring,
cropping, size reduction).
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version, IMG_DIR, SAVEFIGS, SAVE_DIR, FORMAT, DPI,\
OPTPARSE
import pyplis
from matplotlib.pyplot import subplots, close, show
from datetime import datetime
from os.path import join, isfile, basename
from os import listdir
# Check script version
check_version()
# ## RELEVANT DIRECTORIES AND PATHS
OFFSET_FILE = join(IMG_DIR, "EC2_1106307_1R02_2015091607064723_D0L_Etna.fts")
DARK_FILE = join(IMG_DIR, "EC2_1106307_1R02_2015091607064865_D1L_Etna.fts")
if __name__ == "__main__":
close("all")
# ## Get all images in the image path which are FITS files (actually all)
all_paths = [join(IMG_DIR, f) for f in listdir(IMG_DIR) if
isfile(join(IMG_DIR, f)) and f.endswith("fts")]
# Let pyplis know that this is the ECII camera standard, so that it is
# able to import meta information from the FITS header
cam = pyplis.Camera("ecII") # loads default info for ECII camera
# ## Now put them all into an image list
# Note that the files are not separated by filter type, or dark and offset,
# etc. so the list simply contains all images of type fts which were found
# in IMG_DIR
list_all_imgs = pyplis.imagelists.ImgList(all_paths, list_id="all",
camera=cam)
# Split the list by on band file type (which is identified by acronym
# "F01" at 4th position in file name after splitting using delimiter "_")
# creates two new list objects, one containing matches (i.e. on band images
# the other containing the rest)
on_list, rest = list_all_imgs.separate_by_substr_filename(sub_str="F01",
sub_str_pos=4,
delim="_")
# now extract all off band images from the "rest" list
off_list, rest = rest.separate_by_substr_filename(sub_str="F02",
sub_str_pos=4,
delim="_")
# Link the off band list to on band list (the index in the offband list is
# changed whenever the index is changed in the on band list based on acq.
# time stamp). Note, that a Camera class was not defined here (see prev.
# script). Specify the required information in the camera (manually, not
# required if camera was defined beforehand)
on_list.camera.delim = "_"
on_list.camera.time_info_pos = 3
on_list.camera.time_info_str = "%Y%m%d%H%M%S%f"
off_list.camera.delim = "_"
off_list.camera.time_info_pos = 3
off_list.camera.time_info_str = "%Y%m%d%H%M%S%f"
on_list.link_imglist(off_list)
# ## Load dark and offset images and set them in the on-band image list
dark_img = pyplis.image.Img(DARK_FILE,
import_method=cam.image_import_method)
offset_img = pyplis.image.Img(OFFSET_FILE,
import_method=cam.image_import_method)
on_list.add_master_dark_image(dark_img)
on_list.add_master_offset_image(offset_img)
# Go to image number 100 in on-band list
on_list.goto_img(100)
# ## Change image preparation settings (these are applied on image load)
# region of interest
on_list.roi_abs = [100, 100, 1300, 900]
# activate cropping in image list
on_list.crop = 1
# scale down to pyramid level 2
on_list.pyrlevel = 2
# activate automatic dark correction in list
on_list.DARK_CORR_OPT = 1 # see previous script
on_list.darkcorr_mode = True
# blur the on band images
on_list.add_gaussian_blurring(3)
# get current on band image (no. 100 in list), all previously set
# preparation settings are applied to this image
on_img = on_list.current_img()
# get current off band image (no. 100 in list, index was automatically
# changed since it is linked to on band list. Note that this image is
# unedited
off_img = off_list.current_img()
fig, ax = subplots(1, 2, figsize=(18, 6))
on_img.show(ax=ax[0])
on_time_str = datetime.strftime(on_img.start_acq, "%Y-%m-%d %H:%M:%S")
ax[0].set_title("Current img (on-band list): %s" % on_time_str)
# show the current off band image as well (this image is unedited)
off_img.show(ax=ax[1])
off_time_str = datetime.strftime(off_img.start_acq, "%Y-%m-%d %H:%M:%S")
ax[1].set_title("Current img (off-band list): %s" % off_time_str)
on_list.edit_info()
on_list.skip_files = 3
on_list.goto_img(15)
# ## IMPORTANT STUFF FINISHED
if SAVEFIGS:
fig.savefig(join(SAVE_DIR, "ex0_3_out_1.%s" % FORMAT),
format=FORMAT, dpi=DPI)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
npt.assert_array_equal([501, 1, 500, 0],
[on_list.nof + off_list.nof,
on_list.this.edit_log["darkcorr"],
sum(on_list.this.shape),
on_list.gaussian_blurring -
on_list.this.edit_log["blurring"]])
npt.assert_allclose([402.66284],
[off_img.mean() - on_img.mean()],
rtol=1e-7, atol=0)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except Exception:
print("Use option --show 1 if you want the plots to be displayed")
Example 0.4 - Introduction into the Dataset class¶
Automatic image type separation using the Dataset
object and the ECII camera standard.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script 4 - Automatic creation of image lists.
The previous script gave an introduction into the manual creation of
``ImgList`` objects and some basic image preparation features.
In this script, a number of ImgList objects (on, off, dark low gain,
dark high gain, offset low gain, offset high gain) is created automatically
using the Camera class created in example script ex0_2_camera_setup.py (ECII
camera)
Based on the information stored in the Camera class, a MeasSetup class is
created. The latter class collects all meta information relevant for an
emission rate analysis. Apart from the camera specs, this may include source
definitions contains information about the camera specs a the
image base directory (note that in this example, start / stop acq. time stamps
are ignored, i.e. all images available in the specified directory are imported)
"""
from __future__ import (absolute_import, division)
# Imports from SETTINGS.py
from SETTINGS import check_version, IMG_DIR, OPTPARSE
import pyplis
from os.path import basename
# ## IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex0_2_camera_setup import create_ecII_cam_new_filters
# Check script version
check_version()
if __name__ == "__main__":
# create the camera which was
cam = create_ecII_cam_new_filters()
# now throw all this stuff into the BaseSetup objec
stp = pyplis.setupclasses.MeasSetup(IMG_DIR, camera=cam)
# Create a Dataset which creates separate ImgLists for all types (dark,
# offset, etc.)
ds = pyplis.dataset.Dataset(stp)
# The image lists can be accessed in different ways for instance using
# the method "all_lists", which returns a Python list containing all
# ImgList objects that were created within the Dataset
all_imglists = ds.all_lists()
# print some information about each of the lists
for lst in all_imglists:
print("list_id: %s, list_type: %s, number_of_files: %s"
% (lst.list_id, lst.list_type, lst.nof))
# single lists can be accessed using "get_list(<id>)" using a valid ID,
# e.g.:
on_list = ds.get_list("on")
off_list = ds.get_list("off")
on_list.goto_img(50) # this also changes the index in the off band list...
# ... because it is linked to the on band list (automatically set in
# Dataset)
print("\nImgLists linked to ImgList on: %s" % on_list.linked_lists.keys())
print("Current file number on / off list: %d / %d\n" % (on_list.cfn,
off_list.cfn))
# Detected dark and offset image lists are also automatically linked to the
# on and off band image list, such that dark image correction can be
# applied
on_list.darkcorr_mode = True
off_list.darkcorr_mode = True
# the current image preparation settings can be accessed via the
# edit_info method
on_list.edit_info()
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
npt.assert_array_equal([501, 2, 2368, 0, 50],
[on_list.nof + off_list.nof,
on_list.this.is_darkcorr +
off_list.this.is_darkcorr,
sum(on_list.this.shape),
on_list.gaussian_blurring -
on_list.this.edit_log["blurring"],
on_list.cfn])
npt.assert_allclose(actual=[on_list.get_dark_image().mean()],
desired=[190.56119],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
Example 0.5 - Optical flow live view¶
Live view of optical flow calculation using OpticalFlowFarneback
object (requires webcam).
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script 5 - Optical flow Farneback liveview using webcam.
Create an OpticalFlowFarneback object and activate live view (requires webcam)
"""
from __future__ import (absolute_import, division)
# Check script version
from SETTINGS import check_version
import pyplis
check_version()
flow = pyplis.plumespeed.OptflowFarneback()
flow.live_example()
Example 0.6 - Plume cross section lines¶
Creation and orientation of LineOnImage
objects for emission rate retrievals.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script no 6 - LineOnImage objects and their orientation.
This script introduces how to create LineOnImage objects within the image plane
and specify their orientation (i.e. the direction into which the normal vector
of the line points). This is mainly important for emission rate retrievals,
where, for instance velcoity displacement vectors (e.g. from an optical flow
algorithm) have to be multiplied with the normal vector of such a line (using
the dot product).
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version, SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE
from pyplis import LineOnImage
from os.path import join, basename
from matplotlib.pyplot import show, subplots, close
from matplotlib.cm import get_cmap
# Check script version
check_version()
def create_example_lines():
"""Create some exemplary lines."""
lines_r = [] # lines with orientation "right" are stored within this list
lines_l = [] # lines with orientation "left" are stored within this list
cmap = get_cmap("jet")
# horizontal line, normal orientation to the top (0 deg)
lines_r.append(LineOnImage(x0=10, y0=10, x1=90, y1=10,
normal_orientation="right",
color=cmap(0),
line_id="Line 1"))
# horizontal line, normal orientation to the bottom (180 deg)
lines_l.append(LineOnImage(10, 10, 90, 10, normal_orientation="left",
color=cmap(0),
line_id="Line 1"))
# Vertical line normal to the right (90 deg)
lines_r.append(LineOnImage(10, 10, 10, 90, normal_orientation="right",
color=cmap(50),
line_id="Line 2"))
# Vertical line normal to the left (270 deg)
lines_l.append(LineOnImage(10, 10, 10, 90, normal_orientation="left",
color=cmap(50),
line_id="Line 2"))
# Slanted line 45 degrees
lines_r.append(LineOnImage(20, 30, 50, 60, normal_orientation="right",
color=cmap(100),
line_id="Line 3"))
# Slanted line 45 degrees
lines_l.append(LineOnImage(20, 30, 50, 60, normal_orientation="left",
color=cmap(100),
line_id="Line 3"))
# Slanted line 45 degrees
lines_r.append(LineOnImage(90, 10, 10, 90, normal_orientation="right",
color=cmap(150),
line_id="Line 4"))
# Slanted line 45 degrees
lines_l.append(LineOnImage(90, 10, 10, 90, normal_orientation="left",
color=cmap(150),
line_id="Line 4"))
lines_r.append(LineOnImage(60, 20, 80, 90, normal_orientation="right",
color=cmap(200),
line_id="Line 5"))
lines_l.append(LineOnImage(60, 20, 80, 90, normal_orientation="left",
color=cmap(200),
line_id="Line 5"))
lines_r.append(LineOnImage(40, 20, 30, 90, normal_orientation="right",
color=cmap(250),
line_id="Line 6"))
lines_l.append(LineOnImage(40, 20, 30, 90, normal_orientation="left",
color=cmap(250),
line_id="Line 6"))
return lines_r, lines_l
if __name__ == "__main__":
close("all")
fig, ax = subplots(1, 2, figsize=(18, 9))
lines_r, lines_l = create_example_lines()
for k in range(len(lines_r)):
line = lines_r[k]
# print "%d: %s" %(k, line.orientation_info)
normal = line.normal_vector
lbl = "%s" % line.line_id
line.plot_line_on_grid(ax=ax[0], include_normal=1,
include_roi_rot=True, label=lbl)
for k in range(len(lines_l)):
line = lines_l[k]
normal = line.normal_vector
lbl = "%s" % line.line_id
line.plot_line_on_grid(ax=ax[1], include_normal=1,
include_roi_rot=True, label=lbl)
ax[0].set_title("Orientation right")
ax[0].legend(loc="best", fontsize=8, framealpha=0.5)
ax[0].set_xlim([0, 100])
ax[0].set_ylim([100, 0])
ax[1].set_title("Orientation left")
ax[1].legend(loc="best", fontsize=8, framealpha=0.5)
ax[1].set_xlim([0, 100])
ax[1].set_ylim([100, 0])
# ## IMPORTANT STUFF FINISHED
if SAVEFIGS:
fig.savefig(join(SAVE_DIR, "ex0_6_out_1.%s" % FORMAT),
format=FORMAT, dpi=DPI)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
npt.assert_array_equal([sum([x.length() for x in lines_l]),
sum([x.length() for x in lines_r]),
sum([sum(x.coords) for x in lines_r]),
int(sum([sum(x.normal_vector + y.normal_vector)
for x, y in zip(lines_l, lines_r)])),
int(sum([sum(x.center_pix) for x in lines_l])),
],
[459, 459, 1030, 0, 515])
npt.assert_allclose(
actual=[sum([sum(x) for x in lines_l[2].rect_roi_rot])],
desired=[368.79393923],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except Exception:
print("Use option --show 1 if you want the plots to be displayed")
Example 0.7 - Manual cell calibration¶
Manual cell calibration based on a set of background and cell images (on / off).
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis introduction script no. 5: manual cell calibration.
Perform manual cell calibration based on 3 cell images and one background
image. The calibration data consists of 3 cells which were put in front of the
lense successively and a background images both before and after the
cell images.
This script creates an empty CellCalibEngine object in which 6 cell images
are assigned (3 on and 3 off band) with their corresponding SO2 column
densities. Further, 2 background images are assigned for each filter, one
before and one after the cells were put in front of the camera. These are used
to determine tau images from the cell images. Variations in the background
intensities are corrected for (for details see manuscript).
Note, that this calibration does not include a dark correction of the images
before the calibration, therefore, the results are slightly different compared
to the results from ex05_cell_calib_auto.py.
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
import pyplis
from os.path import join, basename
from matplotlib.pyplot import close, show
from time import time
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE, IMG_DIR
# Check script version
check_version()
# SPECIFY IMAGE PATHS FOR EACH CELL AND BACKGROUND IMAGES
BG_BEFORE_ON = "EC2_1106307_1R02_2015091607000845_F01_Etna.fts"
BG_BEFORE_OFF = "EC2_1106307_1R02_2015091607001021_F02_Etna.fts"
A53_ON = "EC2_1106307_1R02_2015091607003032_F01_Etna.fts"
A53_OFF = "EC2_1106307_1R02_2015091607003216_F02_Etna.fts"
A37_ON = "EC2_1106307_1R02_2015091607005847_F01_Etna.fts"
A37_OFF = "EC2_1106307_1R02_2015091607010023_F02_Etna.fts"
A57_ON = "EC2_1106307_1R02_2015091607013835_F01_Etna.fts"
A57_OFF = "EC2_1106307_1R02_2015091607014019_F02_Etna.fts"
BG_AFTER_ON = "EC2_1106307_1R02_2015091607020256_F01_Etna.fts"
BG_AFTER_OFF = "EC2_1106307_1R02_2015091607020440_F02_Etna.fts"
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
close("all")
start = time()
# define camera for ECII custom image import
cam = pyplis.Camera("ecII")
cellcalib = pyplis.cellcalib.CellCalibEngine(cam)
# now add all cell images manually, specifying paths, the SO2 column
# densities of each cell, and the corresponding cell ID as well as image
# type (on, off)
cellcalib.set_cell_images(img_paths=join(IMG_DIR, A53_ON),
cell_gas_cd=4.15e17,
cell_id="a53", filter_id="on")
cellcalib.set_cell_images(img_paths=join(IMG_DIR, A53_OFF),
cell_gas_cd=4.15e17,
cell_id="a53", filter_id="off")
cellcalib.set_cell_images(img_paths=join(IMG_DIR, A37_ON),
cell_gas_cd=8.59e17,
cell_id="a37", filter_id="on")
cellcalib.set_cell_images(img_paths=join(IMG_DIR, A37_OFF),
cell_gas_cd=8.59e17,
cell_id="a37", filter_id="off")
cellcalib.set_cell_images(img_paths=join(IMG_DIR, A57_ON),
cell_gas_cd=1.92e18,
cell_id="a57", filter_id="on")
cellcalib.set_cell_images(img_paths=join(IMG_DIR, A57_OFF),
cell_gas_cd=1.92e18,
cell_id="a57", filter_id="off")
# put the onband background images into a Python list ....
bg_on_paths = [join(IMG_DIR, BG_BEFORE_ON), join(IMG_DIR, BG_AFTER_ON)]
# ... and add them to the calibration object ...
cellcalib.set_bg_images(img_paths=bg_on_paths, filter_id="on")
# ... same for off band background images
bg_off_paths = [join(IMG_DIR, BG_BEFORE_OFF), join(IMG_DIR, BG_AFTER_OFF)]
cellcalib.set_bg_images(img_paths=bg_off_paths, filter_id="off")
# Prepare calibration data (i.e. CellCalibData objets) for on, off and
# AA images. This function determines tau images for each cell using
# background images scaled to the present background intensities at the
# acq. time stamp of each cell using temporal interpolation of the provided
# background images. This results in 3 tau images for each filter (on, off)
# and from that, 3 AA images are determined. Each of these collection of
# 3 tau images (on, off, AA) is then stored within a CellCalibData object
# which can be accessed using e.g. cellcalib.calib_data["aa"]
cellcalib.prepare_calib_data(on_id="on", off_id="off", darkcorr=False)
stop = time()
ax = cellcalib.plot_all_calib_curves()
ax.set_title("Manual cell calibration\n(NO dark correction performed)")
aa_calib = cellcalib.calib_data["aa"]
print("Time elapsed for preparing calibration data: %.4f s"
% (stop - start))
# ## IMPORTANT STUFF FINISHED
if SAVEFIGS:
ax.figure.savefig(join(SAVE_DIR, "ex0_7_out_1.%s" % FORMAT),
format=FORMAT, dpi=DPI)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
slope, offs = aa_calib.calib_coeffs
actual = [cellcalib.tau_stacks["aa"].sum(),
cellcalib.tau_stacks["aa"].mean(),
aa_calib.cd_vec.sum(),
slope,
offs]
npt.assert_allclose(actual=actual,
desired=[1007480.56,
0.24401481,
3.194000052267778e+18,
4.779782684462987e+18,
-2.7244424951657216e+16],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Examples for emission rate analysis¶
The following scripts are directly related for emission rate analysis and build on top of each other ending with a full emission rate analysis in Example 12 - Emission rate analysis (Etna example data).
Example 1 - Creation of analysis setup and Dataset¶
This script introduces the MeasSetup
object and how it can be used to specify all relevant information to create a Dataset
for emission rate analysis.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 1 - Analysis setup for example data set.
In this script an example data set, recorded on the 16/9/15 7:06-7:22 at
Mt. Etna is setup. Most of the following example scripts will use the
information specified here.
A typical analysis setup for plume image data contains information about the
camera (e.g. optics, file naming convention, see also ex0_2_camera_setup.py),
gas source, the measurement geometry and the wind conditions. These information
generally needs to be specified by the user before the analysis. The
information is stored within a MeasSetup object which can be used as a basis
for further analysis.
If not all neccessary information is provided, a MeasSetup object will be
created nonetheless but analysis options might be limited.
Such a MeasSetup object can be used as input for Dataset objects which
creates the analysis environment. The main purpose of Dataset classes
is to automatically separate images by their type (e.g. on / off-band
images, dark, offset) and create ImgList classes from that. ImgList
classes typically contain images of one type. The Dataset also links
relevant ImgList to each other, e.g. if the camera has an off-band
filter, and off-band images can be found in the specified directory,
then it is linked to the list containing on-band images. This means,
that the current image index in the off-band list is automatically
updated whenever the index is changed in the on-band list. If
acquisition time is available in the images, then the index is updated
based on the closest acq. time of the off-band images, based on the
current on-band acq. time. If not, it is updated by index (e.g. on-band
contains 100 images and off-band list 50. The current image number in
the on-band list is 50, then the off-band index is set to 25). Also,
dark and offset images lists are linked both to on and off-band image
lists, such that the image can be correccted for dark and offset
automatically on image load (the latter needs to be activated in the
lists using ``darkcorr_mode=True``).
This script shows how to setup a MeasSetup object and create a Dataset object
from it. As example, the first image of the on-band image time series is
displayed.
The Dataset object created here is used in script ex04_prep_aa_imglist.py
which shows how to create an image list displaying AA images.
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version, IMG_DIR, OPTPARSE
import pyplis as pyplis
from datetime import datetime
from matplotlib.pyplot import show, close
# Check script version
check_version()
# SCRIPT FUNCTION DEFINITIONS
def create_dataset():
"""Initialize measurement setup and creates dataset from that."""
start = datetime(2015, 9, 16, 7, 6, 00)
stop = datetime(2015, 9, 16, 7, 22, 00)
# Define camera (here the default ecII type is used)
cam_id = "ecII"
# the camera filter setup
filters = [pyplis.utils.Filter(type="on", acronym="F01"),
pyplis.utils.Filter(type="off", acronym="F02")]
# camera location and viewing direction (altitude will be retrieved
# automatically)
geom_cam = {"lon": 15.1129,
"lat": 37.73122,
"elev": 20.0,
"elev_err": 5.0,
"azim": 270.0,
"azim_err": 10.0,
"alt_offset": 15.0,
"focal_length": 25e-3} # altitude offset (above topography)
# Create camera setup
# the camera setup includes information about the filename convention in
# order to identify different image types (e.g. on band, off band, dark..)
# it furthermore includes information about the detector specifics (e.g.
# dimension, pixel size, focal length). Measurement specific parameters
# (e.g. lon, lat, elev, azim) where defined in the dictinary above and
# can be passed as additional keyword dictionary using **geom_cam
# Alternatively, they could also be passed directly, e.g.:
# cam = pyplis.setup.Camera(cam_id, filter_list=filters, lon=15.1129,
# lat=37.73122)
cam = pyplis.setupclasses.Camera(cam_id, filter_list=filters,
**geom_cam)
# Load default information for Etna. This information is stored in
# the source_info.txt file of the Pyplis information. You may also access
# information about any volcano via the available online access to the NOAA
# database using the method pyplis.inout.get_source_info_online(source_id).
source = pyplis.setupclasses.Source("etna")
# Provide wind direction
wind_info = {"dir": 0.0,
"dir_err": 1.0}
# "dir_err" : 15.0}
# Create BaseSetup object (which creates the MeasGeometry object)
stp = pyplis.setupclasses.MeasSetup(IMG_DIR, start, stop, camera=cam,
source=source,
wind_info=wind_info)
print(stp.LINK_OFF_TO_ON)
# Create analysis object (from BaseSetup)
# The dataset takes care of finding all vali
return pyplis.Dataset(stp)
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
close("all")
ds = create_dataset()
# get on-band image list
on_list = ds.get_list("on")
on_list.goto_next()
off_list = ds.get_list("off")
# activate dark correction in both lists. Dark and offset image lists are
# automatically assigned to plume on and off-band image lists on initiation
# of the dataset object
on_list.darkcorr_mode = True
off_list.darkcorr_mode = True
print("On-band list contains %d images, current image index: %d"
% (on_list.nof, on_list.cfn))
img = on_list.current_img()
# plume distance image retrieved from MeasGeometry class...
plume_dists = on_list.plume_dists
# ...these may be overwritten or set manually if desired
on_list.plume_dists = 10000
# The same applies for the integration step lengths for emission rate
# retrievals
step_lengths = on_list.integration_step_length
on_list.integration_step_length = 1.8 # m
img_shift = img.duplicate()
# images can be shifted using the scipy.ndimage.interpolation.shift method
# this may be required for image registration in dual camera systems.
# Whether this is supposed to be done automatically can be specified using
# the REG_SHIFT_OFF option in a MeasSetup class. It may also be specified
# directly for your cam in the custom camera definition file cam_info.txt
# using io_opts:REG_SHIFT_OFF=1 (see e.g. defintion of camera with ID
# "usgs"). Also, a default registration offset can be defined here using
#
img_shift.shift(dx_abs=-30, dy_abs=55)
img_shift.show(tit="Shifted")
# Set pixel intensities below 2000 to 0 (method of Img class)
img.set_val_below_thresh(val=0, threshold=2000)
# show modified image
img.show()
print(str(img)) # image object has an informative string representation
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
actual = [plume_dists.mean(), plume_dists.std(),
on_list.get_dark_image().mean()]
npt.assert_allclose(actual=actual,
desired=[10909.873427010458, 221.48844132471388,
190.56119],
rtol=1e-7)
npt.assert_array_equal([418, 2, 2368, 1, 1, 0,
20150916070600,
20150916072200],
[on_list.nof + off_list.nof,
on_list.this.is_darkcorr +
off_list.this.is_darkcorr,
sum(on_list.this.shape),
on_list.cfn,
off_list.cfn,
sum(img.img[img.img < 2000]),
int(ds.setup.start.strftime("%Y%m%d%H%M%S")),
int(ds.setup.stop.strftime("%Y%m%d%H%M%S"))])
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 2 - Measurement Geometry¶
This script introduces the MeasGeometry
object including some basic features.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 2 - Features of the MeasGeometry class.
In this script, some important features of the MeasGeometry class are
introduced. The class itself is automatically created in the MeasSetup
object which was created in example script ex01_analysis_setup.py and was
passed as input for a Dataset object. The relevant MeasGeometry class is stored
within the Dataset object and can be accessed via the ``meas_geometry``
attribute.
As a first feature, the viewing direction of the camera is retrieved from the
image using the position of the south east (SE) crater of Mt.Etna. The result
is then visualized in a 2D map to give an overview of the geometry. The map
further includes the initial viewing direction (see example script
ex01_analysis_setup.py) which was logged in the field using a compass and a
mechanical inclinometer.
Further, the distance to the plume is retrieved on a pixel basis (represented
as image).
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
from geonum import GeoPoint
from matplotlib.pyplot import subplots, show, close
from os.path import join
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex01_analysis_setup import create_dataset
# Check script version
check_version()
# SCRIPT FUNCTION DEFINITIONS
def find_viewing_direction(meas_geometry, draw_result=True):
"""Correct viewing direction using location of Etna SE crater.
Defines location of Etna SE crater within images (is plotted into current
plume onband image of dataset) and uses its geo location to retrieve the
camera viewing direction
:param meas_geometry: :class:`MeasGeometry` object
"""
# Position of SE crater in the image (x, y)
se_crater_img_pos = [806, 736]
# Geographic position of SE crater (extracted from Google Earth)
# The GeoPoint object (geonum library) automatically retrieves the altitude
# using SRTM data
se_crater = GeoPoint(37.747757, 15.002643, name="SE crater")
print("Retrieved altitude SE crater (SRTM): %s" % se_crater.altitude)
# The following method finds the camera viewing direction based on the
# position of the south east crater.
new_elev, new_azim, _, basemap =\
meas_geometry.find_viewing_direction(pix_x=se_crater_img_pos[0],
pix_y=se_crater_img_pos[1],
# for uncertainty estimate
pix_pos_err=100,
geo_point=se_crater,
draw_result=draw_result,
update=True) # overwrite settings
print("Updated camera azimuth and elevation in MeasGeometry, new values: "
"elev = %.1f, azim = %.1f" % (new_elev, new_azim))
return meas_geometry, basemap
def plot_plume_distance_image(meas_geometry):
"""Determine and plot image plume distance and pix-to-pix distance imgs."""
# This function returns three images, the first corresponding to pix-to-pix
# distances in horizontal direction and the second (ignored here) to
# the vertical (in this case they are the same since pixel height and
# pixel width are the same for this detector). The third image gives
# camera to plume distances on a pixel basis
(dist_img, _, plume_dist_img) =\
meas_geometry.compute_all_integration_step_lengths()
fig, ax = subplots(2, 1, figsize=(7, 8))
# Show pix-to-pix distance image
dist_img.show(cmap="gray", ax=ax[0], zlabel="Pix-to-pix distance [m]")
ax[0].set_title("Parameterised pix-to-pix dists")
# Show plume distance image (convert pixel values to from m -> km)
(plume_dist_img / 1000.0).show(cmap="gray",
ax=ax[1], zlabel="Plume distance [km]")
ax[1].set_title("Plume dists")
return fig
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
close("all")
# Create the Dataset object (see ex01)
ds = create_dataset()
# execute function defined above (see above for definition and information)
geom_corr, map_ = find_viewing_direction(ds.meas_geometry)
# execute 2. script function (see above for definition and information)
fig = plot_plume_distance_image(ds.meas_geometry)
# You can compute the plume distance for the camera CFOV pixel column just
# by calling the method plume_dist() without specifying the input azimuth
# angle...
plume_dist_cfov = geom_corr.plume_dist()[0][0]
# ... and the corresponding uncertainty
plume_dist_err_cfov = geom_corr.plume_dist_err()
# You can also retrieve an array containing the camera azimuth angles for
# each pixel column...
all_azimuths = geom_corr.all_azimuths_camfov()
# ... and use this to compute plume distances on a pixel column basis
plume_dists_all_cols = geom_corr.plume_dist(all_azimuths)
# If you want, you can update information about camera, source or
# meteorolgy using either of the following methods (below we apply a
# change in the wind-direction such that the plume propagation direction
# is changed from S to SE )
# geom_corr.update_cam_specs()
# geom_corr.update_source_specs()
geom_corr.update_wind_specs(dir=315)
geom_corr.draw_map_2d() # this figure is only displayed and not saved
# recompute plume distance of CFOV pixel
plume_dist_cfov_new = geom_corr.plume_dist()[0][0]
print("Comparison of plume distances after change of wind direction:\n"
"Previous: %.3f m\n"
"New: %.3f m" % (plume_dist_cfov, plume_dist_cfov_new))
# Using this a
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
map_.ax.figure.savefig(join(SAVE_DIR, "ex02_out_1.%s" % FORMAT),
format=FORMAT, dpi=DPI)
fig.savefig(join(SAVE_DIR, "ex02_out_2.%s" % FORMAT), format=FORMAT,
dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
npt.assert_array_equal([],
[])
# check some propoerties of the basemap (displayed in figure)
# map basemap coordinates to lon / lat values
lon, lat = map_(8335, 9392, inverse=True)
npt.assert_allclose(actual=[lon, lat, map_.delta_lon, map_.delta_lat],
desired=[15.075131135, 37.76678834,
0.149263852982,
0.118462659944],
rtol=1e-7)
# check some basic properties / values of the geometry
npt.assert_allclose(actual=[geom_corr.cam_elev,
geom_corr.cam_elev_err,
geom_corr.cam_azim,
geom_corr.cam_azim_err,
plume_dist_cfov,
plume_dist_err_cfov,
plume_dists_all_cols.mean(),
plume_dist_cfov_new],
desired=[1.547754e+01,
1.064556e+00,
2.793013e+02,
1.065411e+00,
1.073102e+04,
1.645586e+02,
1.076047e+04,
9.961593e+03],
rtol=1e-6)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 3 - Plume background analysis¶
Introduction into PlumeBackgroundModel
object the default modes for retrieval of plume background intensities and plume optical density images.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 3 - Plume background analysis.
This example script introduces features related to plume background modelling
and tau image calculations.
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
import numpy as np
from os.path import join
import pyplis
from matplotlib.pyplot import show, subplots, close
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, OPTPARSE
# Check script version
check_version()
# SCRIPT OPTIONS
# If this is True, then sky reference areas are set in auto mode (note that
# in this case, the tests at the end of the script will fail!)
USE_AUTO_SETTINGS = False
# intensity threshold to init mask for bg surface fit
POLYFIT_2D_MASK_THRESH = 2600
# Choose the background correction modes you want to use
BG_CORR_MODES = [0, # 2D poly surface fit (without sky radiance image)
1, # Scaling of sky radiance image
4, # Scaling + linear gradient correction in x & y direction
6] # Scaling + quadr. gradient correction in x & y direction
# Image file paths relevant for this script
PLUME_FILE = join(IMG_DIR, 'EC2_1106307_1R02_2015091607065477_F01_Etna.fts')
BG_FILE = join(IMG_DIR, 'EC2_1106307_1R02_2015091607022602_F01_Etna.fts')
OFFSET_FILE = join(IMG_DIR, 'EC2_1106307_1R02_2015091607064723_D0L_Etna.fts')
DARK_FILE = join(IMG_DIR, 'EC2_1106307_1R02_2015091607064865_D1L_Etna.fts')
# SCRIPT FUNCTION DEFINITIONS
def init_background_model():
"""Create background model and define relevant sky reference areas."""
# Create background modelling object
m = pyplis.plumebackground.PlumeBackgroundModel()
# Define default gas free areas in plume image
w, h = 40, 40 # width/height of rectangles
m.scale_rect = [1280, 20, 1280 + w, 20 + h]
m.xgrad_rect = [20, 20, 20 + w, 20 + h]
m.ygrad_rect = [1280, 660, 1280 + w, 660 + h]
# Define coordinates of horizontal and vertical profile lines
# row number of profile line for horizontal corretions in the sky
# gradient...
m.xgrad_line_rownum = 40
# ... and start / stop columns for the corrections
m.xgrad_line_startcol = 20
m.xgrad_line_stopcol = 1323
# col number of profile line for vertical corretions in the sky gradient...
m.ygrad_line_colnum = 1300
# ... and start / stop rows for the corrections
m.ygrad_line_startrow = 10
m.ygrad_line_stoprow = 700
# Order of polyonmial fit applied for the gradient correction
m.ygrad_line_polyorder = 2
return m
def load_and_prepare_images():
"""Load images defined above and prepare them for the background analysis.
Returns
-------
- Img, plume image
- Img, plume image vignetting corrected
- Img, sky radiance image
"""
# get custom load method for ECII
fun = pyplis.custom_image_import.load_ecII_fits
# Load the image objects and peform dark correction
plume, bg = pyplis.Img(PLUME_FILE, fun), pyplis.Img(BG_FILE, fun)
dark, offset = pyplis.Img(DARK_FILE, fun), pyplis.Img(OFFSET_FILE, fun)
# Model dark image for tExp of plume image
dark_plume = pyplis.image.model_dark_image(plume.meta["texp"],
dark, offset)
# Model dark image for tExp of background image
dark_bg = pyplis.image.model_dark_image(bg.meta["texp"],
dark, offset)
plume.subtract_dark_image(dark_plume)
bg.subtract_dark_image(dark_bg)
# Blur the images (sigma = 1)
plume.add_gaussian_blurring(1)
bg.add_gaussian_blurring(1)
# Create vignetting correction mask from background image
vign = bg.img / bg.img.max() # NOTE: potentially includes y & x gradients
plume_vigncorr = pyplis.Img(plume.img / vign)
return plume, plume_vigncorr, bg
def autosettings_vs_manual_settings(bg_model):
"""Perform automatic retrieval of sky reference areas.
If you are lazy... (i.e. you dont want to define all these reference areas)
then you could also use the auto search function, a comparison is plotted
here.
"""
auto_params = pyplis.plumebackground.find_sky_reference_areas(plume)
current_params = bg_model.settings_dict()
fig, axes = subplots(1, 2, figsize=(16, 6))
axes[0].set_title("Manually set parameters")
pyplis.plumebackground.plot_sky_reference_areas(plume, current_params,
ax=axes[0])
pyplis.plumebackground.plot_sky_reference_areas(plume, auto_params,
ax=axes[1])
axes[1].set_title("Automatically set parameters")
return auto_params, fig
def plot_pcs_profiles_4_tau_images(tau0, tau1, tau2, tau3, pcs_line):
"""Plot PCS profiles for all 4 methods."""
fig, ax = subplots(1, 1)
tau_imgs = [tau0, tau1, tau2, tau3]
for k in range(4):
img = tau_imgs[k]
profile = pcs_line.get_line_profile(img)
ax.plot(profile, "-", label=r"Mode %d: $\phi=%.3f$"
% (BG_CORR_MODES[k], np.mean(profile)))
ax.grid()
ax.set_ylabel(r"$\tau_{on}$", fontsize=20)
ax.set_xlim([0, pcs_line.length()])
ax.set_xticklabels([])
ax.set_xlabel("PCS", fontsize=16)
ax.legend(loc="best", fancybox=True, framealpha=0.5, fontsize=12)
return fig
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
close("all")
# Create a background model with relevant sky reference areas
bg_model = init_background_model()
# Define exemplary plume cross section line
pcs_line = pyplis.LineOnImage(x0=530,
y0=730,
x1=890,
y1=300,
line_id="example PCS",
color="lime")
plume, plume_vigncorr, bg = load_and_prepare_images()
auto_params, fig0 = autosettings_vs_manual_settings(bg_model)
# Script option
if USE_AUTO_SETTINGS:
bg_model.update(**auto_params)
# Model 4 exemplary tau images
# list to store figures of tau plotted tau images
_tau_figs = []
# mask for corr mode 0 (i.e. 2D polyfit)
mask = np.ones(plume_vigncorr.img.shape, dtype=np.float32)
mask[plume_vigncorr.img < POLYFIT_2D_MASK_THRESH] = 0
# First method: retrieve tau image using poly surface fit
tau0 = bg_model.get_tau_image(plume_vigncorr,
mode=BG_CORR_MODES[0],
surface_fit_mask=mask,
surface_fit_polyorder=1)
# Plot the result and append the figure to _tau_figs
_tau_figs.append(bg_model.plot_tau_result(tau0, PCS=pcs_line))
# Second method: scale background image to plume image in "scale" rect
tau1 = bg_model.get_tau_image(plume, bg, mode=BG_CORR_MODES[1])
_tau_figs.append(bg_model.plot_tau_result(tau1, PCS=pcs_line))
# Third method: Linear correction for radiance differences based on two
# rectangles (scale, ygrad)
tau2 = bg_model.get_tau_image(plume, bg, mode=BG_CORR_MODES[2])
_tau_figs.append(bg_model.plot_tau_result(tau2, PCS=pcs_line))
# 4th method: 2nd order polynomial fit along vertical profile line
# For this method, determine tau on tau off and AA image
tau3 = bg_model.get_tau_image(plume, bg, mode=BG_CORR_MODES[3])
_tau_figs.append(bg_model.plot_tau_result(tau3, PCS=pcs_line))
fig6 = plot_pcs_profiles_4_tau_images(tau0, tau1, tau2, tau3, pcs_line)
if SAVEFIGS:
fig0.savefig(join(SAVE_DIR, "ex03_out_1.%s" % FORMAT), format=FORMAT,
dpi=DPI, transparent=True)
for k in range(len(_tau_figs)):
# _tau_figs[k].suptitle("")
_tau_figs[k].savefig(join(SAVE_DIR, "ex03_out_%d.%s"
% ((k + 2), FORMAT)),
format=FORMAT, dpi=DPI)
fig6.savefig(join(SAVE_DIR, "ex03_out_6.%s" % FORMAT), format=FORMAT,
dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
m = bg_model
# test settings for clear sky reference areas
npt.assert_array_equal([2680, 3960, 160, 6, 1300, 10, 700, 40, 20,
1323, 567584],
[sum(m.scale_rect),
sum(m.ygrad_rect),
sum(m.xgrad_rect),
m.mode,
m.ygrad_line_colnum,
m.ygrad_line_startrow,
m.ygrad_line_stoprow,
m.xgrad_line_rownum,
m.xgrad_line_startcol,
m.xgrad_line_stopcol,
int(m.surface_fit_mask.sum())])
m.update(**auto_params)
# test settings for clear sky reference areas
npt.assert_array_equal([2682, 4142, 1380, 6, 1337, 1, 790, 6, 672,
1343, 567584],
[sum(m.scale_rect),
sum(m.ygrad_rect),
sum(m.xgrad_rect),
m.mode,
m.ygrad_line_colnum,
m.ygrad_line_startrow,
m.ygrad_line_stoprow,
m.xgrad_line_rownum,
m.xgrad_line_startcol,
m.xgrad_line_stopcol,
int(m.surface_fit_mask.sum())])
# test all tau-modelling results
actual = [tau0.mean(), tau1.mean(), tau2.mean(), tau3.mean()]
npt.assert_allclose(actual=actual,
desired=[0.11395558008662043,
0.25279653,
0.13842879832119934,
0.13943940574220634],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 4 - Preparation of AA image list¶
Script showing how to prepare an ImgList
containing on-band plume images, such that aa_mode
can be activated (i.e. images are loaded as AA images).
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 4 - Prepare AA image list (from onband list).
Script showing how to work in AA mode using ImgList object
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
import pyplis
from matplotlib.pyplot import close
from os.path import join
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, OPTPARSE, PCS1
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex01_analysis_setup import create_dataset
from ex02_meas_geometry import find_viewing_direction
# Check script version
check_version()
# SCRIPT FUNCTION DEFINITIONS
def prepare_aa_image_list(bg_corr_mode=6):
"""Get and prepare onband list for aa image mode.
The relevant gas free areas for background image modelling are set
automatically (see also ex. 3 for details)
:return: - on list in AA mode
"""
dataset = create_dataset()
geom, _ = find_viewing_direction(dataset.meas_geometry, False)
# Set plume background images for on and off
# this is the same image which is also used for example script
# ex03_plume_background.py demonstrating the plume background routines
path_bg_on = join(IMG_DIR,
'EC2_1106307_1R02_2015091607022602_F01_Etna.fts')
path_bg_off = join(IMG_DIR,
'EC2_1106307_1R02_2015091607022820_F02_Etna.fts')
# Get on and off lists and activate dark correction
on_lst = dataset.get_list("on")
off_list = dataset.get_list("off")
# Deactivate automatic reload in list while changing some list
# attributes
on_lst.auto_reload = False
off_list.auto_reload = False
on_lst.darkcorr_mode = True
off_list.darkcorr_mode = True
# Prepare on and offband background images
bg_on = pyplis.Img(path_bg_on)
bg_on.subtract_dark_image(on_lst.get_dark_image())
bg_off = pyplis.Img(path_bg_off)
bg_off.subtract_dark_image(off_list.get_dark_image())
# set the background images within the lists
on_lst.set_bg_img(bg_on)
off_list.set_bg_img(bg_off)
# automatically set gas free areas
on_lst.bg_model.set_missing_ref_areas(on_lst.current_img())
# Now update some of the information from the automatically set sky ref
# areas
on_lst.bg_model.xgrad_line_startcol = 20
on_lst.bg_model.xgrad_line_rownum = 25
off_list.bg_model.xgrad_line_startcol = 20
off_list.bg_model.xgrad_line_rownum = 25
# set background modelling mode
on_lst.bg_model.mode = bg_corr_mode
off_list.bg_model.mode = bg_corr_mode
on_lst.aa_mode = True # activate AA mode
off_list.auto_reload = True
on_lst.auto_reload = True
print("INITIATED AA LIST")
on_lst.meas_geometry = geom
return on_lst
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
from matplotlib.pyplot import show
from time import time
close("all")
aa_list = prepare_aa_image_list()
t0 = time()
# Deactivate auto reload while changing some settings (else, after each
# of the following operations the images are reloaded and edited, which)
aa_list.auto_reload = False
aa_list.goto_img(50)
aa_list.add_gaussian_blurring(1)
# aa_list.pyrlevel = 2
aa_list.roi_abs = [300, 300, 1120, 1000]
aa_list.crop = True
# now reactive image reload in list (loads image no. 50 with all changes
# that where set in the previous lines)
aa_list.auto_reload = True
# store some results for tests below
shape_log, mean_log = sum(aa_list.this.shape), aa_list.this.mean()
ax = aa_list.show_current(zlabel=r"$\tau_{AA}$")
print("Elapsed time: %s s" % (time() - t0))
aa_list.crop = False
ax1 = aa_list.bg_model.plot_sky_reference_areas(aa_list.current_img())
fig = aa_list.bg_model.plot_tau_result(aa_list.current_img())
# import plume cross section and computed integrated optical density
# for current image
img = aa_list.this
ax2 = img.show(vmin=-0.1, vmax=0.3)
pcs = PCS1.convert(to_pyrlevel=0)
pcs.plot_line_on_grid(ax=ax2)
pix_steps = aa_list.meas_geometry.compute_all_integration_step_lengths()[0]
integrated_aa = pcs.integrate_profile(img, pix_steps)
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
ax.figure.savefig(join(SAVE_DIR, "ex04_out_1.%s" % FORMAT),
format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
m = aa_list.bg_model
npt.assert_array_equal([2682, 4144, 1380, 6, 1337, 1, 791, 25, 20,
1343],
[sum(m.scale_rect),
sum(m.ygrad_rect),
sum(m.xgrad_rect),
m.mode,
m.ygrad_line_colnum,
m.ygrad_line_startrow,
m.ygrad_line_stoprow,
m.xgrad_line_rownum,
m.xgrad_line_startcol,
m.xgrad_line_stopcol])
actual = [aa_list.meas_geometry.cam_elev,
aa_list.meas_geometry.cam_azim,
aa_list.meas_geometry.plume_dist()[0, 0],
aa_list.this.mean(),
shape_log, mean_log]
npt.assert_allclose(actual=actual,
desired=[15.477542212645357,
279.30130009369515,
10731.024327931776,
0.009083584068527644,
1520,
0.014380159209694215],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 5 - Automatic cell calibration¶
This scripts shows how to perform automatic cell calibration based on a time series of on and off band images containing both, suitable background images and images from different SO2 cells (in both wavelength channels, cf. Example 0.7 - Manual cell calibration).
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 5 - Automatic cell calibration.
Script showing how to use the automatic cell calibration engine which only
requires to specify start / stop time stamps of a calibration window. Based on
that sub time windows for each cell as well as suitable background images are
detected and separated into individual image lists (for all filters, i.e. here
on / off).
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
import pyplis
from datetime import datetime
from time import time
from os.path import join
from matplotlib.pyplot import show, close
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, OPTPARSE
# Check script version
check_version()
# File path for storing cell AA calibration data including sensitivity
# correction mask
AA_CALIB_FILE = join(SAVE_DIR, "ex05_cellcalib_aa.fts")
# SCRIPT FUNCTION DEFINITIONS
def perform_auto_cell_calib():
# Calibration time stamps
start = datetime(2015, 9, 16, 6, 59, 00)
stop = datetime(2015, 9, 16, 7, 3, 00)
# Gas CDs in cells and cell ids
# See supplementary package data about DOAS fit retrieval
calib_cells = {'a37': [8.59e17, 2.00e17],
'a53': [4.15e17, 1.00e17],
'a57': [19.24e17, 3.00e17]}
# the camera used
cam_id = "ecII"
# The camera filter setup is different from the ECII default setup and is
# therefore defined explicitely
filters = [pyplis.utils.Filter(type="on", acronym="F01"),
pyplis.utils.Filter(type="off", acronym="F02")]
# create camera setup, this includes the filename convention for
# image separation
cam = pyplis.setupclasses.Camera(cam_id=cam_id, filter_list=filters)
# Create CellCalibSetup class for initiation of CellCalib object
setup = pyplis.setupclasses.MeasSetup(IMG_DIR, start, stop,
camera=cam,
cell_info_dict=calib_cells)
# Create CellCalibEngine object, read on...
# This is a DataSet object and performs file separation and creation of
# on / off, dark / offset lists for all images in the specified time window
c = pyplis.cellcalib.CellCalibEngine(setup)
# the following high level method calls several functions in the
# CellCalibEngine class, most importantly the method find_cells for on and
# off band image time series, which detects sub time windows for each cell
# and background images. After the individual time windows are detected for
# each cell and filter, the method _assign_calib_specs is called, which
# assigns the SO2 CD amount (specified above in dictionary calib_cells)
# to the detected sub time windows (both for on and off) based on the depth
# of the intensity dip (in the onband) for each sub time window (should
# become clear from the plot produced in this script). Then it creates
# CellImgList objects for each of the cells and for the detected background
# images (i.e. resulting in (M + 1) x 2 lists, with M being the number of
# detected intensity dips, the + 1 is the corresponding background list and
# times 2 for on / off)
c.find_and_assign_cells_all_filter_lists()
return c
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
close("all")
start = time()
c = perform_auto_cell_calib()
# prepare CellCalibData objects for on, off and aa images. These can
# be accessed via c.calib_data[key] where key is "on", "off", "aa"
c.prepare_calib_data(pos_x_abs=None, # change for a specific pix
pos_y_abs=None, # change for a specific pix
radius_abs=1, # radius of retrieval disk
on_id="on", # ImgList ID of onband filter
off_id="off") # ImgList ID of offband filter
stop = time()
# Plot search result of on
ax0 = c.plot_cell_search_result("on", include_tit=False)
ax1 = c.plot_cell_search_result("off", include_tit=False)
# Plot all calibration curves for center pixel and in a radial
# neighbourhood of 20 pixels
ax2 = c.plot_all_calib_curves()
ax2.set_xlim([0, 0.7])
ax2.set_ylim([0, 2.25e18])
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
ax0.figure.savefig(join(SAVE_DIR, "ex05_2_out_1.%s" % FORMAT),
format=FORMAT, dpi=DPI)
ax1.figure.savefig(join(SAVE_DIR, "ex05_2_out_2.%s" % FORMAT),
format=FORMAT, dpi=DPI)
ax2.figure.savefig(join(SAVE_DIR, "ex05_2_out_3.%s" % FORMAT),
format=FORMAT, dpi=DPI)
ax0.set_title("Cell search result on band", fontsize=18)
ax1.set_title("Cell search result off band", fontsize=18)
ax2.set_title("Calibration polynomials", fontsize=18)
# You can explicitely access the individual CellCalibData objects
aa_calib = c.calib_data["aa"]
aa_calib.fit_calib_data()
c.plot_calib_curve("on")
mask = c.get_sensitivity_corr_mask("aa")
aa_calib.save_as_fits(AA_CALIB_FILE)
print("Time elapsed for preparing calibration data: %.4f s"
% (stop - start))
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
calib_reload = pyplis.CellCalibData()
calib_reload.load_from_fits(AA_CALIB_FILE)
calib_reload.fit_calib_data(polyorder=2,
through_origin=True)
# test some basic features of calibraiton dataset (e.g. different
# ImgList classes for on and off and the different cells)
npt.assert_array_equal([c.cell_search_performed,
c.cell_lists["on"]["a37"].nof,
c.cell_lists["on"]["a53"].nof,
c.cell_lists["on"]["a57"].nof,
c.cell_lists["off"]["a37"].nof,
c.cell_lists["off"]["a53"].nof,
c.cell_lists["off"]["a57"].nof,
calib_reload.calib_id,
calib_reload.pos_x_abs,
calib_reload.pos_y_abs],
[True, 2, 3, 3, 2, 3, 3, "aa", 672, 512])
d = c._cell_info_auto_search
vals_approx = [d["a37"][0],
d["a53"][0],
d["a57"][0],
aa_calib.calib_fun(0, *aa_calib.calib_coeffs),
calib_reload.calib_fun(0, *calib_reload.calib_coeffs)]
npt.assert_allclose(actual=vals_approx,
desired=[8.59e+17,
4.15e+17,
1.924e+18,
-1.8313164653590516e+16,
0.0],
rtol=1e-7)
# explicitely check calibration data for on, off and aa (plotted in
# this script)
actual = [c.calib_data["on"].calib_coeffs.mean(),
c.calib_data["off"].calib_coeffs.mean(),
aa_calib.calib_coeffs.mean(),
calib_reload.calib_coeffs.mean()]
npt.assert_allclose(actual=actual,
desired=[1.8926803829028974e+18,
-3.742539080945949e+19,
2.153653759747737e+18,
2.1197681750384312e+18],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 6 - DOAS calibration¶
Introduction into DOAS calibration including FOV search using both, the Pearson and the IFR method.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 6 - DOAS calibration and FOV search.
Script showing how to work with DOAS calibration data
In this script, a stack of plume AA images from the Etna test data is imported
as well as a time series of DOAS SO2-CDs in the plume retrieved using the
software DOASIS (see directory "spectra" in test data folder for corresponding
analysis details, the folder also contains the RAW data and the jscript code
for analysing the spectra). The DOAS result import is performed using the
Python package ``pydoas``.
Based on these data, position and shape of the DOAS FOV within the camera
uimages is identified using both FOV search methods (IFR and Pearson). The
results of the FOV search are plotted as well as the corresponding calibration
curves retrieved for both FOV parametrisations.
Note
------
In case a MemoryError occurs while determining the AA image stack, then the
stack (3D numpy array) is too large for the RAM. In this case, try
increasing script option PYRLEVEL_ROUGH_SEARCH.
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
import pyplis
import pydoas
import numpy.testing as npt
import numpy as np
from datetime import timedelta
from matplotlib.pyplot import close, show, subplots
from os.path import join, exists
from os import remove
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, IMG_DIR, OPTPARSE
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list
# Check script version
check_version()
# SCRIPT OPTONS
# reload and save stack in folder SAVE_DIR, results in increased
# running time due to stack calculation (is automatically switched on if
# the stack is not found at this location)
RELOAD_STACK = False
PYRLEVEL_ROUGH_SEARCH = 2
# Default search settings are at pyramid level 2, the FOV results are upscaled
# to original resolution, if the following option is set 1, then, based on
# the result from pyrlevel=2, another stack is determined at pyrlevel = 0
# (i.e. in full resolution) within a ROI around the center position from the
# search performed at pyrlevel=PYRLEVEL_ROUGH_SEARCH (defined 2 lines below)
DO_FINE_SEARCH = False
# RELEVANT DIRECTORIES AND PATHS
# Directory containing DOAS result files
DOAS_DATA_DIR = join(IMG_DIR, "..", "spectra", "plume_prep", "min10Scans",
"ResultFiles")
STACK_PATH = join(SAVE_DIR, "ex06_aa_imgstack.fts")
# SCRIPT FUNCTION DEFINITIONS
def load_doas_results(lt_to_utc_shift=timedelta(-1. / 12)):
"""Specify DOAS data import from DOASIS fit result files.
In order to perform the DOAS FOV search, as many spectrum datapoints
as possible are needed. Therefore, only 10 spectra were added (to reduce
noise) per plume spectrum. The DOAS fit was performed in a wavelength
range between 314 - 326 nm (fit ID: "f01").
"""
# This dictionary specifies which information is supposed to be imported
# from the DOAS fit result files stored in DOAS_DATA_DIR. In the example
# shown here, only the SO2 fit results are imported from fit scenario
# with ID "f01" (key of dict). The corresponding value of each key is
# a list of format ["header_id", ["fit_id1", "fit_id2", ..., "fit_idN"]]
# specifying the identification string of the species in the result file
# headers and the second entry is a list specifying all fit scenario IDs
# from which this species is supposed to be imported (here only f01)
fit_import_info = {"so2": ["SO2_Hermans_298_air_conv_satCorr1e18",
["f01"]
]}
# Create a result import setup for the DOAS data based on the import
# dictionary and the image base directory of the result files ...
doas_import_setup =\
pydoas.dataimport.ResultImportSetup(DOAS_DATA_DIR,
result_import_dict=fit_import_info)
# ... and create a result dataset from that
doas_dataset = pydoas.analysis.DatasetDoasResults(doas_import_setup)
# get the SO2 fit results from the dataset. Individual results of certain
# species can be accessed using the species ID (key in ``fit_import_info``
# dict) and its fit ID (one of the fit IDs specified for this species, here
# f01).
# Note, that the DOAS data was stored using local time, thus they need to
# be shifted (2h back) to match the camera data time stamps (which are in
# UTC), otherwise the temporal merging of the two datasets (for the DOAS
# calibration) does not work
results_utc = doas_dataset.get_results("so2", "f01").shift(lt_to_utc_shift)
return results_utc
def make_aa_stack_from_list(aa_list, roi_abs=None, pyrlevel=None,
save=True, stack_path=STACK_PATH,
save_dir=SAVE_DIR):
"""Get and prepare onband list for aa image mode."""
# Deactivate auto reload to change some settings (if auto_reload is active
# list images are reloaded whenever a setting is changed in the list. This
# can slow down things, thus, if you intend to change a couple of settings
# you might deactivate auto_reload, adapt the settings and then re-activate
# auto_reload
aa_list.auto_reload = False
if roi_abs is not None:
aa_list.roi_abs = roi_abs
aa_list.crop = True
aa_list.pyrlevel = pyrlevel
aa_list.auto_reload = True
# Stack all images in image list at pyrlevel 2 and cropped using specified
# roi (uncropped if roi_abs=None).
stack = aa_list.make_stack()
if save:
try:
remove(stack_path)
except BaseException:
pass
stack.save_as_fits(save_dir=save_dir,
save_name="ex06_aa_imgstack.fts")
return stack
def get_stack(reload_stack=RELOAD_STACK, stack_path=STACK_PATH,
pyrlevel=PYRLEVEL_ROUGH_SEARCH):
"""Load stack data based on current settings."""
if not exists(stack_path):
reload_stack = 1
if not reload_stack:
stack = pyplis.processing.ImgStack()
stack.load_stack_fits(stack_path)
if stack.pyrlevel != pyrlevel:
reload_stack = True
aa_list = None
if reload_stack:
# import AA image list
aa_list = prepare_aa_image_list()
# Try creating stack
stack = make_aa_stack_from_list(aa_list, pyrlevel=pyrlevel)
return stack, aa_list
# Test functions used at the end of the script
def test_calib_pears_init(calib):
calib.fit_calib_data(polyorder=1)
cc = pyplis.helpers.get_img_maximum(calib.fov.corr_img.img)
assert cc == (124, 159), cc
pyrl = calib.fov.pyrlevel
assert pyrl == 2, pyrl
res_dict = calib.fov.result_pearson
npt.assert_allclose([res_dict['rad_rel'],
np.max(100 * res_dict['corr_curve'].values)],
[3, 95], atol=1)
fov_ext = calib.fov.pixel_extend(abs_coords=True)
(fov_x, fov_y) = calib.fov.pixel_position_center(abs_coords=True)
npt.assert_allclose([fov_ext, fov_x, fov_y],
[res_dict['rad_rel'] * 2**pyrl, 636, 496], atol=1)
npt.assert_allclose(calib.calib_coeffs,
[8.58e+18, 2.71e+17], rtol=1e-1)
def test_calib_pears_fine(calib):
cc = pyplis.helpers.get_img_maximum(calib.fov.corr_img.img)
npt.assert_allclose((186, 180), cc, atol=1)
pyrl = calib.fov.pyrlevel
assert pyrl == 0, pyrl
res_dict = calib.fov.result_pearson
npt.assert_allclose([res_dict['rad_rel'],
np.max(100 * res_dict['corr_curve'].values)],
[6, 95], atol=1)
fov_ext = calib.fov.pixel_extend(abs_coords=True)
(fov_x, fov_y) = calib.fov.pixel_position_center(abs_coords=True)
npt.assert_allclose([fov_ext, fov_x, fov_y],
[6, 630, 493], atol=1)
npt.assert_allclose(calib.calib_coeffs,
[8.38e+18, 2.92e+17], rtol=1e-1)
def test_calib_ifr(calib):
cc = pyplis.helpers.get_img_maximum(calib.fov.corr_img.img)
npt.assert_allclose((123, 157), cc, atol=1)
pyrl = calib.fov.pyrlevel
assert pyrl == 2, pyrl
npt.assert_allclose(calib.fov.result_ifr['popt'][1:5],
[158.6, 122.9, 15.4, 1.5], rtol=1e-1)
(fov_x, fov_y) = calib.fov.pixel_position_center(abs_coords=True)
npt.assert_allclose([fov_x, fov_y, calib_ifr.fov.sigma_x_abs,
calib_ifr.fov.sigma_y_abs],
[635, 492, 61.5, 41.6], atol=2)
npt.assert_allclose(calib.calib_coeffs,
[9.38e+18, 1.75e+17], rtol=1e-1)
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
# close all plots
close("all")
# import DOAS results
doas_time_series = load_doas_results()
# Import script options
(options, args) = OPTPARSE.parse_args()
if options.test:
# if test mode is active, the image stack is always recomputed from
# scratch and the option DO_FINE_SEARCH is activated, since the tests
# are based on this. This will lead to an increased computation time
# Test mode can be activated / deactivated in SETTINGS.py or via
# the script option --test 1 (on) or --test 0 (off)
RELOAD_STACK = True
PYRLEVEL_ROUGH_SEARCH = 2
DO_FINE_SEARCH = True
# reload or create the AA image stack based on current script settings
stack, aa_list = get_stack()
s = pyplis.doascalib.DoasFOVEngine(stack, doas_time_series)
calib_pears = s.perform_fov_search(method="pearson")
calib_ifr = s.perform_fov_search(method="ifr", ifrlbda=4e-3)
# plot the FOV search results
ax0 = calib_pears.fov.plot()
ax1 = calib_ifr.fov.plot()
calib_pears.fit_calib_data()
calib_ifr.fit_calib_data()
fig, ax2 = subplots(1, 1)
calib_pears.plot(add_label_str="Pearson", color="b", ax=ax2)
calib_ifr.plot(add_label_str="IFR", color="g", ax=ax2)
ax2.set_title("Calibration curves Pearson vs. IFR method")
ax2.grid()
ax2.set_ylim([0, 1.8e18])
ax2.set_xlim([0, 0.20])
ax2.legend(loc=4, fancybox=True, framealpha=0.7, fontsize=11)
axes = [ax0, ax1, ax2]
if DO_FINE_SEARCH:
"""Perform FOV search within ROI around result from pearson fov
search at full resolution (pyrlevel=0)
"""
if aa_list is None:
aa_list = prepare_aa_image_list()
# remember some properties of the current image stack that is stored in
# the DoasFOVEngine object (i.e. the merged one in low pixel res.)
num_merge, h, w = s.img_stack.shape
s_fine = s.run_fov_fine_search(aa_list, doas_time_series,
method="pearson")
calib_pears_fine = s_fine.calib_data
calib_pears_fine.plot()
calib_pears_fine.fov.plot()
calib_pears.save_as_fits(save_dir=SAVE_DIR,
save_name="ex06_doascalib_aa.fts")
calib_ifr.save_as_fits(save_dir=SAVE_DIR,
save_name="ex06_doascalib_aa_ifr_method.fts")
# you can also change the order of the calibration polynomial and
# force it to go through the origin
calib_pears.fit_calib_data(polyorder=2, through_origin=True)
calib_pears.plot_calib_fun(add_label_str="Pearson (WRONG,\n"
"2nd order, through origin)",
color="r", ax=ax2)
ax2.legend(loc=4, fancybox=True, framealpha=0.7, fontsize=11)
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
for k in range(len(axes)):
ax = axes[k]
ax.set_title("")
ax.figure.savefig(join(SAVE_DIR, "ex06_out_%d.%s"
% ((k + 1), FORMAT)),
format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
from os.path import basename
num, h, w = stack.shape
num2 = s.img_stack.shape[0] # stack after fine FOV search
prep = stack.img_prep
# check some basic properties of the data used for the different FOV
# searches
npt.assert_array_equal(
[len(doas_time_series), num, num_merge, h, w, stack.pyrlevel,
prep["darkcorr"] * prep["is_tau"] * prep["is_aa"], num2],
[120, 209, 88, 256, 336, 2, 1, 209])
test_calib_pears_init(calib_pears)
test_calib_ifr(calib_ifr)
test_calib_pears_fine(calib_pears_fine)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 7 - AA sensitivity correction masks¶
Combine the results from Example 5 - Automatic cell calibration and Example 6 - DOAS calibration in order to retrieve AA sensitivity correction masks normalised to the position of the DOAS FOV.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 7 - AA sensitivity correction masks.
In this script, cell and DOAS calibration (see previous 2 scripts) of the Etna
test dataset are opposed. Furthermore, it is illustrated, how to create
correction masks for pixel variations in the SO2 sensitivity due to shifts in
the filter transmission windows.
The cell calibration is re-performed (using method ``perform_auto_cell_calib``)
from example script 5. The results from the DOAS calibration
(see prev. example) were stored as a FITS file (including FOV information)
and the results are imported here.
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
import pyplis
from os.path import join, exists
import numpy as np
from matplotlib.pyplot import close, subplots, show
from matplotlib.patches import Circle
import six
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex05_cell_calib_auto import perform_auto_cell_calib
from ex04_prep_aa_imglist import prepare_aa_image_list
# Check script version
check_version()
CELL_AA_CALIB_FILE = join(SAVE_DIR, "ex05_cellcalib_aa.fts")
# RELEVANT DIRECTORIES AND PATHS
# fits file containing DOAS calibration information (from ex6)
DOAS_CALIB_FILE = join(SAVE_DIR, "ex06_doascalib_aa.fts")
# SCRIPT FUNCTION DEFINITIONS
def draw_doas_fov(fov_x, fov_y, fov_extend, ax):
# add FOV position to plot of examplary AA image
c = Circle((fov_x, fov_y), fov_extend, ec="k", fc="lime", alpha=.5)
ax.add_artist(c)
ax.text(fov_x, (fov_y - fov_extend * 1.3), "DOAS FOV")
ax.set_xlim([0, 1343]), ax.set_ylim([1023, 0])
return ax
def plot_pcs_comparison(aa_init, aa_imgs_corr, pcs1, pcs2):
fig, axes = subplots(1, 2, figsize=(18, 6))
p10 = pcs1.get_line_profile(aa_init.img)
p20 = pcs2.get_line_profile(aa_init.img)
num = len(p10)
axes[0].set_title("Line %s" % pcs1.line_id)
axes[1].set_title("Line %s" % pcs2.line_id)
axes[0].plot(p10, "-", label=r"Init $\phi=%.3f$" % (sum(p10) / num))
axes[1].plot(p20, "-", label=r"Init $\phi=%.3f$" % (sum(p20) / num))
for cd, aa_corr in six.iteritems(aa_imgs_corr):
p1 = pcs1.get_line_profile(aa_corr.img)
p2 = pcs2.get_line_profile(aa_corr.img)
axes[0].plot(p1, "-", label=r"Cell CD: %.2e $\phi=%.3f$"
% (cd, sum(p1) / num))
axes[1].plot(p2, "-", label=r"Cell CD: %.2e $\phi=%.3f$"
% (cd, sum(p2) / num))
axes[0].legend(loc='best', fancybox=True, framealpha=0.5, fontsize=10)
axes[1].legend(loc='best', fancybox=True, framealpha=0.5, fontsize=10)
axes[0].grid()
axes[1].grid()
return fig, axes
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
close("all")
if not exists(DOAS_CALIB_FILE):
raise IOError("Calibration file could not be found at specified "
"location:\n %s\nYou might need to run example 6 first")
# Load AA list
aa_list = prepare_aa_image_list()
aa_list.add_gaussian_blurring(2)
# Load DOAS calbration data and FOV information (see example 6)
doascalib = pyplis.doascalib.DoasCalibData()
doascalib.load_from_fits(file_path=DOAS_CALIB_FILE)
doascalib.fit_calib_data()
# Get DOAS FOV parameters in absolute coordinates
fov_x, fov_y = doascalib.fov.pixel_position_center(abs_coords=True)
fov_extend = doascalib.fov.pixel_extend(abs_coords=True)
# Load cell calibration (see example 5)
cellcalib = perform_auto_cell_calib()
# get cell calibration
cellcalib.prepare_calib_data(
pos_x_abs=fov_x, # change if you want it for a specific pix
pos_y_abs=fov_y, # change if you want it for a specific pix
radius_abs=fov_extend, # radius of retrieval disk
on_id="on", # ImgList ID of onband filter
off_id="off") # ImgList ID of offband filter
cell_aa_calib = cellcalib.calib_data["aa"]
# Define lines on image for plume profiles
pcs1 = pyplis.LineOnImage(620, 700, 940, 280,
line_id="center")
pcs2 = pyplis.LineOnImage(40, 40, 40, 600,
line_id="edge")
# Plot DOAS calibration polynomial
ax0 = doascalib.plot(add_label_str="DOAS")
ax0 = cellcalib.calib_data["aa"].plot(ax=ax0, c="r")
ax0.set_title("")
ax0.set_xlim([0, 0.5])
# Get current AA image from image list
aa_init = aa_list.current_img()
# now determine sensitivity correction masks from the different cells
masks = {}
aa_imgs_corr = {}
for cd in cell_aa_calib.cd_vec:
mask = cellcalib.get_sensitivity_corr_mask("aa",
pos_x_abs=fov_x,
pos_y_abs=fov_y,
radius_abs=fov_extend,
cell_cd_closest=cd)
masks[cd] = mask
aa_imgs_corr[cd] = pyplis.Img(aa_init.img / mask.img)
# get mask corresponding to minimum cell CD
mask = list(masks.values())[np.argmin(list(masks.keys()))]
# assing mask to aa_list
aa_list.senscorr_mask = mask
# activate AA sensitivity correction in list
aa_list.sensitivity_corr_mode = True
# set DOAS calibration data in list ...
aa_list.calib_data = doascalib
# ... and activate calibration mode
aa_list.calib_mode = True
ax = aa_list.current_img().show(zlabel=r"$S_{SO2}$ [cm$^{-2}$]")
# plot the two lines into the exemplary AA image
pcs1.plot_line_on_grid(ax=ax, color="r")
pcs2.plot_line_on_grid(ax=ax, color="g")
ax.legend(loc='best', fancybox=True, framealpha=0.5, fontsize=10)
ax = draw_doas_fov(fov_x, fov_y, fov_extend, ax=ax)
fig, _ = plot_pcs_comparison(aa_init, aa_imgs_corr, pcs1, pcs2)
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
ax0.figure.savefig(join(SAVE_DIR, "ex07_out_1.%s" % FORMAT),
format=FORMAT, dpi=DPI)
ax.figure.savefig(join(SAVE_DIR, "ex07_out_2.%s" % FORMAT),
format=FORMAT, dpi=DPI)
fig.savefig(join(SAVE_DIR, "ex07_out_3.%s" % FORMAT), format=FORMAT,
dpi=DPI)
# Save the sensitivity correction mask from the cell with the lowest SO2 CD
so2min = np.min(list(masks.keys()))
mask = masks[so2min]
mask.save_as_fits(SAVE_DIR, "ex07_aa_corr_mask")
# assign mask in doascalib and resave
doascalib.senscorr_mask = mask
doascalib.save_as_fits(DOAS_CALIB_FILE)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
npt.assert_array_equal([],
[])
npt.assert_allclose(actual=[],
desired=[],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 8 - Plume velocity retrieval (Cross correlation)¶
In this script an exemplary plume velocity retrieval is performed using the signal cross correlation algorithm. The velocity is retrieved based on two plume cross sections and a time series of plume AA images (using the AA ImgList
created in Example 4 - Preparation of AA image list).
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 8 - Plume velocity retrieval by cross correlation.
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
from matplotlib.pyplot import close, show, subplots
from os.path import join
from time import time
from pyplis.plumespeed import VeloCrossCorrEngine
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE, LINES
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list
# Check script version
check_version()
# SCRIPT OPTONS
# distance in pixels between two lines used for cross correlation analysis
OFFSET_PIXNUM = 40
RELOAD = 0 # reload AA profile images for PCS lines
# start / stop indices of considered images in image list (only relevant if PCS
# profiles are reloaded, i.e. Opt RELOAD=True)
START_IDX = 10
STOP_IDX = 200
# PCS line for which velocity is supposed to retrieved
PCS = LINES[0] # orange "young_plume" line
# Color of PCS offset line used to perform cross correlation analysis
# (relevant for illustration)
COLOR_PCS_OFFS = "c"
# RELEVANT DIRECTORIES AND PATHS
# the time series of PCS profiles for both lines are stored as
# ProfileTimeSeriesImg objects using the following names. The images
# contain the PCS profiles (y-axis) for each image in the list (y-axis)
# and have thus dimension MXN where M denotes the pixel number of the lines
# and N denotes the total number of images from which the profiles are
# extracted. The images will be stored in SAVE_DIR after this script is run
# once. After that, re-running the script and applying the cross-correlation
# analysis will be much faster, since the profiles are imported from the
# two precomupted images and do not need to be extracted by looping over
# the image list.
PCS_PROFILES_PIC_NAME = "ex08_ica_tseries_pcs.fts"
OFFSET_PROFILES_PIC_NAME = "ex08_ica_tseries_offset.fts"
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
close("all")
axes = []
# prepare the AA image list (see ex4)
aa_list = prepare_aa_image_list()
aa_list.pyrlevel = 1
t0 = time()
cc = VeloCrossCorrEngine(aa_list, PCS)
cc.create_parallel_pcs_offset(offset_pix=40,
color=COLOR_PCS_OFFS,
linestyle="--")
reloaded = False # just a flag for output below
try:
if RELOAD:
raise Exception
cc.load_pcs_profile_img(join(SAVE_DIR, PCS_PROFILES_PIC_NAME),
line_id="pcs")
cc.load_pcs_profile_img(join(SAVE_DIR, OFFSET_PROFILES_PIC_NAME),
line_id="pcs_offset")
except BaseException:
cc.get_pcs_tseries_from_list(start_idx=START_IDX,
stop_idx=STOP_IDX)
cc.save_pcs_profile_images(save_dir=SAVE_DIR,
fname1=PCS_PROFILES_PIC_NAME,
fname2=OFFSET_PROFILES_PIC_NAME)
reloaded = True
t1 = time()
# the run method of the high level VeloCrossCorrEngine class is
# basically a wrapper method for the low-level find_signal_correlation
# function which is part of the plumespeed.py module. Before calling
# the latter, the ICA time-series are extracted from the two
# ProfileTimeSeriesImg objects which were computed above from the
# ImgList class containing AA images, and which are stored as FITS
# files for fast re-computing of this script. The following run
# command passes valid input parameters to the find_signal_correlation
# method.
velo = cc.run(cut_border_idx=10,
reg_grid_tres=100,
freq_unit="L",
sigma_smooth=2,
plot=0)
t2 = time()
fig, ax = subplots(1, 2, figsize=(20, 6))
axes.append(cc.plot_pcs_lines(ax=ax[0]))
cc.plot_ica_tseries_overlay(ax=ax[1])
axes.append(cc.plot_corrcoeff_tseries())
print("Result performance analysis\n"
"Images reloaded from list: %s\n"
"Number of images: %d\n"
"Create ICA images: %.3f s\n"
"Cross-corr analysis: %.3f s"
% (reloaded, (STOP_IDX - START_IDX), (t1 - t0), (t2 - t1)))
print("Retrieved plume velocity of v = %.2f m/s" % velo)
# IMPORTANT STUFF FINISHED
if SAVEFIGS:
for k in range(len(axes)):
axes[k].figure.savefig(join(SAVE_DIR, "ex08_out_%d.%s"
% ((k + 1), FORMAT)),
format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
npt.assert_array_equal([],
[])
npt.assert_allclose(actual=[],
desired=[],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 9 - Plume velocity retrieval (Optical flow Farneback)¶
This script gives an introduction into plume velocity retrievals using the Farneback optical flow algorithm (OpticalFlowFarneback
) and a histogram based post analysis.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 9 - Optical flow Plume velocity retrieval."""
from __future__ import (absolute_import, division)
from os.path import join
import pyplis
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE, LINES
from matplotlib.pyplot import (close, show, subplots, figure, xticks, yticks,
sca, rcParams)
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list
from SETTINGS import check_version
rcParams["font.size"] = 16
PCS1, PCS2 = LINES
# Check script version
check_version()
# SCRIPT OPTIONS
PEAK_SIGMA_TOL = 2
# perform histogram analysis for all images in time series
HISTO_ANALYSIS_ALL = 1
# applies multi gauss fit to retrieve local predominant displacement
# direction, if False, then the latter is calculated from 1. and 2. moment
# of histogram (Faster but more sensitive to additional peaks in histogram)
HISTO_ANALYSIS_MULTIGAUSS = True
HISTO_ANALYSIS_START_IDX = 0
HISTO_ANALYSIS_STOP_IDX = None # 207
# Gauss pyramid level
PYRLEVEL = 1
BLUR = 0
ROI_CONTRAST = [0, 0, 1344, 730]
MIN_AA = 0.05
def analyse_and_plot(lst, lines):
fig = figure(figsize=(14, 8))
ax0 = fig.add_axes([0.01, 0.15, 0.59, 0.8])
# ax0.set_axis_off()
ax1 = fig.add_axes([0.61, 0.15, 0.16, 0.8])
ax2 = fig.add_axes([0.78, 0.15, 0.16, 0.8])
mask = lst.get_thresh_mask(MIN_AA)
fl = lst.optflow
fl.plot(ax=ax0) # , in_roi=True)
for line in lines:
m = mask * line.get_rotated_roi_mask(fl.flow.shape[:2])
line.plot_line_on_grid(ax=ax0, include_normal=1,
include_roi_rot=1)
try:
_, mu, sigma = fl.plot_orientation_histo(pix_mask=m,
apply_fit=True, ax=ax1,
color=line.color)
ax1.legend_.remove()
low, high = mu - sigma, mu + sigma
fl.plot_length_histo(pix_mask=m, ax=ax2, dir_low=low,
dir_high=high, color=line.color)
except BaseException:
pass
# pyplis.helpers.set_ax_lim_roi(roi_disp, ax0)
ax0.get_xaxis().set_ticks([])
ax0.get_yaxis().set_ticks([])
ax0.set_title("")
ymax = max([ax1.get_ylim()[1], ax2.get_ylim()[1]])
ax1.set_title("")
ax1.set_xlabel(r"$\varphi\,[^\circ]$", fontsize=20)
ax1.set_ylim([0, ymax])
ax1.get_yaxis().set_ticks([])
ax2.yaxis.tick_right()
ax2.yaxis.set_label_position("right")
ax2.set_xlabel(r"$|\mathbf{f}|$ [pix]", fontsize=20)
ax2.set_ylabel("Counts / bin", fontsize=20)
ax2.set_ylim([0, ymax])
ax2.set_title("")
ax2.legend_.remove()
sca(ax1)
xticks(rotation=40, ha="right")
sca(ax2)
yticks(rotation=90, va="center")
return fig
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
close("all")
figs = []
# Prepare aa image list (see example 4)
aa_list = prepare_aa_image_list()
# the aa image list includes the measurement geometry, get pixel
# distance image where pixel values correspond to step widths in the plume,
# obviously, the distance values depend on the downscaling factor, which
# is calculated from the analysis pyramid level (PYRLEVEL)
dist_img, _, _ = aa_list.meas_geometry.compute_all_integration_step_lengths( # noqa: E501
pyrlevel=PYRLEVEL)
# set the pyramid level in the list
aa_list.pyrlevel = PYRLEVEL
# add some blurring.. or not (if BLUR = 0)
aa_list.add_gaussian_blurring(BLUR)
# Access to the optical flow module in the image list. If optflow_mode is
# active in the list, then, whenever the list index changes (e.g. using
# list.goto_next(), or list.goto_img(100)), the optical flow field is
# calculated between the current list image and the next one
fl = aa_list.optflow
# (! note: fl is only a pointer, i.e. the "=" is not making a copy of the
# object, meaning, that whenever something changes in "fl", it also does
# in "aa_list.optflow")
# Now activate optical flow calculation in list (this slows down the
# speed of the analysis, since the optical flow calculation is
# comparatively slow
s = aa_list.optflow.settings
s.hist_dir_gnum_max = 10
s.hist_dir_binres = 10
s.hist_sigma_tol = PEAK_SIGMA_TOL
s.roi_rad = ROI_CONTRAST
aa_list.optflow_mode = True
plume_mask = pyplis.Img(aa_list.get_thresh_mask(MIN_AA))
plume_mask.show(tit="AA threshold mask")
figs.append(analyse_and_plot(aa_list, LINES))
figs.append(fl.plot_flow_histograms(PCS1, plume_mask.img))
figs.append(fl.plot_flow_histograms(PCS2, plume_mask.img))
# Show an image containing plume speed magnitudes (ignoring direction)
velo_img = pyplis.Img(fl.to_plume_speed(dist_img))
velo_img.show(vmin=0, vmax=10, cmap="Greens",
tit="Optical flow plume velocities",
zlabel="Plume velo [m/s]")
# Create two objects used to store time series information about the
# retrieved plume properties
plume_props_l1 = pyplis.plumespeed.LocalPlumeProperties(PCS1.line_id)
plume_props_l2 = pyplis.plumespeed.LocalPlumeProperties(PCS2.line_id)
if HISTO_ANALYSIS_ALL:
aa_list.goto_img(HISTO_ANALYSIS_START_IDX)
if HISTO_ANALYSIS_STOP_IDX is None:
HISTO_ANALYSIS_STOP_IDX = aa_list.nof - 1
for k in range(HISTO_ANALYSIS_START_IDX, HISTO_ANALYSIS_STOP_IDX):
plume_mask = aa_list.get_thresh_mask(MIN_AA)
plume_props_l1.get_and_append_from_farneback(
fl, line=PCS1, pix_mask=plume_mask,
dir_multi_gauss=HISTO_ANALYSIS_MULTIGAUSS)
plume_props_l2.get_and_append_from_farneback(
fl, line=PCS2, pix_mask=plume_mask,
dir_multi_gauss=HISTO_ANALYSIS_MULTIGAUSS)
aa_list.goto_next()
# ==============================================================================
# plume_props_l1 = plume_props_l1.interpolate()
# plume_props_l2 = plume_props_l2.interpolate()
# ==============================================================================
fig, ax = subplots(2, 1, figsize=(10, 9))
plume_props_l1.plot_directions(ax=ax[0],
color=PCS1.color,
label="PCS1")
plume_props_l2.plot_directions(ax=ax[0], color=PCS2.color,
label="PCS2")
plume_props_l1.plot_magnitudes(normalised=True, ax=ax[1],
date_fmt="%H:%M:%S", color=PCS1.color,
label="PCS1")
plume_props_l2.plot_magnitudes(normalised=True, ax=ax[1],
date_fmt="%H:%M:%S", color=PCS2.color,
label="PCS2")
ax[0].set_xticklabels([])
# ax[0].legend(loc='best', fancybox=True, framealpha=0.5, fontsize=14)
# ax[0].set_title("Movement direction")
# ax[1].set_title("Displacement length")
figs.append(fig)
# Save the time series as txt
plume_props_l1.save_txt(join(SAVE_DIR,
"ex09_plumeprops_young_plume.txt"))
plume_props_l2.save_txt(join(SAVE_DIR,
"ex09_plumeprops_aged_plume.txt"))
if SAVEFIGS:
for k in range(len(figs)):
figs[k].savefig(join(SAVE_DIR, "ex09_out_%d.%s"
% ((k + 1), FORMAT)),
format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
npt.assert_array_equal([],
[])
npt.assert_allclose(actual=[],
desired=[],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 10 - Import plume background images¶
Create a Dataset
object for a time interval containing only plume background images (on / off).
Note
Stand alone script that is not required for any of the following scripts
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 10 - Create background image dataset."""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
import pyplis
from datetime import datetime
from matplotlib.pyplot import show
# IMPORT GLOBAL SETTINGS
from SETTINGS import IMG_DIR, OPTPARSE
# Check script version
check_version()
# SCRIPT FUNCTION DEFINITIONS
def get_bg_image_lists():
"""Initialize measurement setup and creates dataset from that."""
start = datetime(2015, 9, 16, 7, 2, 0o5)
stop = datetime(2015, 9, 16, 7, 2, 30)
# Define camera (here the default ecII type is used)
cam_id = "ecII"
# the camera filter setup
filters = [pyplis.utils.Filter(type="on", acronym="F01"),
pyplis.utils.Filter(type="off", acronym="F02")]
# create camera setup
cam = pyplis.setupclasses.Camera(cam_id=cam_id, filter_list=filters)
# Create BaseSetup object (which creates the MeasGeometry object)
stp = pyplis.setupclasses.MeasSetup(IMG_DIR, start, stop, camera=cam)
ds = pyplis.dataset.Dataset(stp)
on, off = ds.get_list("on"), ds.get_list("off")
on.darkcorr_mode = True
off.darkcorr_mode = True
return on, off
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
on, off = get_bg_image_lists()
on.show_current()
off.show_current()
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
npt.assert_array_equal([],
[])
npt.assert_allclose(actual=[],
desired=[],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 11 - Image based signal dilution correction¶
This script introduces the image based signal dilution correction including automatic retrieval of terrain distances on a pixel basis.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliß (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script no. 11 - Image based signal dilution correction.
This script illustrates how extinction coefficients can be retrieved from
image data using the DilutionCorr class and by specifying suitable terrain
features in the images.
The extinction coefficients are retrieved from one on an from one off band
image recorded ~15 mins before the dataset used in the other examples. The
latter data is less suited for illustrating the feature since it contains
less terrain (therfore more sky background).
The two example images are then corrected for dilution and the results are
plotted (as comparison of the retrieved emission rate along an exemplary
plume cross section)
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
import pyplis as pyplis
from geonum import GeoPoint
from matplotlib.pyplot import show, close, subplots, Rectangle, plot
from datetime import datetime
from os.path import join, exists
from pyplis.dilutioncorr import DilutionCorr
from pyplis.doascalib import DoasCalibData
# IMPORT GLOBAL SETTINGS
from SETTINGS import IMG_DIR, SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex10_bg_imglists import get_bg_image_lists
# Check script version
check_version()
# SCRIPT OPTONS
# lower boundary for I0 value in dilution fit
I0_MIN = 0.0
AA_THRESH = 0.03
# exemplary plume cross section line for emission rate retrieval (is also used
# for full analysis in ex12)
PCS_LINE = pyplis.LineOnImage(x0=530, y0=586, x1=910, y1=200, line_id="pcs")
# Retrieval lines for dilution correction (along these lines, topographic
# distances and image radiances are determined for fitting the atmospheric
# extinction coefficients)
TOPO_LINE1 = pyplis.LineOnImage(1100, 650, 1000, 900, line_id="flank far",
color="lime",
linestyle="-")
TOPO_LINE2 = pyplis.LineOnImage(1000, 990, 1100, 990, line_id="flank close",
color="#ff33e3",
linestyle="-")
# all lines in this array are used for the analysis
USE_LINES = [TOPO_LINE1, TOPO_LINE2]
# specify pixel resolution of topographic distance retrieval (every nth pixel
# is used)
SKIP_PIX_LINES = 10
# Specify region of interest used to extract the ambient intensity (required
# for dilution correction)
AMBIENT_ROI = [1240, 10, 1300, 70]
# Specify plume velocity (for emission rate estimate)
PLUME_VELO = 4.14 # m/s (result from ex8)
# RELEVANT DIRECTORIES AND PATHS
CALIB_FILE = join(SAVE_DIR, "ex06_doascalib_aa.fts")
# SCRIPT FUNCTION DEFINITIONS
def create_dataset_dilution():
"""Create a :class:`pyplis.dataset.Dataset` object for dilution analysis.
The test dataset includes one on and one offband image which are recorded
around 6:45 UTC at lower camera elevation angle than the time series shown
in the other examples (7:06 - 7:22 UTC). Since these two images contain
more topographic features they are used to illustrate the image based
signal dilution correction.
This function sets up the measurement (geometry, camera, time stamps) for
these two images and creates a Dataset object.
"""
start = datetime(2015, 9, 16, 6, 43, 00)
stop = datetime(2015, 9, 16, 6, 47, 00)
# the camera filter setup
cam_id = "ecII"
filters = [pyplis.utils.Filter(type="on", acronym="F01"),
pyplis.utils.Filter(type="off", acronym="F02")]
geom_cam = {"lon": 15.1129,
"lat": 37.73122,
"elev": 15.0, # from field notes, will be corrected
"elev_err": 5.0,
"azim": 274.0, # from field notes, will be corrected
"azim_err": 10.0,
"focal_length": 25e-3,
"alt_offset": 7} # meters above topography
# create camera setup
cam = pyplis.setupclasses.Camera(cam_id=cam_id, filter_list=filters,
**geom_cam)
# Load default information for Etna
source = pyplis.setupclasses.Source("etna")
# Provide wind direction
wind_info = {"dir": 0.0,
"dir_err": 15.0}
# Create BaseSetup object (which creates the MeasGeometry object)
stp = pyplis.setupclasses.MeasSetup(IMG_DIR, start, stop, camera=cam,
source=source,
wind_info=wind_info)
return pyplis.dataset.Dataset(stp)
def find_view_dir(geom):
"""Perform a correction of the viewing direction using crater in img.
:param MeasGeometry geom: measurement geometry
:param str which_crater: use either "ne" (northeast) or "se" (south east)
:return: - MeasGeometry, corrected geometry
"""
# Use position of NE crater in image
posx, posy = 1051, 605 # pixel position of NE crate in image
# Geo location of NE crater (info from Google Earth)
ne_crater = GeoPoint(37.754788, 14.996673, 3287, name="NE crater")
geom.find_viewing_direction(pix_x=posx, pix_y=posy, pix_pos_err=100,
geo_point=ne_crater, draw_result=True)
return geom
def prepare_lists(dataset):
"""Prepare on and off lists for dilution analysis.
Steps:
1. get on and offband list
#. load background image list on and off (from ex10)
#. set image preparation and assign background images to on / off list
#. configure plume background model settings
:param Dataset dataset: the dilution dataset (see
:func:`create_dataset_dilution`)
:return:
- ImgList, onlist
- ImgList, offlist
"""
onlist = dataset.get_list("on")
offlist = dataset.get_list("off")
# dark_corr_mode already active
bg_onlist, bg_offlist = get_bg_image_lists()
# prepare img pre-edit
onlist.darkcorr_mode = True
onlist.gaussian_blurring = 2
offlist.darkcorr_mode = True
offlist.gaussian_blurring = 2
# prepare background images in lists (when assigning a background image
# to a ImgList, then the blurring amount of the BG image is automatically
# set to the current blurring level of the list, and if the latter is
# zero, then the BG image is blurred using filter width=1)
onlist.bg_img = bg_onlist.current_img()
offlist.bg_img = bg_offlist.current_img()
# prepare plume background modelling setup in both lists
onlist.bg_model.mode = 6
onlist.bg_model.set_missing_ref_areas(onlist.current_img())
onlist.bg_model.xgrad_line_startcol = 10
offlist.bg_model.update(**onlist.bg_model.settings_dict())
return onlist, offlist
def plot_retrieval_points_into_image(img):
"""Plot terrain distance retrieval lines into an image."""
ax = img.show(vmin=-5e16, vmax=5e18, zlabel=r"$S_{SO2}$ [cm$^{-2}$]")
ax.set_title("Retrieval lines")
for line in USE_LINES:
line.plot_line_on_grid(ax=ax, marker="", color=line.color,
lw=2, ls=line.linestyle)
for x, y, dist in dil._add_points:
plot(x, y, " or")
return ax
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
if not exists(CALIB_FILE):
raise IOError("Calibration file could not be found at specified "
"location:\n %s\nYou might need to run example 6 first")
close("all")
pcs_line = PCS_LINE
calib = DoasCalibData()
calib.load_from_fits(CALIB_FILE)
# create dataset and correct viewing direction
ds = create_dataset_dilution()
geom = find_view_dir(ds.meas_geometry)
# get plume distance image
pix_dists, _, plume_dists = geom.compute_all_integration_step_lengths()
# Create dilution correction class
dil = DilutionCorr(USE_LINES, geom, skip_pix=SKIP_PIX_LINES)
# you can also manually add a single pixel position in the image that is
# used to retrieve topographic distances (this function allows you also to
# set the distance to the feature manually, since sometimes the SRTM data-
# set is incomplete or has too low resolution)
dil.add_retrieval_point(700, 930)
dil.add_retrieval_point(730, 607)
# Determine distances to the two lines defined above (every 6th pixel)
for line_id in dil.line_ids:
dil.det_topo_dists_line(line_id)
# Plot the results in a 3D map
basemap = dil.plot_distances_3d(alt_offset_m=10, axis_off=False)
# retrieve pixel distances for pixels on the line
# (for emission rate estimate)
pix_dists_line = pcs_line.get_line_profile(pix_dists)
# get pixel coordinates of PCS center position ...
col, row = pcs_line.center_pix
# ... and get uncertainty in plume distance estimate for the column
pix_dist_err = geom.pix_dist_err(col)
# Prepare on and off-band list for retrieval of extinction coefficients
onlist, offlist = prepare_lists(ds)
# Activate vignetting correction in both lists (required to extract
# measured intensities along the terrain features)
onlist.vigncorr_mode = True
offlist.vigncorr_mode = True
# get the vignetting corrected images
on_vigncorr = onlist.current_img()
off_vigncorr = offlist.current_img()
# estimate ambient intensity for both filters
ia_on = on_vigncorr.crop(AMBIENT_ROI, True).mean()
ia_off = off_vigncorr.crop(AMBIENT_ROI, True).mean()
# perform dilution anlysis and retrieve extinction coefficients (on-band)
ext_on, _, _, ax0 = dil.apply_dilution_fit(img=on_vigncorr,
rad_ambient=ia_on,
i0_min=I0_MIN,
plot=True)
ax0.set_ylabel("Terrain radiances (on band)", fontsize=14)
ax0.set_ylim([0, 2500])
# perform dilution anlysis and retrieve extinction coefficients (off-band)
ext_off, i0_off, _, ax1 = dil.apply_dilution_fit(img=off_vigncorr,
rad_ambient=ia_off,
i0_min=I0_MIN,
plot=True)
ax1.set_ylabel("Terrain radiances (off band)", fontsize=14)
ax1.set_ylim([0, 2500])
# determine plume pixel mask from AA image
onlist.aa_mode = True
plume_pix_mask = onlist.get_thresh_mask(AA_THRESH)
plume_pix_mask[840:, :] = 0 # remove tree in lower part of the image
onlist.aa_mode = False
# assign the just retrieved extinction coefficients to the respective
# image lists
onlist.ext_coeffs = ext_on
offlist.ext_coeffs = ext_off
# save the extinction coefficients into a txt file (re-used in example
# script 12). They are stored as pandas.Series object in the ImgList
onlist.ext_coeffs.to_csv(join(SAVE_DIR, "ex11_ext_scat_on.txt"))
offlist.ext_coeffs.to_csv(join(SAVE_DIR, "ex11_ext_scat_off.txt"))
# now activate automatic dilution correction in both lists
# get dilution corrected on and off-band image
onlist.dilcorr_mode = True
offlist.dilcorr_mode = True
# get current dilution corrected raw images (i.e. in intensity space)
on_corr = onlist.this
off_corr = offlist.this
# now activate tau mode (note that dilution correction mode is still
# active)
onlist.tau_mode = True
offlist.tau_mode = True
# extract tau images
tau_on_corr = onlist.this
tau_off_corr = offlist.this
# determine corrected SO2-CD image from the image lists
so2_img_corr = calib(tau_on_corr - tau_off_corr)
so2_img_corr.edit_log["is_tau"] = True # for plotting
so2_cds_corr = pcs_line.get_line_profile(so2_img_corr)
(phi_corr,
phi_corr_err) = pyplis.fluxcalc.det_emission_rate(
cds=so2_cds_corr,
velo=PLUME_VELO,
pix_dists=pix_dists_line,
cds_err=calib.err(),
pix_dists_err=pix_dist_err)
# determine uncorrected so2-CD image from the image lists
offlist.dilcorr_mode = False
onlist.dilcorr_mode = False
# the "this" attribute returns the current list image (same as
# method "current_img()")
so2_img_uncorr = calib(onlist.this - offlist.this)
# Retrieve column density profile along PCS in uncorrected image
so2_cds_uncorr = pcs_line.get_line_profile(so2_img_uncorr)
# Calculate flux and uncertainty
(phi_uncorr,
phi_uncorr_err) = pyplis.fluxcalc.det_emission_rate(
cds=so2_cds_uncorr,
velo=PLUME_VELO,
pix_dists=pix_dists_line,
cds_err=calib.err(),
pix_dists_err=pix_dist_err)
# IMPORTANT STUFF FINISHED (below follow some plots)
ax2 = plot_retrieval_points_into_image(so2_img_corr)
pcs_line.plot_line_on_grid(ax=ax2, ls="-", color="g")
ax2.legend(loc="best", framealpha=0.5, fancybox=True, fontsize=20)
ax2.set_title("Dilution corrected AA image", fontsize=12)
ax2.get_xaxis().set_ticks([])
ax2.get_yaxis().set_ticks([])
x0, y0, w, h = pyplis.helpers.roi2rect(AMBIENT_ROI)
ax2.add_patch(Rectangle((x0, y0), w, h, fc="none", ec="c"))
fig, ax3 = subplots(1, 1)
ax3.plot(so2_cds_uncorr, ls="-", color="#ff33e3",
label=r"Uncorr: $\Phi_{SO2}=$%.2f (+/- %.2f) kg/s"
% (phi_uncorr / 1000.0, phi_uncorr_err / 1000.0))
ax3.plot(so2_cds_corr, "-g", lw=3,
label=r"Corr: $\Phi_{SO2}=$%.2f (+/- %.2f) kg/s"
% (phi_corr / 1000.0, phi_corr_err / 1000.0))
ax3.set_title("Cross section profile", fontsize=12)
ax3.legend(loc="best", framealpha=0.5, fancybox=True, fontsize=12)
ax3.set_xlim([0, len(pix_dists_line)])
ax3.set_ylim([0, 5e18])
ax3.set_ylabel(r"$S_{SO2}$ [cm$^{-2}$]", fontsize=14)
ax3.set_xlabel("PCS", fontsize=14)
ax3.grid()
# also plot plume pixel mask
ax4 = pyplis.Img(plume_pix_mask).show(cmap="gray", tit="Plume pixel mask")
if SAVEFIGS:
ax = [ax0, ax1, ax2, ax3, ax4]
for k in range(len(ax)):
ax[k].set_title("") # remove titles for saving
ax[k].figure.savefig(join(SAVE_DIR, "ex11_out_%d.%s"
% (k, FORMAT)),
format=FORMAT, dpi=DPI)
basemap.ax.set_axis_off()
basemap.ax.view_init(15, 345)
basemap.ax.figure.savefig(join(SAVE_DIR, "ex11_out_5.%s" % FORMAT),
format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
npt.assert_array_equal([],
[])
npt.assert_allclose(actual=[],
desired=[],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")
Example 12 - Emission rate analysis (Etna example data)¶
Perform emission rate analysis for the example data. The analysis is performed along one plume cross section (in the image center) and using three different plume velocity retrievals.
Code
# -*- coding: utf-8 -*-
#
# Pyplis is a Python library for the analysis of UV SO2 camera data
# Copyright (C) 2017 Jonas Gliss (jonasgliss@gmail.com)
#
# This program is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License a
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
"""Pyplis example script 12 - Etna emission rate retrieval.
This example import results from the previous examples, for instance the AA
image list including measurement geometry (ex 4), the DOAS calibration
information (which was stored as FITS file, see ex. 6) and the AA sensitivity
correction mask retrieved from the cell calibration and normalised to the
position of the DOAS FOV (ex 7). The emission rates are retrieved for three
different plume velocity retrievals: 1. using the global velocity vector
retrieved from the cross correlation algorithm (ex8), 2. using the raw output
of the optical flow Farneback algorithm (``flow_raw``) and 3. using the
histogram based post analysis of the optical flow field (``flow_histo``).
The analysis is performed using the EmissionRateAnalysis class which basically
checks the AA list and activates ``calib_mode`` (-> images are loaded as
calibrated gas CD images) and loops over all images to retrieve the emission
rates for the 3 velocity modes. Here, emission rates are retrieved along 1
exemplary plume cross section. This can be easily extended by adding additional
PCS lines in the EmissionRateAnalysis class using ``add_pcs_line``.
The results for each velocity mode and for each PCS line are stored within
EmissionRateResults classes.
"""
from __future__ import (absolute_import, division)
from SETTINGS import check_version
from os.path import join, exists
from matplotlib.pyplot import close, show, GridSpec, figure, rc_context
from matplotlib.cm import get_cmap
import pyplis
# IMPORT GLOBAL SETTINGS
from SETTINGS import SAVEFIGS, SAVE_DIR, FORMAT, DPI, OPTPARSE, LINES
# IMPORTS FROM OTHER EXAMPLE SCRIPTS
from ex04_prep_aa_imglist import prepare_aa_image_list
rc_context({'font.size': '18'})
# Check script version
check_version()
PCS = LINES[0]
# If false, then only the working environment is initalised
DO_EVAL = True
# Dilution correction
DILCORR = True
# You can specify here if you only want a certain number of images analysed
START_INDEX = 0
STOP_INDEX = None # 20
# SCRIPT OPTONS
PYRLEVEL = 1
PLUME_VELO_GLOB = 4.29 # m/s
PLUME_VELO_GLOB_ERR = 1.5
# applies multi gauss fit to retrieve local predominant displacement
# direction, if False, then the latter is calculated from 1. and 2. moment
# of histogram (Faster but more sensitive to additional peaks in histogram)
HISTO_ANALYSIS_MULTIGAUSS = True
# molar mass of SO2
MMOL = 64.0638 # g/mol
# minimum required SO2-CD for emission-rate retrieval
CD_MIN = 5e16
# activate background check mode, if True, emission rates are only
# retrieved for images showing SO2-CDs within specified interval around
# zero in BG reference rectangle LOG_ROI_SKY (see above). This can be
# used to ensure that significant systematic errors are induced in case
# the plume background retrieval failed. The latter could, for instance
# happen, if, for instance a cloud moves through one of the background
# reference areas used to model the background (cf. example script 3)
REF_CHECK_LOWER = -5e16
REF_CHECK_UPPER = 5e16
REF_CHECK_MODE = True
# the following ROI is in the upper right image corner, where no gas occurs in
# the time series. It is used to log mean, min and max for each analysed image
# this information can be used to check, whether the plume background retrieval
# worked well
LOG_ROI_SKY = [530, 30, 600, 100] # correspond to pyrlevel 1
# RELEVANT DIRECTORIES AND PATHS
# DOAS calibration results from example script 6
CALIB_FILE = join(SAVE_DIR, "ex06_doascalib_aa.fts")
# Scattering extinction coeffcients from example script 11 (stored as txt)
EXT_ON = join(SAVE_DIR, "ex11_ext_scat_on.txt")
EXT_OFF = join(SAVE_DIR, "ex11_ext_scat_off.txt")
# AA sensitivity correction mask retrieved from cell calib in script 7
CORR_MASK_FILE = join(SAVE_DIR, "ex07_aa_corr_mask.fts")
# time series of predominant displacement vector from histogram analysis of
# optical flow field in ROI around the PCS line "young_plume" which is used
# here for the emission rate retrieval. These information is optional, and is
# calculated during the evaluation if not provided
RESULT_PLUMEPROPS_HISTO = join(SAVE_DIR, "ex09_plumeprops_young_plume.txt")
# SCRIPT FUNCTION DEFINITIONS
def plot_and_save_results(ana, line_id="young_plume", date_fmt="%H:%M"):
# plot colors for different optical flow retrievals
cmap = get_cmap("Oranges")
c_optflow_hybrid = cmap(255)
c_optflow_histo = cmap(175)
c_optflow_raw = cmap(100)
fig = figure(figsize=(16, 12))
gs = GridSpec(4, 1, height_ratios=[.6, .2, .2, .2], hspace=0.05)
ax3 = fig.add_subplot(gs[3])
ax0 = fig.add_subplot(gs[0], sharex=ax3)
ax1 = fig.add_subplot(gs[1], sharex=ax3)
ax2 = fig.add_subplot(gs[2], sharex=ax3)
ax1.yaxis.tick_right()
ax1.yaxis.set_label_position("right")
ax3.yaxis.tick_right()
ax3.yaxis.set_label_position("right")
# Get emission rate results for the PCS line
res0 = ana.get_results(line_id=line_id, velo_mode="glob")
res1 = ana.get_results(line_id=line_id, velo_mode="flow_raw")
res2 = ana.get_results(line_id=line_id, velo_mode="flow_histo")
res3 = ana.get_results(line_id=line_id, velo_mode="flow_hybrid")
res0.save_txt(join(SAVE_DIR, "ex12_flux_velo_glob.txt"))
res1.save_txt(join(SAVE_DIR, "ex12_flux_flow_raw.txt"))
res2.save_txt(join(SAVE_DIR, "ex12_flux_flow_histo.txt"))
res3.save_txt(join(SAVE_DIR, "ex12_flux_flow_hybrid.txt"))
# Plot emission rates for the different plume speed retrievals
res0.plot(yerr=True, date_fmt=date_fmt, ls="-", ax=ax0,
color="c", ymin=0, alpha_err=0.08)
res1.plot(yerr=False, ax=ax0, ls="-", color=c_optflow_raw, ymin=0)
res2.plot(yerr=False, ax=ax0, ls="--", color=c_optflow_histo, ymin=0)
res3.plot(yerr=True, ax=ax0, lw=3, ls="-", color=c_optflow_hybrid, ymin=0)
# ax[0].set_title("Retrieved emission rates")
ax0.legend(loc='best', fancybox=True, framealpha=0.5, fontsize=12)
ax0.grid()
# Plot effective velocity retrieved from optical flow histogram analysis
res3.plot_velo_eff(ax=ax1, date_fmt=date_fmt, color=c_optflow_hybrid)
# ax[1].set_title("Effective plume speed
# (from optflow histogram analysis)")
ax1.set_ylim([0, ax1.get_ylim()[1]])
# Plot time series of predominant plume direction (retrieved from optical
# flow histogram analysis and stored in object of type LocalPlumeProperties
# which is part of plumespeed.py module
ana.pcs_lines[line_id].plume_props.plot_directions(ax=ax2,
date_fmt=date_fmt,
color=c_optflow_hybrid)
ax2.set_ylim([-180, 180])
pyplis.helpers.rotate_xtick_labels(ax=ax2)
ax0.set_xticklabels([])
ax1.set_xticklabels([])
ax2.set_xticklabels([])
# tight_layout()
ax3 = ana.plot_bg_roi_vals(ax=ax3, date_fmt="%H:%M")
# gs.tight_layout(fig, h_pad=0)#0.03)
gs.update(hspace=0.05, top=0.97, bottom=0.07)
return fig
# SCRIPT MAIN FUNCTION
if __name__ == "__main__":
close("all")
figs = []
if not exists(CALIB_FILE):
raise IOError("Calibration file could not be found at specified "
"location:\n%s\nPlease run example 6 first")
if not exists(CORR_MASK_FILE):
raise IOError("Cannot find AA correction mask, please run example"
"script 7 first")
# convert the retrieval line to the specified pyramid level (script option)
pcs = PCS.convert(to_pyrlevel=PYRLEVEL)
# now try to load results of optical flow histogram analysis performed for
# this line in script no. 9. and assign them to the pcs line. This has the
# advantage, that missing velocity vectors (i.e. from images where optical
# flow analysis failed) can be interpolated. It is, however, not
# necessarily required to do this in advance. In the latter case the
# emission rates show gaps at all images, where the optical flow was
# considered not reliable
try:
p = pyplis.LocalPlumeProperties()
p.load_txt(RESULT_PLUMEPROPS_HISTO)
p = p.to_pyrlevel(PYRLEVEL)
fig = p.plot(color="r")
# p.interpolate()
# p = p.apply_significance_thresh(0.2).interpolate()
# p = p.apply_median_filter(3).apply_gauss_filter(2)
fig = p.plot(date_fmt="%H:%M", fig=fig)
pcs.plume_props = p
except BaseException:
print("Local plume properties could not be loaded and will be "
"calculated during the emission rate analysis")
# Load AA list
# includes viewing direction corrected geometry
aa_list = prepare_aa_image_list()
aa_list.pyrlevel = PYRLEVEL
if DILCORR:
aa_list.import_ext_coeffs_csv(EXT_ON)
aa_list.get_off_list().import_ext_coeffs_csv(EXT_OFF)
# Load DOAS calbration data and FOV information (see example 6)
doascalib = pyplis.doascalib.DoasCalibData()
doascalib.load_from_fits(file_path=CALIB_FILE)
doascalib.fit_calib_data()
# Load AA corr mask and set in image list(is normalised to DOAS FOV see
# ex7)
aa_corr_mask = pyplis.Img(CORR_MASK_FILE)
aa_list.senscorr_mask = aa_corr_mask
# set DOAS calibration data in image list
aa_list.calib_data = doascalib
ana = pyplis.EmissionRateAnalysis(
imglist=aa_list,
bg_roi=LOG_ROI_SKY,
pcs_lines=pcs,
velo_glob=PLUME_VELO_GLOB,
velo_glob_err=PLUME_VELO_GLOB_ERR,
ref_check_lower_lim=REF_CHECK_LOWER,
ref_check_upper_lim=REF_CHECK_UPPER,
velo_dir_multigauss=HISTO_ANALYSIS_MULTIGAUSS,
senscorr=True,
dilcorr=DILCORR)
ana.settings.ref_check_mode = REF_CHECK_MODE
ana.settings.velo_modes["flow_raw"] = 1
ana.settings.velo_modes["flow_histo"] = True
ana.settings.velo_modes["flow_hybrid"] = 1
ana.settings.min_cd = CD_MIN
# plot all current PCS lines into current list image (feel free to define
# and add more PCS lines above)
ax = ana.plot_pcs_lines(
vmin=-
5e18,
vmax=6e18,
tit="Dilution corr: %s" %
DILCORR)
ax = ana.plot_bg_roi_rect(ax=ax, to_pyrlevel=PYRLEVEL)
figs.append(ax.figure)
if not DO_EVAL:
aa_list.dilcorr_mode = not DILCORR
aa_list.show_current(
vmin=-
5e18,
vmax=6e18,
tit="Dilution corr: %s" %
(not DILCORR))
# you can check the settings first
print(ana.settings)
# check if optical flow works
ana.imglist.optflow_mode = True
aa_mask = ana.imglist.get_thresh_mask(CD_MIN)
ana.imglist.optflow.plot_flow_histograms(line=pcs, pix_mask=aa_mask)
else:
ana.run_retrieval(start_index=START_INDEX,
stop_index=STOP_INDEX)
figs.append(plot_and_save_results(ana))
# the EmissionRateResults class has an informative string
# representation
print(ana.get_results("young_plume", "flow_histo"))
if SAVEFIGS:
for k in range(len(figs)):
figs[k].savefig(join(SAVE_DIR, "ex12_out_%d.%s" % (k + 1, FORMAT)),
format=FORMAT, dpi=DPI)
# IMPORTANT STUFF FINISHED (Below follow tests and display options)
# Import script options
(options, args) = OPTPARSE.parse_args()
# If applicable, do some tests. This is done only if TESTMODE is active:
# testmode can be activated globally (see SETTINGS.py) or can also be
# activated from the command line when executing the script using the
# option --test 1
if int(options.test):
import numpy.testing as npt
from os.path import basename
npt.assert_array_equal([],
[])
npt.assert_allclose(actual=[],
desired=[],
rtol=1e-7)
print("All tests passed in script: %s" % basename(__file__))
try:
if int(options.show) == 1:
show()
except BaseException:
print("Use option --show 1 if you want the plots to be displayed")