Title

This notebook is an AC SAF notebook prepared by BIRA-IASB (Gaia Pinardi, Hugues Brenot), on illustrating the use of volcanic SO2 GOME-2 data:

• reading and plotting of the AC SAF GOME-2 L2 data from the GDP 4.8/4.9 data using the xarray tools
• this notebook focus on SO2 VCD over one orbit that passes either on the Norilsk and Kivu region and plot and calculate the SO2 mass around it.

A sample zip file is provided with this notebook and instructions on how to register and download more data from the AC SAF archive are also provided below.

AC SAF recomend the usage of HARP for reading, filtering, and gridding, as shown in a companion notebook (also available here: https://acsaf.org/notebooks.php).

Table of Contents

This notebook contains several sections:

• Data description, with documentation, downloading informations and data file format
• The different imports
• Definition of a number of useful variables/region bounding boxes/colorbar and some functions
• Start of the routine itself
• Selection of a specific case to look to
• Illustration on how to read the GOME-2 SO2 data
• Including filter the GOME-2 SO2 data (flag, which SO2 column, backscans, plot region)
• Then the plot of these data with Cartopy and the calculation of the SO2 mass

Data description

Documentation

The GOME-2/MetOp total column products are generated at DLR in the framework of EUMETSAT’s AC SAF.

For information and access to all AC SAF products, please refer to the AC SAF web page and help desk: AC SAF Web page: http://acsaf.org AC SAF Helpdesk: helpdesk@acsaf.org

Additional information can be found in EUMETSAT's product Navigator accessible through the following link: http://navigator.eumetsat.int

Current and historical quick look images as well as value added products generated by DLR are available at: http://atmos.eoc.dlr.de/gome2

The ATBD (Algorithm Theoretical Basis Documents) and PUM (PRODUCT USER MANUAL) are available here: https://acsaf.org/docs/atbd/Algorithm_Theoretical_Basis_Document_NTO_OTO_Nov_2019.pdf and https://acsaf.org/docs/pum/Product_User_Manual_NTO_OTO_Nov_2019.pdf

Data download is described here: https://acsaf.org/offline_access.php and rely on accessing an ftp server.

Downloading test data

For this example, the data are included in the zip folder available with this notebook: SO2.zip

save then the zip folder into the data folder, at the same level than this notebook

Data files

The GOME-2 GDP 4.8/4.9 L2 reprocessed datafiles are structures as one file per orbit, ie 15 files per day:

HDF5 Product Format Description

The product content comprises the following groups:

Several gases are provided into the same orbital file, as mentionned in the file name:

i.e., GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201134102_051_METOPA_62878_DLR_06.HDF5 will include O3, NO2, NO2Tropo, BrO, SO2, H2O and HCHO data.

In the TOTAL_COLUMNS, the final VCD values are given for all the gases, in a order specified in the file name and in the variable /META_DATA/MainSpecies:

In addition, each gas has his own specific variables within the DETAILED_RESULTS group.

SO2 specific information

GOME-2 reports SO2 VCDs under several plume-height hypotheses, since the retrieved column depends strongly on this a-priori assumption. The 4 standard assumed layer heights are 15 km, 6 km, 2 km and 1 km.

The highest altitudes are usually relevant for the volcanic cases, while the lower altitudes are targetting anthropogenic emissions. The recommendation is to interpolate between the SO2 VCD provided at the different altitudes, for the altitude that best match the situation of interest.

In order to be able to do this, we recommend to read the VCD in the DETAILED_RESULTS/SO2 group:

The order of the altitudes is provided in the PlumeHeights variable in the same SO2 group.

Imports

Standard library / numerics / IO imports for what follows: h5py/netCDF4/xarray to read the HDF5 file with three different approaches, numpy for arrays, cf_units for unit handling, datetime/warnings as usual helpers (glob for multi-file patterns, matplotlib/cartopy for the maps built later).

These imports support several parallel approaches shown later in the notebook: xarray/h5py to read the HDF5 file directly, cartopy/matplotlib (with the PolyCollection/GeoAxes patch) to draw the satellite pixels as true geolocated polygons (not just points), numpy.ma and scipy.ndimage.label for the noise/background filtering of the SO2 signal, and math/calendar/datetime for the date and solar-angle calculations used to build the SZA quality flag. warnings.filterwarnings("ignore") silences library warnings (convenient here, but worth removing if you need to debug something upstream).

#
# pip install cartopy
import sys
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from matplotlib.axes import Axes
from matplotlib.collections import PolyCollection
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from math import radians, cos, sin, asin, sqrt, acos, atan2
from numpy import ma
import calendar
import h5py
import glob 
from datetime import datetime
from scipy.ndimage import label
import warnings
warnings.filterwarnings("ignore")

Define some useful variables

constants/unit conversions

In order to later turn an SO2 column (in DU) into an SO2 mass in kilotonnes once combined with a pixel area, we need a few physical conversion factors: Dobson Units (DU) to molecules/m², to mass (kt, kg/m²), using Avogadro's number and the molar mass of SO2 (64.066 g/mol). We provide the unit conversions here

######## constants
DU_to_m2 = 2.68668e20     #  molecule/m2
Avogadro = 6.022140857e23 #  molecule/mol
molSO2 = 64.066           #  g/mol
DU_to_kT = 1e-9 * DU_to_m2 * molSO2 / Avogadro  # DU to mass kT/m2
molecperm2_to_DU = 2241.3962 # conversion_factor_DU
cm2tom2 = 1e4
DU_to_GperM2 = DU_to_m2 * molSO2 / Avogadro    # 0.02863
molecperm2_to_GperM2 = molSO2 / Avogadro

fillValue = -999.

day1950 = 2433282.5
molecPERcm2_to_DU = 2.68668e16

DU_TO_KG_M2 = 2.6867e20 / 6.022e23 * 64.066 / 1000  # ~ 2.85e-5 kg/m² per DU

nice colorbar for SO2

We build an SO2-specific colour scale (white → violet → cyan → green → yellow → orange → red → dark red) with a non-linear, piecewise normalization (PiecewiseNorm). The colorbar bounds (bounds = [0, 0.5, 1, 1.5, 2, 3, 4, 5, 6, 8, 10]) are irregular, so a standard linear normalization would crowd out the low values; PiecewiseNorm stretches every interval to the same visual width regardless of its actual size in DU, so both the background noise and concentrated plumes stay visible on the same map.

###### colorbar 1
colors = [
    "#FFFFFF","#fcecfd","#f1bcfd",
    "#8ba3fa","#27f8fb","#34fb27",
    "#ffe94a","#ffa24a","#fe0001",
    "#8d0101","#5d0003"
]
###### colorbar continous linear
# # Interpolate discrete colors into a smooth gradient
# cmap_so2_smooth = mcolors.LinearSegmentedColormap.from_list("SO2_smooth", colors, N=256)
# # No boundary norm needed for smooth
# norm_so2_smooth = None
# cmap_so2 = cmap_so2_smooth
# norm_so2 = norm_so2_smooth

###### colorbar continous non-linear
bounds = [0.0,0.5,1,1.5,2,3,4,5,6,8,10]

# --- Continuous non-linear normalization ---
class PiecewiseNorm(mcolors.Normalize):
    """
    Normalize data to 0..1 according to piecewise bounds.
    Each interval is mapped linearly but interval lengths are proportional to their bounds.
    """
    def __init__(self, boundaries, clip=False):
        super().__init__(vmin=boundaries[0], vmax=boundaries[-1], clip=clip)
        self.boundaries = np.array(boundaries)
        # Create cumulative lengths of intervals
        self._cum_lengths = np.linspace(0, 1, len(boundaries))

    def __call__(self, value, clip=None):
        value = np.asarray(value)
        result = np.interp(value, self.boundaries, self._cum_lengths)
        return np.ma.masked_array(result)

# Colormap
cmap_so2 = mcolors.LinearSegmentedColormap.from_list("SO2_cont", colors, N=256)
norm_so2 = PiecewiseNorm(bounds)

define some functions

We define here some functions that we will use to calculate SO2 mass, visualize footprints, calculate date from julian date, ...

Utility functions used further down:

• SZAlimitGOME2 — for a given day-of-year, computes the acceptable solar zenith angle (SZA) limit for the northern and southern hemisphere, interpolated between seasons, because NO2/SO2 retrieval quality degrades close to the terminator (high SZA).
• haversine — great-circle distance (m) between two lat/lon points.
• triangle_area / pixel_area — area of a satellite pixel (a quadrilateral defined by its 4 corners), computed by splitting it into 2 triangles and applying Heron's formula, since a GOME-2 ground pixel is not a regular lat/lon rectangle.
• calcSO2massPXL — converts an SO2 column (DU) and a pixel area (m²) into an SO2 mass (kt) using DU_TO_KG_M2; clips negative values (not physically meaningful).
• visualize_footprints — the main mapping function: draws each satellite pixel as a true polygon (via PolyCollection), not a scatter point; auto-selects a polar projection near the poles; adds coastlines/borders; applies the colorbar bounds defined above; builds a dynamic title (SO2 mass, area, max value and its location); and plot the figure.

####### define some useful functions

def SZAlimitGOME2(year, DOY):
    """ 
    SZAnorth, SZAsouth = SZAlimitGOME2(year, DOY)
    compute, for a given day-of-year, the acceptable solar zenith angle (SZA) limit for the northern and southern hemisphere, 
    interpolated between seasons, because NO2/SO2 retrieval quality degrades close to the terminator (high SZA).
    """
    if year % 4 == 0:
        NBdays = 366
        DOYspring =  81.
        DOYsummer = 173.
        DOYautumn = 265.
        DOYwinter = 356.
    else:
        NBdays = 365.
        DOYspring =  80.
        DOYsummer = 172.
        DOYautumn = 264.
        DOYwinter = 355.

    #   North SZA limit estimate
    SZAnorthSPRING = 75.
    SZAnorthSUMMER = 70.
    SZAnorthAUTUMN = 75.
    SZAnorthWINTER = 75.
    SZAnorth = 75.

    if DOY >= DOYspring and DOY <= DOYsummer:
        SZAnorth = SZAnorthSPRING + (SZAnorthSUMMER - SZAnorthSPRING) * (DOY - DOYspring) / (DOYsummer -  DOYspring)

    if DOY > DOYsummer and DOY <= DOYautumn:
        SZAnorth = SZAnorthSUMMER + (SZAnorthAUTUMN - SZAnorthSUMMER) * (DOY - DOYsummer) / (DOYautumn -  DOYsummer)

    #   South SZA limit estimate
    SZAsouthSPRING = 80.
    SZAsouthSUMMER = 75.
    SZAsouthAUTUMN = 80.
    SZAsouthWINTER = 70.

    if DOY < DOYspring:
        SZAsouth = SZAsouthWINTER + (SZAsouthSPRING - SZAsouthWINTER) * ((NBdays - DOYwinter) + DOY) / ((NBdays - DOYwinter) + DOYspring)

    if DOY >= DOYspring and DOY <= DOYsummer:
        SZAsouth = SZAsouthSPRING + (SZAsouthSUMMER - SZAsouthSPRING) * (DOY - DOYspring) / (DOYsummer -  DOYspring)

    if DOY > DOYsummer and DOY <= DOYautumn:
        SZAsouth = SZAsouthSUMMER + (SZAsouthAUTUMN - SZAsouthSUMMER) * (DOY - DOYsummer) / (DOYautumn -  DOYsummer)

    if DOY > DOYautumn and DOY <= DOYwinter:
        SZAsouth = SZAsouthAUTUMN + (SZAsouthWINTER - SZAsouthAUTUMN) * (DOY - DOYautumn) / (DOYwinter -  DOYautumn)

    if DOY > DOYwinter:
        SZAsouth = SZAsouthWINTER + (SZAsouthSPRING - SZAsouthWINTER) * (NBdays - DOY) / ((NBdays - DOYwinter) + DOYspring)

    return SZAnorth, SZAsouth


def haversine(lat1, lon1, lat2, lon2):
    """ dist=haversine(lat1, lon1, lat2, lon2)
    calculate the great-circle distance (m) between two lat/lon points.
    """
    R = 6371000.0  # meters
    lat1 = np.radians(lat1)
    lat2 = np.radians(lat2)
    dlat = lat2 - lat1
    dlon = np.radians(lon2 - lon1)
    a = np.sin(dlat/2)**2 + np.cos(lat1)*np.cos(lat2)*np.sin(dlon/2)**2
    c = 2 * np.arcsin(np.sqrt(a))
    return R * c

def triangle_area(a, b, c):
    """ area=triangle_area(a, b, c))
    calculate the area of a triangle.
    """
    s = 0.5 * (a + b + c)
    return np.sqrt(np.maximum(s * (s-a) * (s-b) * (s-c), 0))

def pixel_area(lonNW, lonNE, lonSE, lonSW, latNW, latNE, latSE, latSW):
    """ areaPix=pixel_area(lonNW, lonNE, lonSE, lonSW, latNW, latNE, latSE, latSW)
    calculate the area of a satellite pixel (a quadrilateral defined by its 4 corners), computed by splitting it into 2 triangles 
    and applying Heron's formula, since a GOME-2 ground pixel is not a regular lat/lon rectangle.
    """
    # triangle 1: NW-NE-SE
    a1 = haversine(latNW, lonNW, latNE, lonNE)
    b1 = haversine(latNE, lonNE, latSE, lonSE)
    c1 = haversine(latSE, lonSE, latNW, lonNW)
    area1 = triangle_area(a1, b1, c1)
    # triangle 2: NW-SE-SW
    a2 = haversine(latNW, lonNW, latSE, lonSE)
    b2 = haversine(latSE, lonSE, latSW, lonSW)
    c2 = haversine(latSW, lonSW, latNW, lonNW)
    area2 = triangle_area(a2, b2, c2)
    return area1 + area2   # m²

def calcSO2massPXL(so2DU, lonNW, lonNE, lonSE, lonSW, latNW, latNE, latSE, latSW):
    """ mass_kt, area_km2 = calcSO2massPXL(so2DU, lonNW, lonNE, lonSE, lonSW, latNW, latNE, latSE, latSW)
    converts an SO2 column (DU) and a pixel area (m²) into an SO2 mass (kt) using `DU_TO_KG_M2`; 
    clips negative values (not physically meaningful).
    """
    area = pixel_area(lonNW, lonNE, lonSE, lonSW, latNW, latNE, latSE, latSW)
    DU_TO_KG_M2 = 2.85e-5
    mass_kg = so2DU * area * DU_TO_KG_M2
    mass_kt = mass_kg * 1e-6
    # clean negatives
    mass_kt = np.where(mass_kt < 0, 0, mass_kt)
    area = np.where(area < 0, 0, area)
    return mass_kt, area * 1e-6  # km²

def visualize_footprints(data_array,
                        loncor,
                        latcor,
                        period,
                        projection=ccrs.PlateCarree(),
                        color_scale='viridis',
                        cmap=None,
                        norm=None,
                        unit='',
                        long_name='',
                        vmin=None,
                        vmax=None,
                        set_global=False,
                        latmin=None,
                        latmax=None,
                        lonmin=None,
                        lonmax=None,
                        title_pad=25,
                        coastline_res='10m',
                        metop=None,
                        orbit=None,
                        SO2mass=None,
                        SO2area=None,
                        SO2max=None,
                        SO2max_lat=None,
                        SO2max_lon=None):
    """
    main mapping function, that draws each satellite pixel as a true polygon.
    """

    # --- Ensure numpy arrays ---
    data_array = np.asarray(data_array).flatten()
    loncor = np.asarray(loncor)
    latcor = np.asarray(latcor)

    # --- Decide projection automatically ---
    if projection is None:
        if latmin is not None:
            if latmin >= 60:
                projection = ccrs.NorthPolarStereo()
            elif latmax is not None and latmax <= -60:
                projection = ccrs.SouthPolarStereo()
            else:
                projection = ccrs.PlateCarree()
        else:
            projection = ccrs.PlateCarree()

    # --- Build polygons ---
    polygons = [list(zip(loncor[i], latcor[i])) for i in range(len(data_array))]

    # --- Plot ---
    fig = plt.figure(figsize=(10, 5))
    ax = plt.axes(projection=projection)

    # Gridlines
    gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.7, linestyle='--')
    gl.top_labels = False
    gl.right_labels = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER

    # Coastlines
    ax.add_feature(cfeature.OCEAN, facecolor='lightgrey', zorder=0)

    # Land (countries)
    ax.add_feature(cfeature.LAND, facecolor='darkgrey', zorder=0)

    # Optional: borders and coastlines on top
    ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='black')
    ax.coastlines(resolution=coastline_res, linewidth=0.8, color='black')


    # Polygons
    poly = PolyCollection(
        polygons,
        array=data_array,
        cmap=cmap if cmap is not None else color_scale,
        norm=norm,
        edgecolors='none',
        transform=ccrs.PlateCarree()
    )
    if norm is None:
        poly.set_clim(vmin, vmax)

    ax.add_collection(poly)

    # Colorbar
    # bounds = [0.0,0.5,1,1.5,2,3,4,5,10,20,50]
    bounds = [0.0,0.5,1,1.5,2,3,4,5,6,8,10]
    # bounds = [0,5,10,20,50]

    # Colorbar with ticks exactly at bounds
    cbar = plt.colorbar(poly, ax=ax, fraction=0.046, pad=0.08, ticks=bounds)
    cbar.set_label(unit)
    cbar.ax.set_yticklabels([str(b) for b in bounds])

    # Build title
    title_str = long_name
    # Format SO2
    title_str = title_str.replace("SO2", r"$\mathrm{SO_2}$")
    # Add METOP and orbit if provided
    if metop is not None:
        title_str += f" - GOME2-{metop}"
    title_str += f" - {period[0:4]}/{period[4:6]}/{period[6:8]}"
    if orbit is not None:
        title_str += f" - Orbit {orbit}"
    # Add SO2 mass and area
    if SO2mass is not None:
        title_str += f" - Mass: {SO2mass:.2f} kt"
    if SO2area is not None:
        title_str += f", Area: {SO2area:.1f} km²"
    # Add SO2 max at coordinates
    if SO2max is not None and SO2max_lat is not None and SO2max_lon is not None:
        title_str += f" - Max: {SO2max:.2f} DU at ({SO2max_lat:.2f}°E°, {SO2max_lon:.2f}°N°)"
    # Set the title
    # ax.set_title(title_str, fontsize=12, pad=title_pad)
    ax.set_title(title_str, fontsize=9, pad=title_pad)
    # ax.set_title(title_str, fontsize=8, pad=title_pad)

    # Extent
    if set_global:
        ax.set_global()
    else:
        if all(val is not None for val in [latmin, latmax, lonmin, lonmax]):
            ax.set_extent([lonmin, lonmax, latmin, latmax], crs=ccrs.PlateCarree())

    plt.tight_layout()
    return fig, ax

start the routine

Reading the data

Here we points at a some sample orbit files of GOME-2 data on the downloaded and unzipped folder, builds the selected instrument and period and then lists the matching files on disk as a quick existence check before opening anything.

pathDATA='./data/SO2/' 

# ls -lh ./data/SO2/GOME_*.HDF5 # this only works for linux
G2files = glob.glob(pathDATA + 'GOME_*.HDF5')

There are several possibilities in python to read an HDF5 file.

We refer here to illustrations presented in the NO2 notebook (https://acsaf.org/notebooks.php) both using HARP or xarray. In the following, we will show the use of both without many explanation.

select instrument and period

Defines the parameters of this run: which satellite (metop), which orbit/time (period, format YYYYMMDDHHMMSS), which SO2 layer height to assume (SO2LH, in km — lower for anthropogenic sources, higher for volcanic plumes), and which region (region). The values are hard-coded here, with several alternative example cases left commented out for quick switching.

Norilsk and Kivu/DR Congo are the two pre-configured cases for anthropogenic emission and volcanic emissions, respectively.

# comment and uncomment both region/SO2LH together, and then metop and period together

# SO2LH = 1  # Anthropogenic
# # SO2LH = 5  # Anthropogenic
# region='Norilsk'

# metop='B'
# period='20250324062345' #-> ok for both
# period='20250404055437' #--> this one gives empty SO2_vect for Kivu METOPB
# period='20250406065347' #-> ok or both
# period='20250407063246'  #-> ok or both

metop='C'
period='20250324071642'
# period='20250404064714'
# period='20250406074625'
# period='20250407072519'

SO2LH = 10 # Volcanic
region='Kivu'

get the data from the path and separate the file name into its different parts

# files = glob.glob(pathDATA+'GOME*SO2*'+period+'*.HDF5')
files = glob.glob(pathDATA+'GOME*SO2*'+period+'*'+metop+'*.HDF5')
print(files)

# extract the orbit
parts=files[0].split('_')
orbit = parts[-3]

['./data/SO2/GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20250324071642_051_METOPC_33092_DLR_04.HDF5']

read data (with xarray)

Reading approach: xarray: opens the HDF5 file with xr.open_dataset, targeting in turn the DETAILED_RESULTS, GEOLOCATION, CLOUD_PROPERTIES and DETAILED_RESULTS/SO2 groups.

gome2 = xr.open_dataset(files[0])

gome2_res = xr.open_dataset(files[0], group='DETAILED_RESULTS')
gome2_cor = xr.open_dataset(files[0], group='GEOLOCATION')
gome2_cld = xr.open_dataset(files[0], group='CLOUD_PROPERTIES')
gome2_so2 = xr.open_dataset(files[0], group='DETAILED_RESULTS/SO2')

and get all the needed variables

We extract the SO2 column at the 4 standard assumed layer heights (15 km, 6 km, 2 km, 1 km) and convert them from molecules/cm² to DU. These will later be used to interpolate to the height of interest. We also retrieve the 4 pixel-corner coordinates (A/B/C/D → SW/SE/NE/NW) and the SZA. loncen is wrapped into [-180, 180] (the source files sometimes use [0, 360]).

gome2_so2_15km = gome2_so2['VCDCorrected'][:,0] / molecPERcm2_to_DU
gome2_so2_6km  = gome2_so2['VCDCorrected'][:,1] / molecPERcm2_to_DU
gome2_so2_2km  = gome2_so2['VCDCorrected'][:,2] / molecPERcm2_to_DU
gome2_so2_1km  = gome2_so2['VCDCorrected'][:,3] / molecPERcm2_to_DU

gome2_so2_15km_err = gome2_so2['VCDCorrected_Error'][:,0] / molecPERcm2_to_DU
gome2_so2_6km_err  = gome2_so2['VCDCorrected_Error'][:,1] / molecPERcm2_to_DU
gome2_so2_2km_err  = gome2_so2['VCDCorrected_Error'][:,2] / molecPERcm2_to_DU
gome2_so2_1km_err  = gome2_so2['VCDCorrected_Error'][:,3] / molecPERcm2_to_DU

flagSO2 = gome2_so2['SO2_Flag']

latcen = gome2_cor['LatitudeCentre']
loncen = gome2_cor['LongitudeCentre']
loncen = (((loncen + 180) % 360) - 180)

latSW = gome2_cor['LatitudeA']
latSE = gome2_cor['LatitudeB']
latNE = gome2_cor['LatitudeD']
latNW = gome2_cor['LatitudeC']

lonSW = gome2_cor['LongitudeA']
lonSE = gome2_cor['LongitudeB']
lonNE = gome2_cor['LongitudeD']
lonNW = gome2_cor['LongitudeC']

sza = gome2_cor['SolarZenithAngleCentre']

Filter the data

We want to estimate the best SZA limit, which depend on the day of the year. For this, we compute the decimal day-of-year (integer DOY plus the fraction of the day already elapsed) from period, which is then used by the SZAlimitGOME2 function (defined earlier) to derive the SZA quality limit for this specific pixel time.

# Ground pixel date/time
dt = datetime.strptime(period, "%Y%m%d%H%M%S")
# integer DOY
doy = dt.timetuple().tm_yday
# seconds since midnight
seconds = dt.hour*3600 + dt.minute*60 + dt.second
# decimal DOY
decimal_doy = doy + seconds / 86400.0
# deriuve year froim the period:
YYYY=period[0:4]
MM=period[4:6]
DD=period[6:8]
(SZAnorth,SZAsouth) = SZAlimitGOME2(float(YYYY), decimal_doy)

flag due to high sza

We now build a flagVCD = 1 where the pixel's SZA is below the seasonal limit (north or south, chosen from the scene's mean latitude), 0 otherwise. We then cleans the data: negative SO2 columns (not physically possible, noise artifacts) are set to 0, and pixels with flagVCD < 0 are set to NaN.

flagVCD = np.empty(latcen.size)
if ma.mean(latcen) < 0:
    flagVCD[sza>SZAsouth] = 0
    flagVCD[sza<=SZAsouth] = 1
else:
    flagVCD[sza>SZAnorth] = 0
    flagVCD[sza<=SZAnorth] = 1

gome2_so2_15km_err[gome2_so2_15km<0] = 0.
gome2_so2_6km_err[gome2_so2_6km<0]   = 0.
gome2_so2_2km_err[gome2_so2_2km<0]   = 0.
gome2_so2_1km_err[gome2_so2_1km<0]   = 0.

gome2_so2_15km[gome2_so2_15km<0] = 0.
gome2_so2_6km[gome2_so2_6km<0]   = 0.
gome2_so2_2km[gome2_so2_2km<0]   = 0.
gome2_so2_1km[gome2_so2_1km<0]   = 0.

# mask = (flagVCD == 0) | (flagSO2 == 0) | (flagSolarEclipse > 0)
# mask = (flagVCD == 0) | (flagSO2 == 0)
mask = (flagVCD < 0)

gome2_so2_15km[mask] = np.nan
gome2_so2_6km[mask]  = np.nan
gome2_so2_2km[mask]  = np.nan
gome2_so2_1km[mask]  = np.nan

gome2_so2_15km_err[mask] = np.nan
gome2_so2_6km_err[mask]  = np.nan
gome2_so2_2km_err[mask]  = np.nan
gome2_so2_1km_err[mask]  = np.nan

select SO2 column (depending on the SO2 layer height mean values)

As discussed above, the idea is to linearly interpolates between the 2 standard layer heights bracketing SO2LH (the parameter chosen at the top of the notebook) to estimate the SO2 column at the actually assumed plume height. E.g. with SO2LH=8 (between 6 and 15 km), it interpolates between gome2_so2_6km and gome2_so2_15km.

# estimation of SO2 column using SO2 layer height mean values
gome2_so2_volc = gome2_so2_1km

if SO2LH >=1   and SO2LH < 2.5:
    gome2_so2_volc = gome2_so2_1km + (gome2_so2_2km - gome2_so2_1km) * (SO2LH - 1.)  / (2.5 - 1)
if SO2LH >=2.5 and SO2LH < 6.:
    gome2_so2_volc = gome2_so2_2km + (gome2_so2_6km - gome2_so2_2km) * (SO2LH - 2.5) / (6. - 2.5)
if SO2LH >=6.  and SO2LH < 15.:
    gome2_so2_volc = gome2_so2_6km + (gome2_so2_15km - gome2_so2_6km) * (SO2LH - 6) / (15. - 6.)
if SO2LH >= 15.:
    gome2_so2_volc = gome2_so2_15km

criteria to filter backscan

Typically, the "backscan" pixels are geometrically larger and of lower quality, so we compute each pixel's area (calcSO2massPXL/pixel_area) and excludes pixels whose area exceeds 9000 km². The mask is applied consistently to every geometry-related variable (corners, center, flag) and the corner arrays are then re-stacked (latcor, loncor) as needed for the polygon plot.

# criteria to filter backscan
(massTMP, area) = calcSO2massPXL(gome2_so2_volc, lonNW, lonNE, lonSE, lonSW, latNW, latNE, latSE, latSW)
mask = (area < 9000)

############
gome2_so2_volc = gome2_so2_volc[mask]
latSW = latSW[mask]
latSE = latSE[mask]
latNE = latNE[mask]
latNW = latNW[mask]
lonSW = lonSW[mask]
lonSE = lonSE[mask]
lonNE = lonNE[mask]
lonNW = lonNW[mask]

latcen = latcen[mask]
loncen = loncen[mask]

latcor = np.stack([latSW, latSE, latNE, latNW], axis=-1)
loncor = np.stack([lonSW, lonSE, lonNE, lonNW], axis=-1)

flagVCD = flagVCD[mask]
############

criteria for the lat/lon box used for mapping

We Define lat/lon boxes for the selected region: first a larger box (lonmin/max, latmin/max) used for the map display, and a tighter box (lonmini/maxi, latmini/maxi) that is used to select the pixels that go into the SO2 mass estimate. Norilsk (anthropogenic emission) and Kivu/DR Congo (volcanic) are the two pre-configured cases.

# define some limits for the lat/lon box used for mapping

# Anthropogenic emission in Norilsk
if region == "Norilsk":
  lonmin=70
  lonmax=100
  latmin=60
  latmax=80
  lonmini=78
  lonmaxi=90
  latmini=68
  latmaxi=74

# Volcanic emission in DR of Congo
if region == "Kivu":
  lonmin=10
  lonmax=40
  latmin=-10
  latmax=10
  lonmini=21
  lonmaxi=29.5
  latmini=-2.6
  latmaxi=2

We now restrict the data to the small lat/lon box (latmini/maxi, lonmini/maxi), which is the precise area where the SO2 mass should be estimated (e.g. just around the Norilsk/Kivu plume, not the whole display region).

These variables are then converted to plain numpy arrays (.values) for the calculations that follow.

# criteria for the lat/lon box used for mapping
mask = ((latcen < latmaxi) & (latcen > latmini) &
        (loncen < lonmaxi) & (loncen > lonmini))

##### Selection of SO2 column data for SO2 mass estimates
SO2_vect = gome2_so2_volc.values[mask.values]
lat_vect = latcen.values[mask.values]
lon_vect = loncen.values[mask.values]

latSW_vect = latSW.values[mask.values]
latSE_vect = latSE.values[mask.values]
latNE_vect = latNE.values[mask.values]
latNW_vect = latNW.values[mask.values]
lonSW_vect = lonSW.values[mask.values]
lonSE_vect = lonSE.values[mask.values]
lonNE_vect = lonNE.values[mask.values]
lonNW_vect = lonNW.values[mask.values]

flagVCD_vect = flagVCD[mask.values]

latcor_vect = np.stack([latSW_vect, latSE_vect, latNE_vect, latNW_vect], axis=-1)
loncor_vect = np.stack([lonSW_vect, lonSE_vect, lonNE_vect, lonNW_vect], axis=-1)

calculate SO2 max and do plot of the region

Now we follow a 5-step processing chain to calculate the maximum SO2 and prepare everything for the plot:

1. Quality filter — keep only pixels with flagVCD != 0 (valid SZA).
2. Noise estimate — standard deviation of values below 1 DU (presumed signal-free area).
3. Background removal — computes a "background" value (median of values under the noise/detection threshold) and zeroes out anything below it.
4. Thresholding — boolean mask of pixels above the threshold (3× noise, or the 0.45 DU detection limit).
5. Connected components (scipy.ndimage.label) — removes small isolated clusters under 20 connected pixels, filtering out residual noise that happens to cross the threshold in isolation, keeping only genuinely contiguous plumes.

The total SO2 mass is then computed (calcSO2massPXL, summed over positive pixels), the colour-scale max (vmax) is chosen dynamically based on the plume's intensity, a summary is printed (mass, area, max value and its location), the output PNG filename is built, and in the cell after the visualize_footprints function is run to produce the map.

# SO2 max
if SO2_vect.size > 0:
    INDmax = int(np.argmax(SO2_vect))
    SO2max = round(float(SO2_vect[INDmax]),2) 
    LONmax = round(float(lon_vect[INDmax]),2)
    LATmax = round(float(lat_vect[INDmax]),2)

    print('SO2 max: ', SO2max, 'DU  (', LONmax,'°E, ',LATmax,'°N )')

    # SO2 detection limit
    # if SO2LH < 15:
    #   SO2limit = round(0.45 * (15 / (15 - SO2LH)) * 15 / (15 - SO2LH), 3)
    # else:
    #   SO2limit = 0.45
    # print('SO2 limit: ', SO2limit, 'DU     (SO2 layer height of ', SO2LH,'km)')
    SO2limit = 0.45

    # 1. quality flags
    SO2_vect = np.where(flagVCD_vect != 0, SO2_vect, np.nan)

    # 2. noise estimate
    noise = np.nanstd(SO2_vect[SO2_vect < 1])
    threshold = np.max((SO2limit, 3 * noise))

    # 3. background removal
    # Filter values below threshold and non-zero
    filtered = SO2_vect[(SO2_vect < threshold) & (SO2_vect != 0)]
    background = np.median(filtered)
    # SO2_vect = SO2_vect - background
    SO2_vect = np.where(SO2_vect < background, 0, SO2_vect)

    # 4. threshold
    mask = SO2_vect > threshold

    # 5. connected components
    labeled, _ = label(mask)
    sizes = np.bincount(labeled.ravel())
    clean = sizes[labeled] > 20
    SO2_vect = np.where(clean, SO2_vect, 0)

    #  SO2 mass calculation
    (SO2mass, area) = calcSO2massPXL(SO2_vect, lonNW_vect, lonNE_vect, lonSE_vect, lonSW_vect, latNW_vect, latNE_vect, latSE_vect, latSW_vect)

    SO2mass[SO2mass != SO2mass] = 0
    area[area != area] = 0

    SO2mass = round(ma.sum(SO2mass[SO2_vect>0.]), 4) # in kt
    area = int(round(ma.sum(area[SO2_vect>0.]), 0)) # in km2

    print('SO2 mass: ', round(SO2mass,3), 'kt    area: ',area, 'm²')

    vmax = 10
    if SO2max > 20:
        vmax = 15
    if SO2max > 50:
        vmax = 40 
    if SO2max > 100:
        vmax = 50
    print('SO2 legend: from 0 to', vmax,'DU')

    print(YYYY, MM, DD, '12 00', f"{round(SO2mass*1000,1):9.1f}", f"{area:7d}",  f"{SO2max:6.1f}", f"{LATmax:6.2f}", f"{LONmax:7.2f}")

    mask = ~np.isnan(SO2_vect)
    lon_vect = lon_vect[mask]
    lat_vect = lat_vect[mask]
    SO2_vect = SO2_vect[mask]

    print(period[0:8]+'_GOME2-'+metop+'_SO2_orb-'+orbit+'_'+region)

SO2 max:  7.33 DU  ( 27.1 °E,  -2.2 °N )
SO2 mass:  3.451 kt    area:  197534 m²
SO2 legend: from 0 to 10 DU
2025 03 24 12 00    3451.0  197534    7.3  -2.20   27.10
20250324_GOME2-C_SO2_orb-33092_Kivu

Now we Call visualize_footprints with the cleaned, region-cropped SO2 field and its pixel corners, plus the mass/area/max statistics just computed, so they get included into the map title.

if SO2_vect.size > 0:
    fig, ax = visualize_footprints(
        gome2_so2_volc,
        loncor,
        latcor,
        period,
        color_scale='',
        cmap=cmap_so2,
        norm=norm_so2,
        unit='$\mathrm{SO_2}$ (DU)',
        long_name='SO2',
        vmin=0,
        vmax=vmax,
        latmin=latmin,
        latmax=latmax,
        lonmin=lonmin,
        lonmax=lonmax,
        metop=metop,
        orbit=orbit,
        SO2mass=SO2mass,
        SO2area=area,
        SO2max=SO2max,
        SO2max_lat=LATmax,
        SO2max_lon=LONmax,
      )
else:
    print('SO2_vect is empty...for %s, %s'% (region, period))
    fig=plt.figure()

We can now save the figure as a PNG file:

out_png = period[0:8]+'_GOME2-'+metop+'_SO2_orb-'+orbit+'_'+region+'_SO2LH-'+str(SO2LH)+'.png'
fig.savefig(out_png, dpi=150, bbox_inches='tight')
print(f"Saved: {out_png}")

Saved: 20250324_GOME2-C_SO2_orb-33092_Kivu_SO2LH-10.png

end of the notebook