Skip to content

Sentinel-5P and Python

  • blog
Sentinel_5P

Sentinel_5P

The data for Sentinel 5P was made available on 11th July 2018.

This is a ground breaking sensor and is offering unprecedented measurements of our atmosphere. If you want to know more about this sensor, I have previously written here about why Sentinel-5P is important.

Ever since its launch on the 13th October 2017 the data that has been released has looked incredible:

and,

and,

Stunning images and, behind it, science that informs and excites. All with the free and open policy of Copernicus; this is a golden age of satellites.

Data

The data is now available to download from the scihub here. It is currently in the ‘Pre Ops’ phase so we can expect improvements shortly. At the time of writing there are no quicklooks on the data available, but if you are already familiar with downloading Sentinel 1,2&3 data, navigation is identical.

Sentinel-5P data is supplied in .nc format. This means netcdf 4.

I downloaded an NO2 level 2 file from the sci hub that covered the UK. The file is an image from 16th July 2018.

Viewing the data

I tried to load the .nc in QGIS 3.2 but it would not work for me. The only ‘easy’ viewer I found was Panoply. More information below:

https://www.giss.nasa.gov/tools/panoply/

I was able to display the data in panoply as shown below:

I should add that the ever impressive SNAP toolbox was also able to open and view this data. However, I wondered if I could access it in Python. This is new data from a new sensor and I couldn’t find much about doing this so I wondered if I could inspect and plot the data. It turns out it is simple to do and I’ll show you how next.

Using Sentinel-5P data in Python

The netcdf files can be read in Python using the netcdf4 library. Using anaconda you can install this using the following command:

conda install -c anaconda netcdf4

Test to see if you can import the library; if you get no errors you are good to go.

import netCDF4 as nc4

Next, let’s read the file in:

from netCDF4 import Dataset
import numpy as np

my_example_nc_file = 'S5P_NRTI_L2__NO2____20180716T122509_20180716T123009_03918_01_010002_20180716T131625.nc'
fh = Dataset(my_example_nc_file, mode='r')
print (fh)

Print (fh) to see if the file has been loaded. This should return a lot of information about the file, such as when it was acquired, what sensor and its version, among many other things.

Next, let’s drill into the data.

print (fh.groups)

print (fh.groups['PRODUCT'])

By printing the fh.groups we can find the groups in the data, which are called ‘PRODUCT’ and ‘METADATA’. By printing fh.groups[‘PRODUCT’] we find the product has data about the dimension and variables.

print (fh.groups['PRODUCT'].variables.keys())

By printing the .variables.keys() we can see that this .nc file has the following variables:

‘scanline’, ‘ground_pixel’, ‘time’, ‘corner’, ‘polynomial_exponents’, ‘intensity_offset_polynomial_exponents’, ‘layer’, ‘vertices’, ‘latitude’, ‘longitude’, ‘delta_time’, ‘time_utc’, ‘qa_value’, ‘nitrogendioxide_tropospheric_column’, ‘nitrogendioxide_tropospheric_column_precision’, ‘nitrogendioxide_tropospheric_column_precision_kernel’, ‘averaging_kernel’, ‘air_mass_factor_troposphere’, ‘air_mass_factor_total’, ‘tm5_tropopause_layer_index’, ‘tm5_constant_a’, ‘tm5_constant_b’

If we select one of these variables, say the ‘nitrogendioxide_tropospheric_column_precision’ using this command…

print (fh.groups['PRODUCT'].variables['nitrogendioxide_tropospheric_column_precision'])

…we will finally get information about the data itself. You will see that the data has a value called current shape which we can convert into a numpy array.

I am indebted to the following post for the remaining code: http://joehamman.com/2013/10/12/plotting-netCDF-data-with-Python/ I have made only a few slight changes, otherwise I replicate it below.

The main change is that the above blog assumes a 2D array and the Sentinel-5P data output is a 3D array so we need to use the [0,:,:] code to slice it to a 2D array. If in doubt print the shape.

lons = fh.groups['PRODUCT'].variables['longitude'][:][0,:,:]
lats = fh.groups['PRODUCT'].variables['latitude'][:][0,:,:]
no2 = fh.groups['PRODUCT'].variables['nitrogendioxide_tropospheric_column_precision'][0,:,:]
print (lons.shape)
print (lats.shape)
print (no2.shape)


no2_units = fh.groups['PRODUCT'].variables['nitrogendioxide_tropospheric_column_precision'].units

The final line gets the units of NO2 so I can add to the plot later on. I am using Matplotlib to make the plots in a Jupyter Notebook, so I need to use the inline command at the start of the script; the final import helps plot the basemap using the mpl_toolkits.basemap.

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from mpl_toolkits.basemap import Basemap

The final part of the code is taken entirely from the blog post listed above.

lon_0 = lons.mean()
lat_0 = lats.mean()

m = Basemap(width=5000000,height=3500000,
            resolution='l',projection='stere',\
            lat_ts=40,lat_0=lat_0,lon_0=lon_0)

xi, yi = m(lons, lats)

# Plot Data
cs = m.pcolor(xi,yi,np.squeeze(no2),norm=LogNorm(), cmap='jet')

# Add Grid Lines
m.drawparallels(np.arange(-80., 81., 10.), labels=[1,0,0,0], fontsize=10)
m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10)

# Add Coastlines, States, and Country Boundaries
m.drawcoastlines()
m.drawstates()
m.drawcountries()

# Add Colorbar
cbar = m.colorbar(cs, location='bottom', pad="10%")
cbar.set_label(no2_units)

# Add Title
plt.title('NO2 in atmosphere')
plt.show()

This creates the following plot from the nitrogendioxide_tropospheric_column_precision variable.

Sentinel_5P

Europe is being imaged everyday now so it should be posssible to build a timeseries or to zoom into an area of interest, for example the south of the UK.

I am very excited about this sensor, I can see huge potential.

The code for this walkthrough is available on my Github page here.

Recap

  • Sentinel-5P data is now available to download
  • There are tools out there, namely panoply and SNAP that can display the data
  • Using the netcdf4 library Python can read, inspect and finally plot the data in matplotlib.

I don’t claim to be an atmospheric scientist, this is just a guide as to how to access the data. I hope it is helpful.

 

I am a freelancer able to help you with your projects. I offer consultancy, training and writing. I’d be delighted to hear from you. Please check out the books I have written on QGIS 3.4

https://www.packtpub.com/application-development/learn-qgis-fourth-edition

https://www.packtpub.com/application-development/qgis-quick-start-guide

 

I have grouped all my previous blogs (technical stuff / tutorials / opinions / ideas) at http://gis.acgeospatial.co.uk.

Feel free to connect or follow me; I am always keen to talk about Earth Observation.

I am @map_andrew on twitter

Tags: