madrid_banner.png

Comparison of OMI and TROPOMI spatial resolution of NO2 observations at city scale¶

This use case compares the spatial resolution of the Ozone Monitoring Instrument (OMI) and the TROPOspheric Monitoring Instrument nitrogen dioxide (NO2) observations. OMI was launched onboard NASA's Aura satellite in 2004, and since that it has been providing observations on atmospheric composition now for nearly 18 years. TROPOMI was launched in 2017 onboard ESA Copernicus Sentinel 5p satellite. TROPOMI, as the successor of the OMI instrument, is very similar in design but has e.g. much higher spatial resolution. The spatial resolution of OMI NO2 data is 13km x 24 km at nadir, whereas for TROPOMI the resolution (since August 2019) is 5.5 km × 3.5 km. Both instruments are on an afternoon orbit and their local overpassing times are very close to each other. It should be noted that OMI has been suffering from a so called row anomaly that affects certain pixels of the orbit, and hence data can not be obtained for those pixels affected.

This notebook illustrates how to compare the spatial resolution of these two instruments at a city scale. The new features introduced in this notebook are:

  • Reading NO2 data from two different satellite instruments
  • Differences in quality filtering of NO2 data for TROPOMI and OMI
  • Derivation of a new parameter from the original data for pixel area
  • Plotting pixel-based data on a map

Table of contents¶

  1. Importing Python modules
  2. Obtaining data files for OMI and TROPOMI
  3. HARP import for OMI and TROPOMI NO2
  4. Creating GeoDataFrames for plotting
  5. Plotting data on a map
  6. References

1. Importing Python modules ¶

In this notebook new modules "folium", "branca", and "geopandas" are used. Folium is used as a basemap for plotting, and geopandas are used to modify the output data from HARP import.

import harp
import numpy as np
import matplotlib.pyplot as plt
import branca
import folium
import pyproj
import geopandas as gpd
from shapely.geometry import Polygon
import sentinelsat
import avl

2. Obtaining data files for OMI and TROPOMI ¶

The TROPOMI data file can be downloaded automatically using the sentinelsat library as in use case 2. The TROPOMI file used in this exercise is:

S5P_OFFL_L2__NO2____20210713T113312_20210713T131441_19424_02_020200_20210715T052115.nc

The OMI file is downloaded from a cache of example data that is used for use cases like this. The original OMI data can be found on NASA GES DISC. The specific data set can be found by searching "OMNO2". To download the data one needs to have registered to NASA Earthdata. Note that currently there are two data sets for Level 2 OMI NO2 data available, make sure that the data file is named OMNO2. The OMI file that is used in this example is:

OMI-Aura_L2-OMNO2_2021m0713t1219-o90393_v003-2021m0824t153120.he5

filename_tropomi = "S5P_OFFL_L2__NO2____20210713T113312_20210713T131441_19424_02_020200_20210715T052115.nc"
filename_omi = "OMI-Aura_L2-OMNO2_2021m0713t1219-o90393_v003-2021m0824t153120.he5"
api = sentinelsat.SentinelAPI('s5pguest', 's5pguest', 'https://s5phub.copernicus.eu/dhus')
result = api.download_all(api.query(filename=filename_tropomi))
result = avl.download(filename_omi)

3. HARP import for OMI and TROPOMI NO2 ¶

The advantage of HARP is that for the most part same operations can be used to import both TROPOMI and OMI NO2 products because the naming conventions are similar even though the original file structures are different. Hence, for both data sets operations in HARP import are very similar, the main difference is quality filtering which is different for OMI and TROPOMI. For both data set the parameter we are interested in is called "tropospheric_NO2_column_number_density". For TROPOMI the data quality for each pixel is indicated with a qa value varying from 0 to 100 (or 0 to 1), whereas for OMI the data quality is indicated in bitwise manner.

To follow the recommendations of the algorithm teams, the following quality filters should be included in operations:

  • TROPOMI NO2: "tropospheric_NO2_column_number_density_validity>75"
  • OMI NO2: "validity !& 1"

For this example we are interested only data in the vicinity of Madrid, therefore the imported data is restricted to area lat 39.5N...41N, lon -5W...-2W. In the operations we also derive a new parameter, pixel area as km^2, that can be derived from latitude and longitude bounds:

  • "derive(area {time} [km2])"; where area is the name of the parameter, time is the dimension of the imported data and km2 are the units.
operations_trop = ";".join([
    "tropospheric_NO2_column_number_density_validity>75",
    "latitude>39.5",
    "latitude<41",
    "longitude<-2",
    "longitude>-5",
    "keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density)",
    "derive(tropospheric_NO2_column_number_density [Pmolec/cm2])",
    "derive(area {time} [km2])",
])
tropomi_NO2 = harp.import_product(filename_tropomi, operations=operations_trop)
print(tropomi_NO2)
source product = 'S5P_OFFL_L2__NO2____20210713T113312_20210713T131441_19424_02_020200_20210715T052115.nc'
history = "2023-02-24T11:26:30Z [harp-1.17] harp.import_product('S5P_OFFL_L2__NO2____20210713T113312_20210713T131441_19424_02_020200_20210715T052115.nc',operations='tropospheric_NO2_column_number_density_validity>75;latitude>39.5;latitude<41;longitude<-2;longitude>-5;keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density);derive(tropospheric_NO2_column_number_density [Pmolec/cm2]);derive(area {time} [km2])')"

float latitude_bounds {time=778, 4} [degree_north]
float longitude_bounds {time=778, 4} [degree_east]
double tropospheric_NO2_column_number_density {time=778} [Pmolec/cm2]
double area {time=778} [km2]

operations_omi = ";".join([
    "validity !& 1",
    "latitude>39.5",
    "latitude<41",
    "longitude<-2",
    "longitude>-5",
    "keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density)",
    "derive(tropospheric_NO2_column_number_density [Pmolec/cm2])",
    "derive(area {time} [km2])",
])
omi_NO2 = harp.import_product(filename_omi, operations=operations_omi)
print(omi_NO2)
source product = 'OMI-Aura_L2-OMNO2_2021m0713t1219-o90393_v003-2021m0824t153120.he5'
history = "2023-02-24T11:26:30Z [harp-1.17] harp.import_product('OMI-Aura_L2-OMNO2_2021m0713t1219-o90393_v003-2021m0824t153120.he5',operations='validity !& 1;latitude>39.5;latitude<41;longitude<-2;longitude>-5;keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density);derive(tropospheric_NO2_column_number_density [Pmolec/cm2]);derive(area {time} [km2])')"

double longitude_bounds {time=100, 4} [degree_east]
double latitude_bounds {time=100, 4} [degree_north]
double tropospheric_NO2_column_number_density {time=100} [Pmolec/cm2]
double area {time=100} [km2]

4. Creating GeoDataFrames for plotting ¶

This use case uses the folium module as a basemap for plotting. Folium makes it easy to visualize data that’s been manipulated in Python on an interactive Leaflet map. For plotting purposes the HARP-imported OMI and TROPOMI data are put into GeoDataFrames.

TROPOMI GeoDataFrame:

First step is to create an empty GeoDataFrame for TROPMI data. For GeoDataFrame you need to define the The Coordinate Reference System (crs) that defines how the coordinates relate to places on the Earth. For that first an empty column named "geometry" is created to the data frame, and then the crs is set.

crs = 'epsg:4326'
tropomi = gpd.GeoDataFrame()
tropomi['geometry'] = None
tropomi.crs = crs

The loop over imported TROPOMI data creates a polygon from each individual pixel using latitude and longitude bounds. Each pixel is represented as a polygon having a corresponding value for tropospheric NO2 and pixel area, as well as pixel id.

# index for the for loop to go through all TROPOMI pixels
idx = tropomi_NO2.latitude_bounds.data.shape[0]

# loop over pixels
for x in range(idx):
    lat_c = tropomi_NO2.latitude_bounds.data[x]
    lon_c = tropomi_NO2.longitude_bounds.data[x]
    poly = Polygon(zip(lon_c, lat_c))
    tropomi.loc[x, 'pixel_id'] = x
    tropomi.loc[x, 'geometry'] = poly
    tropomi.loc[x, 'NO2'] = tropomi_NO2.tropospheric_NO2_column_number_density.data[x]
    tropomi.loc[x, 'pixelarea'] = tropomi_NO2.area.data[x]

# print the tropomi GeoDataFrame
print(tropomi)    
                                              geometry  pixel_id       NO2  \
0    POLYGON ((-2.13517 39.46452, -2.04550 39.50186...       0.0  1.185136   
1    POLYGON ((-2.04550 39.50186, -1.95687 39.53862...       1.0  0.617515   
2    POLYGON ((-2.24938 39.47373, -2.15859 39.51170...       2.0  1.217406   
3    POLYGON ((-2.15859 39.51170, -2.06888 39.54907...       3.0  1.040929   
4    POLYGON ((-2.06888 39.54907, -1.98020 39.58585...       4.0  0.983635   
..                                                 ...       ...       ...   
773  POLYGON ((-4.86402 40.87101, -4.75064 40.92185...     773.0  1.258333   
774  POLYGON ((-4.75064 40.92185, -4.63892 40.97171...     774.0  1.009654   
775  POLYGON ((-5.00508 40.86570, -4.88993 40.91758...     775.0  1.393554   
776  POLYGON ((-4.88993 40.91758, -4.77650 40.96844...     776.0  1.096100   
777  POLYGON ((-5.03110 40.91214, -4.91590 40.96405...     777.0  0.750597   

     pixelarea  
0    48.723557  
1    48.099241  
2    49.369719  
3    48.727886  
4    48.102719  
..         ...  
773  61.653107  
774  60.652073  
775  62.607019  
776  61.577735  
777  62.632193  

[778 rows x 4 columns]

The GeoDataFrame for OMI is created similar way:

omi = gpd.GeoDataFrame()
omi['geometry'] = None
omi.crs = crs

idx2 = omi_NO2.latitude_bounds.data.shape[0]

# loop over pixels
for x in range(idx2):
    lat_c = omi_NO2.latitude_bounds.data[x]
    lon_c = omi_NO2.longitude_bounds.data[x]
    poly = Polygon(zip(lon_c, lat_c))
    omi.loc[x, 'pixel_id'] = x
    omi.loc[x, 'geometry'] = poly
    omi.loc[x, 'NO2'] = omi_NO2.tropospheric_NO2_column_number_density.data[x]
    omi.loc[x, 'pixelarea'] = omi_NO2.area.data[x]

# print the tropomi GeoDataFrame
print(omi)    
                                             geometry  pixel_id       NO2  \
0   POLYGON ((-2.30579 39.45745, -1.99512 39.52248...       0.0  0.012496   
1   POLYGON ((-3.00057 39.43277, -2.66938 39.50672...       1.0  1.179602   
2   POLYGON ((-2.66938 39.50672, -2.34890 39.57602...       2.0  1.327874   
3   POLYGON ((-2.34890 39.57602, -2.03774 39.64116...       3.0 -0.913196   
4   POLYGON ((-3.38886 39.47174, -3.04487 39.55107...       4.0  0.537193   
..                                                ...       ...       ...   
95  POLYGON ((-4.68468 40.70724, -4.30332 40.80120...      95.0  0.616325   
96  POLYGON ((-4.30332 40.80120, -3.93834 40.88830...      96.0 -0.096588   
97  POLYGON ((-3.93834 40.88830, -3.58763 40.96935...      97.0 -1.720481   
98  POLYGON ((-5.13347 40.72291, -4.73266 40.82484...      98.0  0.528854   
99  POLYGON ((-4.73266 40.82484, -4.35069 40.91897...      99.0 -0.373166   

     pixelarea  
0   378.216118  
1   405.216951  
2   390.802544  
3   378.228071  
4   421.727840  
..         ...  
95  462.282766  
96  440.643172  
97  421.762061  
98  487.117973  
99  462.255545  

[100 rows x 4 columns]

Next we can compare the pixel areas of TROPOMI and OMI by plotting the "pixelarea" variable

fig = plt.figure(figsize=(20, 10))
plt.subplot(121)
plt.hist(tropomi.pixelarea, bins=range(25, 550, 15))
plt.xticks(np.arange(0,550, step=50))
plt.xlabel('Pixel area [km^2]')
plt.title('TROPOMI')

plt.subplot(122)
plt.hist(omi.pixelarea, bins=range(25, 550, 15))
plt.xticks(np.arange(0,550, step=50))
plt.xlabel('Pixel area [km^2]')
plt.title('OMI')

plt.show()

5. Plotting data on a map ¶

In this use case we are illustrating the satellite data resolution at a city scale, and therefore we are using folium module for plotting. For the folium map the center location is given to Madrid and zoom_start level is set to 10. In the first plot both omi and tropomi pixels are plotted as polygons. The styles (namely color of pixel edges) for each data (omi=plum, tropomi=purple) are defined with style function. The map can be zoomed.

madrid_map = folium.Map(location=[40.416775, -3.703790], zoom_start=8, tiles='Stamen Terrain')

style1 = {'fillColor': 'none', 'color': '#dd3497'}
style2 = {'fillColor': 'none', 'color': '#4a1486'}

folium.GeoJson(tropomi, style_function=lambda x:style1).add_to(madrid_map)
folium.GeoJson(omi, style_function=lambda x:style2).add_to(madrid_map)
display(madrid_map)
Make this Notebook Trusted to load map: File -> Trust Notebook