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:
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
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)
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:
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:
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]
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()
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)