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:
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 below. Other classical options (i.e., xarray) are shown in a companion notebook (also available here: https://acsaf.org/notebooks.php).
We thanks MariLiza Koukouli (AUTH) for discussions and input.
This notebook contains several sections:
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.
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
The GOME-2 GDP 4.8/4.9 L2 reprocessed datafiles are structures as one file per orbit, ie 15 files per day:
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.
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.
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")
HARP is a toolkit for reading, processing and inter-comparing satellite remote sensing data, model data, in-situ data, and ground based remote sensing data, and it contains interface descriptions to AC SAF data files. The toolkit is composed of a set of command line tools, a C library of analysis functions, and direct interfaces to Python, R, IDL, and MATLAB.
Install HARP (Python interface)
Typical install in a conda environment and check that it works:
conda install -c conda-forge harp
# conda install -c conda-forge harp
Imports the (just-installed) package
import harp
import coda
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.
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
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)
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
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.
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_20250324062345_051_METOPB_64934_DLR_04.HDF5']
You can find here https://stcorp.github.io/harp/doc/html/ingestions/index.html all the ingestions definitions. More particularly, for the GOME-2 AC SAF trace gases product: https://stcorp.github.io/harp/doc/html/ingestions/GOME2_L2_O3MOTO.html
you will find the common HARP conventions, and the mapping from the original product:
harp.version()
'1.23'
You can apply operations during import (very powerful feature).
keep only some of the variables (the SO2 column, but also the geolocations, SZA, flag...)derive operation For SO2, the various SO2 columns at different altitudes are available through ingestion options
The interface for harp.import_product (see here: https://stcorp.github.io/harp/doc/html/python.html) is as follows:
harp.import_product(filename, operations='', options='', reduce_operations='', post_operations='')
options expects a string with Ingestion module specific options that should be specified as a semi-colon separated string of key=value pairs
The SO2_Flag is translated in HARP into the variable "SO2_column_number_density_validity", and the filtering can be done at the ingestion level, keeping only the values not equal to 255 (we can further filter later on).
Other variables can be used to filter at the ingestion level, ie SZA ≤ 89°, forward scan only, ...
# SO2 at 15km, including the SO2 flag, latitudes and longitudes of pixel center and of the corners, sza, and a few more,
# as illustration on how to keep additional variables of interest (not used here below)
product15 = harp.import_product(files[0], options = 'detailed_results=SO2;corrected=true;so2_column=15km', \
operations = 'solar_zenith_angle <= 89; scan_direction_type==0; SO2_column_number_density_validity!& 255;\
derive(SO2_column_number_density [DU]);derive(SO2_slant_column_number_density [DU]);\
derive(SO2_column_number_density_uncertainty [DU]); \
keep(latitude,longitude, solar_zenith_angle,\
latitude_bounds, longitude_bounds, \
SO2_column_number_density_validity, \
SO2_column_number_density, SO2_column_number_density_uncertainty,\
SO2_slant_column_number_density, \
cloud_fraction,SO2_column_number_density_validity, index)')
print(product15)
source product = 'GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20250324062345_051_METOPB_64934_DLR_04.HDF5'
history = "2026-09-28T09:20:43Z [harp-1.23] harp.import_product('./data/SO2/GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20250324062345_051_METOPB_64934_DLR_04.HDF5',operations='solar_zenith_angle <= 89; scan_direction_type==0; SO2_column_number_density_validity!& 255; derive(SO2_column_number_density [DU]);derive(SO2_slant_column_number_density [DU]); derive(SO2_column_number_density_uncertainty [DU]); keep(latitude,longitude, solar_zenith_angle, latitude_bounds, longitude_bounds, SO2_column_number_density_validity, SO2_column_number_density, SO2_column_number_density_uncertainty, SO2_slant_column_number_density, cloud_fraction,SO2_column_number_density_validity, index)',options='detailed_results=SO2;corrected=true;so2_column=15km')"
double longitude {time=9325} [degree_east]
double latitude {time=9325} [degree_north]
double longitude_bounds {time=9325, 4} [degree_east]
double latitude_bounds {time=9325, 4} [degree_north]
double solar_zenith_angle {time=9325} [degree]
double SO2_column_number_density {time=9325} [DU]
double SO2_column_number_density_uncertainty {time=9325} [DU]
int SO2_column_number_density_validity {time=9325} []
double SO2_slant_column_number_density {time=9325} [DU]
double cloud_fraction {time=9325} []
long index {time=9325}
check the loaded SO2 values:
product15['SO2_column_number_density'].data, product15['SO2_column_number_density_uncertainty'].data, product15['SO2_column_number_density_validity'].data
(array([ 1.61267452, -0.30240387, -2.52558647, ..., -0.76954981,
-1.10726417, 0.21395172]),
array([1.22589741, 3.5754426 , 1.29847765, ..., 1.50114247, 1.39016521,
0.45483026]),
array([0, 0, 0, ..., 0, 0, 0], dtype=int16))
np.unique(product15['SO2_column_number_density_validity'].data)
array([ 0, 256], dtype=int16)
product15.index.data
array([ 613, 645, 677, ..., 14241, 14242, 14273], dtype=int32)
In case we would like to add other variables that are not in the HARP ingestion definition, these need to be ingested by CODA. We give below and illustration on how to do this, but we will not need these variables in the rest of the notebook
# if we would like extra variables that are not in the HARP ingestion definition, these need to be ingested by CODA
# Open the HDF5 file in read mode
with h5py.File(files[0], "r") as f:
main_species = f["/META_DATA/MainSpecies"][()] # Read the dataset
# Convert byte strings to normal strings (UTF-8)
main_species = np.array([s.decode("utf-8") for s in main_species])
# main_species = main_species.decode("utf-8") if isinstance(main_species, bytes) else main_species
# Find indices where the value is "SO2"
so2_indices = np.where(main_species == "SO2")[0]
# Handle the case where "SO2" is missing
if so2_indices.size > 0:
so2_index = so2_indices[0] # Get the first (and only) index
# open the file with coda, fetch(ie, read) the variables we are interested in, and then add it as a HARP variable
# (depending on the time coordinates, and only some indexes), in the previously HARP-loaded list of products:
with coda.open(files[0]) as pf:
# variables
# 0 is the first height, 15km
# so2_index is the SO2 species in the sequence
product15.chisq = harp.Variable(coda.fetch(pf, "DETAILED_RESULTS/FittingChiSquare")[product15.index.data,0], ["time"])
ESC = harp.Variable(coda.fetch(pf, "DETAILED_RESULTS/ESC")[product15.index.data,so2_index], ["time"])
ESC_Error = harp.Variable(coda.fetch(pf, "DETAILED_RESULTS/ESC_Error")[product15.index.data,so2_index], ["time"])
product15.SO2_slant_column_number_density_uncertainty = harp.Variable(ESC_Error.data * 0.01 * ESC.data,["time"])
print(product15)
source product = 'GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20250324062345_051_METOPB_64934_DLR_04.HDF5'
history = "2026-09-28T09:20:43Z [harp-1.23] harp.import_product('./data/SO2/GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20250324062345_051_METOPB_64934_DLR_04.HDF5',operations='solar_zenith_angle <= 89; scan_direction_type==0; SO2_column_number_density_validity!& 255; derive(SO2_column_number_density [DU]);derive(SO2_slant_column_number_density [DU]); derive(SO2_column_number_density_uncertainty [DU]); keep(latitude,longitude, solar_zenith_angle, latitude_bounds, longitude_bounds, SO2_column_number_density_validity, SO2_column_number_density, SO2_column_number_density_uncertainty, SO2_slant_column_number_density, cloud_fraction,SO2_column_number_density_validity, index)',options='detailed_results=SO2;corrected=true;so2_column=15km')"
double longitude {time=9325} [degree_east]
double latitude {time=9325} [degree_north]
double longitude_bounds {time=9325, 4} [degree_east]
double latitude_bounds {time=9325, 4} [degree_north]
double solar_zenith_angle {time=9325} [degree]
double SO2_column_number_density {time=9325} [DU]
double SO2_column_number_density_uncertainty {time=9325} [DU]
int SO2_column_number_density_validity {time=9325} []
double SO2_slant_column_number_density {time=9325} [DU]
double cloud_fraction {time=9325} []
long index {time=9325}
float chisq {time=9325}
float SO2_slant_column_number_density_uncertainty {time=9325}
now load in a similar way the SO2 VCD at 6km, 2.5km and 1km: we keep the same list of operations (to filter in the same way) and only change the soé_column of interest in the options:
operationsSO2= 'solar_zenith_angle <= 89; scan_direction_type==0; SO2_column_number_density_validity!& 255;\
derive(SO2_column_number_density [DU]);derive(SO2_slant_column_number_density [DU]);\
derive(SO2_column_number_density_uncertainty [DU]); \
keep(latitude,longitude, solar_zenith_angle,\
latitude_bounds, longitude_bounds, \
SO2_column_number_density_validity, \
SO2_column_number_density, SO2_column_number_density_uncertainty,\
SO2_slant_column_number_density, \
cloud_fraction,SO2_column_number_density_validity, index)'
# SO2 at 6km
product6 = harp.import_product(files[0], options = 'detailed_results=SO2;corrected=true;so2_column=6km',operations = operationsSO2)
# SO2 at 2.5km
product2p5 = harp.import_product(files[0], options = 'detailed_results=SO2;corrected=true;so2_column=2.5km',
operations =operationsSO2)
# SO2 at 1km
product1 = harp.import_product(files[0], options = 'detailed_results=SO2;corrected=true;so2_column=1km',
operations = operationsSO2)
print(product1)
source product = 'GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20250324062345_051_METOPB_64934_DLR_04.HDF5'
history = "2026-09-28T09:20:46Z [harp-1.23] harp.import_product('./data/SO2/GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20250324062345_051_METOPB_64934_DLR_04.HDF5',operations='solar_zenith_angle <= 89; scan_direction_type==0; SO2_column_number_density_validity!& 255; derive(SO2_column_number_density [DU]);derive(SO2_slant_column_number_density [DU]); derive(SO2_column_number_density_uncertainty [DU]); keep(latitude,longitude, solar_zenith_angle, latitude_bounds, longitude_bounds, SO2_column_number_density_validity, SO2_column_number_density, SO2_column_number_density_uncertainty, SO2_slant_column_number_density, cloud_fraction,SO2_column_number_density_validity, index)',options='detailed_results=SO2;corrected=true;so2_column=1km')"
double longitude {time=9325} [degree_east]
double latitude {time=9325} [degree_north]
double longitude_bounds {time=9325, 4} [degree_east]
double latitude_bounds {time=9325, 4} [degree_north]
double solar_zenith_angle {time=9325} [degree]
double SO2_column_number_density {time=9325} [DU]
double SO2_column_number_density_uncertainty {time=9325} [DU]
int SO2_column_number_density_validity {time=9325} []
double SO2_slant_column_number_density {time=9325} [DU]
double cloud_fraction {time=9325} []
long index {time=9325}
product15['longitude_bounds'].data
array([[119.95105743, 119.34223938, 117.55440521, 118.13777924],
[118.77355194, 118.11597443, 116.34420013, 116.97520447],
[117.60547638, 116.89801025, 115.14362335, 115.82321167],
...,
[ 32.41490555, 30.41965675, 30.30853653, 32.34132004],
[ 30.63964272, 28.69121933, 28.55161285, 30.53241539],
[ 32.37281036, 30.35896683, 30.24716377, 32.29969025]])
product15['latitude_bounds'].data
array([[ 72.52399445, 73.31822968, 73.20911407, 72.41389465],
[ 72.45469666, 73.24177551, 73.12216949, 72.3339386 ],
[ 72.378479 , 73.16152191, 73.03159332, 72.24719238],
...,
[-65.60639954, -65.57637787, -66.10956573, -66.17212677],
[-65.58220673, -65.52998352, -66.0357132 , -66.11894226],
[-65.95180511, -65.92266846, -66.45606232, -66.51761627]])
Extracts the SO2 column at the 4 standard assumed layer heights (15 km, 6 km, 2 km, 1 km) — GOME-2 reports SO2 under several plume-height hypotheses, since the retrieved column depends strongly on this a-priori assumption — and converts from molecules/cm² to DU. Also retrieves 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]).
Note that HARP group and re-arrange the order of the corners coordinates:
as described here: https://stcorp.github.io/harp/doc/html/ingestions/GOME2_L2_O3MOTO.html
# adapt here, from the harp reading, I need all these variables:
gome2_so2_15km = product15['SO2_column_number_density'].data
gome2_so2_6km = product6['SO2_column_number_density'].data
gome2_so2_2km = product2p5['SO2_column_number_density'].data
gome2_so2_1km = product1['SO2_column_number_density'].data
gome2_so2_15km_err = product15['SO2_column_number_density_uncertainty'].data
gome2_so2_6km_err = product6['SO2_column_number_density_uncertainty'].data
gome2_so2_2km_err = product2p5['SO2_column_number_density_uncertainty'].data
gome2_so2_1km_err = product1['SO2_column_number_density_uncertainty'].data
flagSO2 =product15['SO2_column_number_density_validity'].data
latcen = product15['latitude'].data #['LatitudeCentre']
loncen = product15['longitude'].data #['LongitudeCentre']
loncen = (((loncen + 180) % 360) - 180)
latSW = product15['latitude_bounds'].data[:,3] #['LatitudeA']
latSE = product15['latitude_bounds'].data[:,0] #['LatitudeB']
latNE = product15['latitude_bounds'].data[:,1] # ['LatitudeD']
latNW = product15['latitude_bounds'].data[:,2] #['LatitudeC']
lonSW = product15['longitude_bounds'].data[:,3] #['LongitudeA']
lonSE = product15['longitude_bounds'].data[:,0] #['LongitudeB']
lonNE = product15['longitude_bounds'].data[:,1] #['LongitudeD']
lonNW = product15['longitude_bounds'].data[:,2] #['LongitudeC']
sza = product15['solar_zenith_angle'].data #['SolarZenithAngleCentre']
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)
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
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
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.
In principle, this is already done within the HARP ingestion. We leave the calculation here for coherence with the SO2 xArray notebook.
# 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]
############
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]
SO2_vect = gome2_so2_volc[mask]
lat_vect = latcen[mask]
lon_vect = loncen[mask]
latSW_vect = latSW[mask]
latSE_vect = latSE[mask]
latNE_vect = latNE[mask]
latNW_vect = latNW[mask]
lonSW_vect = lonSW[mask]
lonSE_vect = lonSE[mask]
lonNE_vect = lonNE[mask]
lonNW_vect = lonNW[mask]
flagVCD_vect = flagVCD[mask]
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)
Now we follow a 5-step processing chain to calculate the maximum SO2 and prepare everything for the plot:
flagVCD != 0 (valid SZA).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: 5.09 DU ( 27.2 °E, -2.22 °N ) SO2 mass: 0.871 kt area: 58785 m² SO2 legend: from 0 to 10 DU 2025 03 24 12 00 871.1 58785 5.1 -2.22 27.20 20250324_GOME2-B_SO2_orb-64934_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-B_SO2_orb-64934_Kivu_SO2LH-10.png