Title

This notebook is an AC SAF notebook prepared by BIRA-IASB (Gaia Pinardi), on illustrating the use of:

• reading and plotting of the AC SAF GOME-2 L2 data from the GDP 4.8/4.9 offline data using the HARP tools
• this notebook focus on NO2 and BrO VCD over one or several orbits.

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).

Table of Contents

This notebook contains several sections:

• Data description, with documentation, downloading informations and data file format
• The different imports and indications on how to set-up HARP
• Illustration on how to read the GOME-2 NO2 data with HARP
• Then how to filter the GOME-2 NO2 data with HARP and several ingestions options

• Then how to plot these data with Cartopy

• Illustration of a second example focusing on BrO and O3 for all the orbits of one day

• And finally some general considerations about HARP

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: GOME2A_201812.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:

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).

import sys,h5py,os,cf_units
import numpy as np
import datetime as dt
import warnings
#makes the plots interactive ...
#%matplotlib widget 

import glob

from netCDF4 import Dataset
# import xarray as xr
import datetime

from matplotlib import pyplot as plt #can also be written as "import matplotlib.pyplot as plt"

import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

How to setup HARP

**What is HARP?** HARP (https://atmospherictoolbox.org/harp/) is **one of the best tools to read GOME-2 ACSAF Level-2 HDF5 products**, because it already contains built-in ingestion definitions for many satellite products (including GOME-2, OMI, TROPOMI). **HARP converts the native product structure into a standardized format so that different instruments have the same variable names and dimensions.**

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

Command-line sanity check (harpcheck, HARP's CLI tool): run it on any HDF5 product to confirm the ingestion definition recognizes the file — replace file.HDF5 with an actual GOME-2 file path.

Reading the data

Here we points at a sample day of GOME-2 data on the downloaded and unzipped folder, builds the exact file name for that day/orbit and then lists the matching files on disk as a quick existence check before opening anything.

G2path='./data/GOME2A_201812/' 

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

# G2file='%s/GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201134102_051_METOPA_62878_DLR_04.HDF5' % G2path
G2file='%s/GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5' % G2path

There are several possibilities in python to read an HDF5 file. We show here a few classical ones (with h5py and netCDF4) that preserve the original file structure and then illustate the use of HARP.

First option: inspecting the data with h5py

import h5py

f = h5py.File(G2file, "r")
print(list(f.keys()))

['CLOUD_PROPERTIES', 'DETAILED_RESULTS', 'GEOLOCATION', 'META_DATA', 'TOTAL_COLUMNS']

Lists the sub-groups (variable names) inside DETAILED_RESULTS and GEOLOCATION, the two groups used throughout the rest of this notebook, as a quick way to see what's available without opening a separate HDF5 viewer.

print(list(f["DETAILED_RESULTS"].keys()))

['AAI', 'AMFToCloudTop', 'AMFToCloudTop_Error', 'AMFToGround', 'AMFToGround_Error', 'AMFTotal', 'AMFTotal_Error', 'BrO', 'Clouds', 'ESC', 'ESC_Error', 'FittingChiSquare', 'FittingGoodness', 'FittingNumberOfIterations', 'FittingRMS', 'H2O', 'HCHO', 'IntensityWeightedCloudFraction', 'NO2', 'O3', 'QualityFlags', 'SO2', 'SurfaceAlbedo', 'SurfaceConditionFlags', 'SurfaceHeight', 'SurfacePressure', 'VCD', 'VCD_Error']

print(list(f["GEOLOCATION"].keys()))

['EarthRadius', 'IndexInScan', 'InternalPixelId', 'LatitudeA', 'LatitudeB', 'LatitudeC', 'LatitudeCentre', 'LatitudeD', 'LineOfSightZenithAngleCentre', 'LineOfSightZenithAngleE', 'LineOfSightZenithAngleG', 'LineOfSightZenithAngleSatCentre', 'LineOfSightZenithAngleSatE', 'LineOfSightZenithAngleSatG', 'LongitudeA', 'LongitudeB', 'LongitudeC', 'LongitudeCentre', 'LongitudeD', 'OrbitalMode', 'PMDCenterLat', 'PMDCenterLon', 'PMDPixelMap', 'RelativeAzimuthCentre', 'RelativeAzimuthE', 'RelativeAzimuthG', 'RelativeAzimuthSatCentre', 'RelativeAzimuthSatE', 'RelativeAzimuthSatG', 'SAAFlag', 'SatelliteAltitude', 'SolarEclipseFlag', 'SolarZenithAngleCentre', 'SolarZenithAngleE', 'SolarZenithAngleG', 'SolarZenithAngleSatCentre', 'SolarZenithAngleSatE', 'SolarZenithAngleSatG', 'SubPixelInScan', 'SubSatellitePointLatitude', 'SubSatellitePointLongitude', 'SwathMode', 'Time', 'ViewMode']

Second option: inspecting the data with netCDF4

file = Dataset(G2file, diskless=True, persist=False)
groups = file.groups.keys()
print(groups)

dict_keys(['CLOUD_PROPERTIES', 'DETAILED_RESULTS', 'GEOLOCATION', 'META_DATA', 'TOTAL_COLUMNS'])

Reused from an EUMETSAT example notebook (credited in the comment): pretty-prints every group's variable names, 3 per line, so the whole file structure can be scanned at a glance via the netCDF4 interface.

# from https://gitlab.eumetsat.int/eumetlab/atmosphere/atmosphere/-/blob/master/20_data_exploration/211_Metop-A_GOME-2_NO2Tropo_L2_load_browse.ipynb
for group in groups:
    print('{}:'.format(group))
    # variables in a group
    variables = sorted(file[group].variables.keys())
    # pretty print of variables in 3 columns
    for v1, v2, v3 in zip(variables[::3], variables[1::3], variables[2::3]):
        print('   - {:<30}- {:<30}- {:<}'.format(v1, v2, v3))
    if (len(variables) % 3) == 1:
        print('   - {:<30}'.format(variables[-1]))
    elif (len(variables) % 3) == 2:
        print('   - {:<30}- {:<}'.format(*variables[-2:]))
    # white line between one line and another
    print()

CLOUD_PROPERTIES:
   - CloudFraction                 - CloudFraction_Error           - CloudMode
   - CloudOpticalThickness         - CloudOpticalThickness_Error   - CloudTopAlbedo
   - CloudTopAlbedo_Error          - CloudTopHeight                - CloudTopHeight_Error
   - CloudTopPressure              - CloudTopPressure_Error        - CloudType
   - PMDCloudFraction              

DETAILED_RESULTS:
   - AAI                           - AMFToCloudTop                 - AMFToCloudTop_Error
   - AMFToGround                   - AMFToGround_Error             - AMFTotal
   - AMFTotal_Error                - ESC                           - ESC_Error
   - FittingChiSquare              - FittingGoodness               - FittingNumberOfIterations
   - FittingRMS                    - IntensityWeightedCloudFraction- QualityFlags
   - SurfaceAlbedo                 - SurfaceConditionFlags         - SurfaceHeight
   - SurfacePressure               - VCD                           - VCD_Error

GEOLOCATION:
   - EarthRadius                   - IndexInScan                   - InternalPixelId
   - LatitudeA                     - LatitudeB                     - LatitudeC
   - LatitudeCentre                - LatitudeD                     - LineOfSightZenithAngleCentre
   - LineOfSightZenithAngleE       - LineOfSightZenithAngleG       - LineOfSightZenithAngleSatCentre
   - LineOfSightZenithAngleSatE    - LineOfSightZenithAngleSatG    - LongitudeA
   - LongitudeB                    - LongitudeC                    - LongitudeCentre
   - LongitudeD                    - OrbitalMode                   - PMDCenterLat
   - PMDCenterLon                  - PMDPixelMap                   - RelativeAzimuthCentre
   - RelativeAzimuthE              - RelativeAzimuthG              - RelativeAzimuthSatCentre
   - RelativeAzimuthSatE           - RelativeAzimuthSatG           - SAAFlag
   - SatelliteAltitude             - SolarEclipseFlag              - SolarZenithAngleCentre
   - SolarZenithAngleE             - SolarZenithAngleG             - SolarZenithAngleSatCentre
   - SolarZenithAngleSatE          - SolarZenithAngleSatG          - SubPixelInScan
   - SubSatellitePointLatitude     - SubSatellitePointLongitude    - SwathMode
   - ViewMode                      

META_DATA:
   - FWLowerBound                  - FWName                        - FWUpperBound
   - MainSpecies                   - SolarShifts                   - VCDQualityIndicator

TOTAL_COLUMNS:
   - BrO                           - BrO_Error                     - H2O
   - H2O_Error                     - HCHO                          - HCHO_Error
   - NO2                           - NO2Tropo                      - NO2Tropo_Error
   - NO2_Error                     - O3                            - O3_Error
   - SO2                           - SO2_Error

Suggested option: reading the data with HARP

import harp

product = harp.import_product(G2file)
print(product)

source product = 'GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5'

double datetime {time=15127} [seconds since 2000-01-01]
double longitude {time=15127} [degree_east]
double latitude {time=15127} [degree_north]
double longitude_bounds {time=15127, 4} [degree_east]
double latitude_bounds {time=15127, 4} [degree_north]
double sensor_solar_zenith_angle {time=15127} [degree]
double solar_zenith_angle {time=15127} [degree]
double viewing_zenith_angle {time=15127} [degree]
double relative_azimuth_angle {time=15127} [degree]
double BrO_column_number_density {time=15127} [molec/cm^2]
double BrO_column_number_density_uncertainty {time=15127} [molec/cm^2]
byte BrO_column_number_density_validity {time=15127} []
double H2O_column_density {time=15127} [kg/m^2]
double H2O_column_density_uncertainty {time=15127} [kg/m^2]
byte H2O_column_number_density_validity {time=15127} []
double HCHO_column_number_density {time=15127} [molec/cm^2]
double HCHO_column_number_density_uncertainty {time=15127} [molec/cm^2]
int HCHO_column_number_density_validity {time=15127} []
double NO2_column_number_density {time=15127} [molec/cm^2]
double NO2_column_number_density_uncertainty {time=15127} [molec/cm^2]
byte NO2_column_number_density_validity {time=15127} []
double tropospheric_NO2_column_number_density {time=15127} [molec/cm^2]
double tropospheric_NO2_column_number_density_uncertainty {time=15127} [molec/cm^2]
double O3_column_number_density {time=15127} [molec/cm2]
double O3_column_number_density_uncertainty {time=15127} [molec/cm2]
byte O3_column_number_density_validity {time=15127} []
double SO2_column_number_density {time=15127} [molec/cm2]
double SO2_column_number_density_uncertainty {time=15127} [molec/cm2]
int SO2_column_number_density_validity {time=15127} []
double cloud_fraction {time=15127} []
double cloud_fraction_uncertainty {time=15127} []
double cloud_top_pressure {time=15127} [hPa]
double cloud_top_pressure_uncertainty {time=15127} [hPa]
double cloud_top_height {time=15127} [km]
double cloud_top_height_uncertainty {time=15127} [km]
double cloud_top_albedo {time=15127} []
double cloud_top_albedo_uncertainty {time=15127} []
double cloud_optical_depth {time=15127} []
double cloud_optical_depth_uncertainty {time=15127} []
double absorbing_aerosol_index {time=15127} []
double surface_heigth {time=15127} [km]
double surface_pressure {time=15127} [hPa]
byte scan_subindex {time=15127}
byte scan_direction_type {time=15127}
long index {time=15127}

harp.import_product() reads the satellite file and converts it to the HARP data model, returning variables as numpy arrays. Notice that all pixels are flattened into a single time dimension — this is normal HARP behavior.

check the HARP version:

harp.version()

'1.23'

# ! harpconvert --version #check on unix server a related harp tool version!

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:

select only the variables I need

You can apply operations during import (very powerful feature). One of these is to keep only some of the variables:

operations = "keep(latitude,longitude,tropospheric_NO2_column_number_density, datetime,solar_zenith_angle, cloud_fraction, NO2_column_number_density_validity,NO2_column_number_density )"

product = harp.import_product(G2file, operations)

Displays the trimmed-down HARP product, containing only the variables listed in the keep(...) operation above.

Note that the operations string may contain several instructions, separated by ";":

operations = ";".join([
    "keep(latitude,longitude,datetime,"
    "NO2_column_number_density,"
    "tropospheric_NO2_column_number_density,"
    "solar_zenith_angle,"
    "cloud_fraction,"
    "NO2_column_number_density_validity)"
])

product = harp.import_product(G2file, operations)

print(operations)

keep(latitude,longitude,datetime,NO2_column_number_density,tropospheric_NO2_column_number_density,solar_zenith_angle,cloud_fraction,NO2_column_number_density_validity)

print(product)

source product = 'GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5'
history = "2026-09-28T09:00:10Z [harp-1.23] harp.import_product('./data/GOME2A_201812//GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5',operations='keep(latitude,longitude,datetime,NO2_column_number_density,tropospheric_NO2_column_number_density,solar_zenith_angle,cloud_fraction,NO2_column_number_density_validity)')"

double datetime {time=15127} [seconds since 2000-01-01]
double longitude {time=15127} [degree_east]
double latitude {time=15127} [degree_north]
double solar_zenith_angle {time=15127} [degree]
double NO2_column_number_density {time=15127} [molec/cm^2]
byte NO2_column_number_density_validity {time=15127} []
double tropospheric_NO2_column_number_density {time=15127} [molec/cm^2]
double cloud_fraction {time=15127} []

the HARP version is also visible with:

product.history

"2026-09-28T09:00:10Z [harp-1.23] harp.import_product('./data/GOME2A_201812//GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5',operations='keep(latitude,longitude,datetime,NO2_column_number_density,tropospheric_NO2_column_number_density,solar_zenith_angle,cloud_fraction,NO2_column_number_density_validity)')"

product

<Product variables=odict_keys(['datetime', 'longitude', 'latitude', 'solar_zenith_angle', 'NO2_column_number_density', 'NO2_column_number_density_validity', 'tropospheric_NO2_column_number_density', 'cloud_fraction'])>

check the loaded NO2 values:

product['NO2_column_number_density'].data, product['tropospheric_NO2_column_number_density'].data

(array([1.37736890e+15, 1.37126656e+15, 1.41298465e+15, ...,
        5.79471995e+15, 5.86677554e+15, 5.90898701e+15]),
 array([nan, nan, nan, ..., nan, nan, nan]))

'tropospheric_NO2_column_number_density' data are not given for all pixels, as it is not always valid. In that case, a nan is returned.

As discussed in the PUM:

"The tropospheric NO2 column is only available for measurements with latitude < 70º and SZA < 80º. The error in the retrieved tropospheric NO2 columns can become (very) large for measurements with cloud fractions > 20% (cloud radiance fraction > 50%), and these tropospheric NO2 measurements should normally not be used. Therefore, the NO2Tropo entry in the TOTAL_COLUMNS group is not available for these cloudy conditions. "

This is what we see above, with values put to nan for these cases. These case be filtered out using the NO2_Flag

( For expert users and to facilitate detailed data analyses, the tropospheric NO2 columns for cloudy conditions remain available via the VCDTropo entry in DETAILED_RESULTS/NO2. As discussed in the PUM, several variables are reported in several places in the files: The NO2Tropo entry in the TOTAL_COLUMNS group contains the same value as the VCDTropo entry in DETAILED_RESULTS/NO2. These two entries contain the tropospheric NO2 column density. )

NO2_Flag:

NO2 Flag is degined in the PUM as:

Filtering the data

**The data need to be filtered** The data need to be filtered for clouds and for quality of the tropopsheric NO2 variable. We use the `NO2_Flag`: The tropospheric NO2 column and pollution corrected total NO2 column should only be used for measurements with **NO2Tropo_Flag set to 1 or 2 (as discussed in the [PUM](https://acsaf.org/docs/pum/Product_User_Manual_NTO_OTO_Nov_2019.pdf))**. We can apply this quality filtering during HARP ingestion. **HARP performs the filtering before returning the data, which makes it very efficient.**

quality filtering during HARP ingestion

Some variables are optional:

filter on "validity":

from PUM: keep NO2_Flag=1 or 2. This is translated in HARP into the variable "NO2_column_number_density_validity", and the filtering can be done at the ingestion level:

operations_sel = ";".join([
    "NO2_column_number_density_validity <=2",
    "keep(latitude,longitude,datetime,"
    "NO2_column_number_density,"
    "tropospheric_NO2_column_number_density,"
    "solar_zenith_angle,"
    "cloud_fraction,"
    "NO2_column_number_density_validity)"
])

product_sel = harp.import_product(G2file, operations_sel)

other ingestion operations:

We create a third, stricter product: on top of the validity filter, we also drops pixels with SZA ≥ 75° or cloud fraction ≥ 0.2 at ingestion time:

operations2 = ";".join([
    "NO2_column_number_density_validity <= 2",     #see above, from PUM: keep NO2_Flag=1 or 2
    "solar_zenith_angle < 85",
    "cloud_fraction < 0.2",
    "keep(latitude,longitude,datetime,"
    "NO2_column_number_density,"
    "tropospheric_NO2_column_number_density,"
    "solar_zenith_angle, cloud_fraction, "
    "NO2_column_number_density_validity)"
])

product_sel2 = harp.import_product(G2file, operations2)

We can now compare the three products built so far side by side: the basic loading (the unfiltered product) and the distinct validity-flag values with different extra options. We show the pixel counts that each filter level keeps and confirm e.g. that product_sel/product_sel2 only contain flags ≤ 2.

print(product)

source product = 'GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5'
history = "2026-09-28T09:00:10Z [harp-1.23] harp.import_product('./data/GOME2A_201812//GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5',operations='keep(latitude,longitude,datetime,NO2_column_number_density,tropospheric_NO2_column_number_density,solar_zenith_angle,cloud_fraction,NO2_column_number_density_validity)')"

double datetime {time=15127} [seconds since 2000-01-01]
double longitude {time=15127} [degree_east]
double latitude {time=15127} [degree_north]
double solar_zenith_angle {time=15127} [degree]
double NO2_column_number_density {time=15127} [molec/cm^2]
byte NO2_column_number_density_validity {time=15127} []
double tropospheric_NO2_column_number_density {time=15127} [molec/cm^2]
double cloud_fraction {time=15127} []

print(product_sel)

source product = 'GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5'
history = "2026-09-28T09:00:11Z [harp-1.23] harp.import_product('./data/GOME2A_201812//GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5',operations='NO2_column_number_density_validity <=2;keep(latitude,longitude,datetime,NO2_column_number_density,tropospheric_NO2_column_number_density,solar_zenith_angle,cloud_fraction,NO2_column_number_density_validity)')"

double datetime {time=15051} [seconds since 2000-01-01]
double longitude {time=15051} [degree_east]
double latitude {time=15051} [degree_north]
double solar_zenith_angle {time=15051} [degree]
double NO2_column_number_density {time=15051} [molec/cm^2]
byte NO2_column_number_density_validity {time=15051} []
double tropospheric_NO2_column_number_density {time=15051} [molec/cm^2]
double cloud_fraction {time=15051} []

print(product_sel2)

source product = 'GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5'
history = "2026-09-28T09:00:12Z [harp-1.23] harp.import_product('./data/GOME2A_201812//GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5',operations='NO2_column_number_density_validity <= 2;solar_zenith_angle < 85;cloud_fraction < 0.2;keep(latitude,longitude,datetime,NO2_column_number_density,tropospheric_NO2_column_number_density,solar_zenith_angle, cloud_fraction, NO2_column_number_density_validity)')"

double datetime {time=6197} [seconds since 2000-01-01]
double longitude {time=6197} [degree_east]
double latitude {time=6197} [degree_north]
double solar_zenith_angle {time=6197} [degree]
double NO2_column_number_density {time=6197} [molec/cm^2]
byte NO2_column_number_density_validity {time=6197} []
double tropospheric_NO2_column_number_density {time=6197} [molec/cm^2]
double cloud_fraction {time=6197} []

len(product.tropospheric_NO2_column_number_density.data), len(product_sel.tropospheric_NO2_column_number_density.data), len(product_sel2.tropospheric_NO2_column_number_density.data)

(15127, 15051, 6197)

np.unique(product.NO2_column_number_density_validity.data)

array([0, 4], dtype=int8)

np.unique(product_sel.NO2_column_number_density_validity.data)

array([0], dtype=int8)

np.unique(product_sel2.NO2_column_number_density_validity.data)

array([0], dtype=int8)

Define the 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

operations

'keep(latitude,longitude,datetime,NO2_column_number_density,tropospheric_NO2_column_number_density,solar_zenith_angle,cloud_fraction,NO2_column_number_density_validity)'

operations2

'NO2_column_number_density_validity <= 2;solar_zenith_angle < 85;cloud_fraction < 0.2;keep(latitude,longitude,datetime,NO2_column_number_density,tropospheric_NO2_column_number_density,solar_zenith_angle, cloud_fraction, NO2_column_number_density_validity)'

# Define the ingestion options
ingest_options = ";".join(['detailed_results=NO2', 'corrected=true'])
# Import the product with options
productT = harp.import_product(G2file, operations=operations,options=ingest_options)
print(productT)

source product = 'GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5'
history = "2026-09-28T09:00:13Z [harp-1.23] harp.import_product('./data/GOME2A_201812//GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5',operations='keep(latitude,longitude,datetime,NO2_column_number_density,tropospheric_NO2_column_number_density,solar_zenith_angle,cloud_fraction,NO2_column_number_density_validity)',options='detailed_results=NO2;corrected=true')"

double datetime {time=15127} [seconds since 2000-01-01]
double longitude {time=15127} [degree_east]
double latitude {time=15127} [degree_north]
double solar_zenith_angle {time=15127} [degree]
double NO2_column_number_density {time=15127} [molec/cm^2]
byte NO2_column_number_density_validity {time=15127} []
double tropospheric_NO2_column_number_density {time=15127} [molec/cm^2]
double cloud_fraction {time=15127} []

As before, we can do a sanity-checks productT (built with the detailed_results=NO2;corrected=true ingestion options)

np.unique(productT.NO2_column_number_density_validity.data), len(productT.tropospheric_NO2_column_number_density.data), 

(array([0, 4], dtype=int8), 15127)

var2plot=productT.tropospheric_NO2_column_number_density.data[:]
print(var2plot)
np.nanmin(var2plot), np.nanmax(var2plot), np.nanmean(var2plot), var2plot.mean()

[nan nan nan ... nan nan nan]

(-3.198635876758323e+16, 3.735216234220749e+16, 282959913121168.6, nan)

Compare several of the loaded data (different options):

product['NO2_column_number_density'].data, product['tropospheric_NO2_column_number_density'].data

(array([1.37736890e+15, 1.37126656e+15, 1.41298465e+15, ...,
        5.79471995e+15, 5.86677554e+15, 5.90898701e+15]),
 array([nan, nan, nan, ..., nan, nan, nan]))

product_sel['NO2_column_number_density'].data, product_sel['tropospheric_NO2_column_number_density'].data

(array([1.37736890e+15, 1.37126656e+15, 1.41298465e+15, ...,
        5.79471995e+15, 5.86677554e+15, 5.90898701e+15]),
 array([nan, nan, nan, ..., nan, nan, nan]))

product_sel2['NO2_column_number_density'].data, product_sel2['tropospheric_NO2_column_number_density'].data

(array([1.77375089e+15, 1.81522364e+15, 1.81812824e+15, ...,
        5.53822235e+15, 5.54880086e+15, 5.20241282e+15]),
 array([nan, nan, nan, ..., nan, nan, nan]))

productT['NO2_column_number_density'].data, productT['tropospheric_NO2_column_number_density'].data

(array([nan, nan, nan, ..., nan, nan, nan]),
 array([nan, nan, nan, ..., nan, nan, nan]))

Plotting of the data

To plot the tropospheric NO2 data, we create a function that works on arrays. We provide as input some arguments such as the type of projection, the colorbar type, the colorbar min and max, ...

conversion_factor=1e15
projection=ccrs.PlateCarree()
vmin=0
vmax=1
point_size=20
cmap='viridis' #'jet' #

# provide some initial arrays, to use with the plotting function
lat = productT.latitude.data
lon = productT.longitude.data
var2plot = productT.tropospheric_NO2_column_number_density.data
title=' (from HARP)'
unit='molec/cm2'

# initialize to nan lat, lon and var2plot to avoid problems:

# provide some initial arrays, to use with the plotting function
lat = productT.latitude.data
lon = productT.longitude.data
var2plot = productT.tropospheric_NO2_column_number_density.data

create a fct for the plotting

def plot_map_GOME2(lon=lon, lat=lat,var2plot=var2plot,conversion_factor=1e15,unit='molec/cm2',cmap='viridis',
                   projection=ccrs.PlateCarree(),point_size=20,vmin=0,vmax=1,title='GOME-2 AC SAF L2'):
    """
    fig = plot_map_GOME2(lon=lon, lat=lat,var2plot=var2plot,conversion_factor=1e-15,unit='molec/cm2',cmap='viridis',
                    projection=ccrs.PlateCarree(),point_size=20,vmin=0,vmax=1,title='GOME-2 AC SAF L2')
        plot a global map of the L2 points (as scatter points)
    """
    fig=plt.figure(figsize=(20, 10))
    ax = plt.axes(projection=projection)
    ax.coastlines()
    ax.set_global()
    # plot pixel positions
    img = ax.scatter( lon, lat, c=var2plot/conversion_factor,
        cmap=plt.cm.get_cmap(cmap),
        marker='o',
        s=point_size,
        transform=ccrs.PlateCarree(),
        vmin=vmin,
        vmax=vmax
    )
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    plt.xlabel("Longitude", fontsize=16)
    plt.ylabel("Latitude", fontsize=16)
    cbar = fig.colorbar(img, ax=ax, orientation='horizontal', fraction=0.04, pad=0.1)
    cbar.set_label(str(conversion_factor) + ' ' + unit, fontsize=16)
    cbar.ax.tick_params(labelsize=14)
    ax.set_title(title, fontsize=20, pad=20.0)
    plt.show()
    
    return fig

plot from the HARP product:

We first pulls lat/lon/var2plot from productT and then we call plot_map_GOME2; the next two cells inspect the plotted array and its summary statistics.

lat = productT.latitude.data
lon = productT.longitude.data
var2plot = productT.tropospheric_NO2_column_number_density.data
title=' (from HARP)'

fig = plot_map_GOME2(lon=lon, lat=lat,var2plot=var2plot,conversion_factor=1e15,unit=unit,cmap='viridis',projection=ccrs.PlateCarree(),
                   point_size=20,vmin=0,vmax=1,title='GOME-2 AC SAF L2 NO2 '+title)

var2plot

array([nan, nan, nan, ..., nan, nan, nan])

np.nanmin(var2plot), np.nanmax(var2plot), np.nanmean(var2plot), var2plot.mean()

(-3.198635876758323e+16, 3.735216234220749e+16, 282959913121168.6, nan)

now do a plot of the filtered tropo NO2 data:

lat = product_sel2.latitude.data
lon = product_sel2.longitude.data
var2plot = product_sel2.tropospheric_NO2_column_number_density.data
title=' (HARP filtered)'

fig = plot_map_GOME2(lon=lon, lat=lat,var2plot=var2plot,conversion_factor=1e15,unit=unit,cmap='viridis',projection=ccrs.PlateCarree(),
                   point_size=20,vmin=0,vmax=1,title='GOME-2 AC SAF L2 NO2 '+title)
var2plot.min(), var2plot.max(), np.nanmean(var2plot), var2plot.mean()

(nan, nan, 416438653841538.94, nan)

We now have a map of tropospheric NO2 valid points for one GOME-2A orbit !

HARP loading of BrO, O3, NO2 (for all the daily orbits)

As mentionned above, the GOME-2 file contains information about several gases, and we can use similar approach to plot, e.g. BrO VCDs, but also keep O3 and NO2 values.

We can also load all the orbits of that specific day at once:

filesD = glob.glob('%s/*.HDF5'%G2path)

Defines the HARP ingestion for BrO (keeping O3 and NO2 alongside it, all read from DETAILED_RESULTS via the detailed_results=BrO;corrected=true options) and imports every file matching filesD at once, so product_bro ends up holding the whole day's orbits concatenated; the next cells sanity-check the BrO mean and plot it (unfiltered) on a map.

specific ingestion options: e.g., BrO

There are some ingestions options that are important in our BrO case:

for GOME-2 ACSAF NO₂ products, HARP already knows the product structure, but if you want to select the tropospheric NO₂ column, you must use the correct ingestion options (because the file contains several NO₂ quantities).

operations_bro = ";".join([
#     "solar_zenith_angle < 75",
#     "cloud_fraction < 0.2",
    "keep(latitude,longitude,datetime,"
    "solar_zenith_angle, cloud_fraction, "
    "BrO_column_number_density, BrO_column_number_density_validity,"
    "NO2_column_number_density, NO2_column_number_density_validity,"
    "O3_column_number_density, O3_column_number_density_validity)"
])

ingest_options_bro = ";".join(['detailed_results=BrO', 'corrected=true'])
# Import the product with options
product_bro = harp.import_product(filesD, operations=operations_bro,options=ingest_options_bro)
print(product_bro)

history = "2026-09-28T09:00:23Z [harp-1.23] harp.import_product('./data/GOME2A_201812/GOME_O3-NO2-NO2Tropo-BrO-SO2-H2O-HCHO_L2_20181201033250_051_METOPA_62872_DLR_04.HDF5',operations='keep(latitude,longitude,datetime,solar_zenith_angle, cloud_fraction, BrO_column_number_density, BrO_column_number_density_validity,NO2_column_number_density, NO2_column_number_density_validity,O3_column_number_density, O3_column_number_density_validity)',options='detailed_results=BrO;corrected=true')"

double datetime {time=226797} [seconds since 2000-01-01]
double longitude {time=226797} [degree_east]
double latitude {time=226797} [degree_north]
double solar_zenith_angle {time=226797} [degree]
double BrO_column_number_density {time=226797} [molec/cm^2]
byte BrO_column_number_density_validity {time=226797} []
double NO2_column_number_density {time=226797} [molec/cm^2]
byte NO2_column_number_density_validity {time=226797} []
double O3_column_number_density {time=226797} [molec/cm2]
byte O3_column_number_density_validity {time=226797} []
double cloud_fraction {time=226797} []

check the BrO values (and the O3 and NO2 ones):

np.nanmean(product_bro.BrO_column_number_density.data)

43263799303718.89

np.nanmean(product_bro.NO2_column_number_density.data), np.nanmean(product_bro.O3_column_number_density.data)

(2720505026029188.0, 7.703957739201935e+18)

do the plot, with the plotting function create above:

lat = product_bro.latitude.data
lon = product_bro.longitude.data
v_bro = product_bro.BrO_column_number_density.data
v_o3 = product_bro.O3_column_number_density.data
v_no2 = product_bro.NO2_column_number_density.data

fig = plot_map_GOME2(lon=lon, lat=lat,var2plot=v_bro,conversion_factor=1,unit=unit,cmap='viridis',projection=ccrs.PlateCarree(),
                   point_size=20,vmin=0,vmax=1e14,title='GOME-2 AC SAF L2 BrO VCD (no filters)')

do a map for the O3 columns. O3 is typically shown in Dobson Units, so convert the molec/cm2 to DU

1 Dobson Unit (DU) equals 2.687 × 10¹⁶ molecules of ozone per square centimeter (cm²)

conversion_factorDU=2.687e16  # 1DU = 2.6867e16 molecules/cm².

fig = plot_map_GOME2(lon=lon, lat=lat,var2plot=v_o3,conversion_factor=conversion_factorDU,unit=unit,cmap='viridis',projection=ccrs.PlateCarree(),
                   point_size=20,vmin=100,vmax=400,title='GOME-2 AC SAF L2 O3 VCD')

use HARP to do the unit conversion (O3 in DU):

HARP can be used to do the conversion (see here):

harp.convert_unit(from_unit, to_unit, values)

product_bro.O3_column_number_density.unit

'molec/cm2'

from_unit = product_bro.O3_column_number_density.unit
to_unit = 'DU'
v_o3_DU = harp.convert_unit(from_unit, to_unit, product_bro.O3_column_number_density.data)

fig = plot_map_GOME2(lon=lon, lat=lat,var2plot=v_o3_DU,conversion_factor=1,unit=to_unit,cmap='viridis',projection=ccrs.PlateCarree(),
                   point_size=20,vmin=100,vmax=400,title='GOME-2 AC SAF L2 O3 VCD')

you can change the Cartopy map projection type to better see the poles, and focus on the Southern Pole:

see the varius projections here: https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html

print("NumPy:", np.__version__)
import cartopy
print("Cartopy:", cartopy.__version__)

NumPy: 1.26.4
Cartopy: 0.18.0

In our servers we have a relative recent numpy version with an older cartopy version, and with projections different than ccrs.PlateCarree() we have a problem with the deprecated/old np.float syntax that Cartopy is using, which was removed from newer NumPy versions. This can be fixed locally with a local compatibility fix before calling the cartopy projection:

# projection2=ccrs.PlateCarree(),
# Compatibility fix for older Cartopy versions with recent NumPy
if not hasattr(np, 'float'):
    np.float = float  #temporary restore np.float, that we have in some Cartopy projections!

projection2=ccrs.Orthographic(central_latitude=-90.0) #map centered on the South Pole

fig = plot_map_GOME2(lon=lon, lat=lat,var2plot=v_o3_DU,conversion_factor=1,unit=to_unit,cmap='viridis',projection=projection2,
                   point_size=20,vmin=100,vmax=400,title='GOME-2 AC SAF L2 O3 VCD')

v_o3_DU.min(), v_o3_DU.max(), 

(76.74784851074219, 637.6615600585938)

# End of Notebook

raise RuntimeError("end of the notebook")

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-63-38e915c95553> in <module>
----> 1 raise RuntimeError("end of the notebook")

RuntimeError: end of the notebook

Considerations about HARP

import multiple orbits at once

harp will automatically concatenate them

files = [ "file1.HDF5", "file2.HDF5", "file3.HDF5" ]

product = harp.import_product(files)

Below, some reminders for discovering, for a given ingestion module, which options and variables are available. This us useful when you don't already know the exact operation/option names HARP expects.

# How to check available ingestion options:
# ! harpcheck file.HDF5
# ! harpcheck --list-variables file.HDF5

Define the 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

# Define the ingestion options
ingest_options = ";".join(['detailed_results=NO2', 'corrected=true'])
# Import the product with options
productT = harp.import_product(G2file, operations=operations,options=ingest_options)
print(productT)

gridding L2 to L3

One big advantage of HARP is gridding L2 pixels into L3 grids:

operations = ";".join([
    "bin_spatial(180,-90,1,360,-180,1)",
    "derive(latitude {latitude})",
    "derive(longitude {longitude})",
    "keep(latitude,longitude,tropospheric_NO2_column_number_density)"
])

product = harp.import_product(filename, operations)

Compared with reading HDF5 directly (xarray/h5py):

HARP automatically:

• renames variables consistently
• converts units
• flattens pixel dimensions
• handles QA flags
• merges files
• grids L2 → L3

This is exactly why it is used for multi-sensor intercomparisons (OMI / GOME-2 / TROPOMI).