Sentinel 5p raster creation

[luisgpaixao]

Hello.

I’ve been using ESA’s Sentinel data for some time, but now i’m interested in using the data of the new Sentinel 5p. I’ve been using different software (ArcGIS, QGis, Snap, Panoply, Gdal) to accomplish my goal, but no luck so far. What i want is to be able to extract a geo-referenced raster (Geotiff) from the sentinel 5p products so i can use it in any GIS software. The best i can do is export a raster with the variable i want, but not geo-referenced correctly.

Does anyone has similar goals, and knows how i can achieve mine?

Much appreciated

[marpet]

Because of its special characteristics, SNAP does not support Sentinel 5p.
On the Sentinel Online page are tools mentioned which can handle the data. Probably you want to try the BEAT Toolbox.

[luisgpaixao]

Hello again and thank for the reply.
Already tried the BEAT Tollbox (VISAN Package), but in all the documentation there is no mention of Sentinel 5p products being supported. In fact, this is what i get when i try the beatl2.ingest() function on a Level 2 S5p .nc file:

Beatl2Error: BEAT error (-2106): product file of this type not supported

Which is odd because ESA website mention the BEAT Toolbox “supports extraction, visualization and analysis of Sentinel-5 Precursor products”.
Have you worked with BEAT before? Can you give some tips or point me in the right direction on how i can achieve my goal?

Regards

[luisgpaixao]

Just to let you know i got a solution, finally. Not ideal but it works.

First i open the .nc file on snap, then i export it to NetCDF4-CF. I then open this last file in a tool called Weather and Climate Toolkit (WCT), and export the variable, in this case the Ozone column, to a GeoTiff.

Nevertheless, i sent an e-mail to ESA Support asking about the BEAT Toolbox. If i get more news i’ll update this topic.

Thanks and regards

[marpet]

No, I haven’t worked with BEAT before. So I can’t help, sorry.

But I think you can shorten your current workflow. You don’t need to export to NetCDF4-CF, you can export directly to GeoTIFF from SNAP.

[luisgpaixao]

I’ve tried that before, but the GeoTiff comes out in a wrong geolocation.

[marpet]

My export has the same geolocation as the S5p product. Maybe you want to have it in a different CRS? You can reproject it before doing the subset.
However, as long as you found a way to achieve your results.

[luisgpaixao]

The CRS is fine. WGS84. The geolocation is the problem, also the GeoTiff comes as an RGB image.
Maybe i’m doing something wrong in Snap. Can you tell me the steps you made to export to GeoTiff?

[marpet]

The problem is probably that the data is provided with lat/lon values per pixel.
That’s why you have an RGB with 3 channels. Two of them are lat and lon.

I have now also noticed (SNAP-945) that when doing a reprojection on the original data the northern part is cut-off.

The red shape represents the reprojected data and the white one is the original data.

There is still a SNAP-only way. But I’m not sure if it is better than your solution. It is a bit strange. Obviously, something is not working correctly. Might be because of the S5p data, or maybe there is a bug in SNAP. But on the other side, the behaviour is correct. SNAP tries to preserve the best geolocation which is available.

However, here is the SNAP-only solution:
The trick is to do the export two times.

First open the S5p product, then export only the ozone band to GeoTiff.

Answer the question regarding the flag dataset with no.

Now open the exported file in SNAP and do the same. Now you have to exclude the tie-point grids on the corresponding tab too.

After this, you should get the result you want.

[luisgpaixao]

Thanks Marco. The two export method, produces in fact a raster with one band. But if i open that raster in ArcGIS or QGIS it’s not in the correct location.
This is the SNAP view of that final raster (left) and the view in ArcGIS (right).

Now if i do the method i explained before, with the Weather and Climate Toolkit, this is what i get.

Not sure what is wrong, but SNAP doesn’t seem to work for me.

[jmendrok]

Hi Luis,

I had similar issues with S-3 OLCI data. Reprojecting to Geographical Lat/Lon (WGS84) in SNAP before exporting helped. See here.

wishes,
Jana

[luisgpaixao]

Thanks Jana. I now understand what Marco was telling about the north part being cut-off. Not just the north part but also a little bit on the east part as well. Every part of the original data that exceeds the footprint of the reprojected data is cut-off.

[marpet]

It seems that an issue is back that was resolved almost a year ago. I have reopened SNAP-839.
We are working on a more general solution to better handle pixel-based geocoding. Unfortunately, this needs some time. Not a simple task.

[svniemeijer]

I am one of the main developers of the ESA Atmospheric Toolbox (BEAT).

The toolbox definitely supports S5P, but you will have to use the HARP component to do this (https://github.com/stcorp/harp).
The beatl2 interfaces are obsolete and will be phased out in the near future.

The Python interfaces of HARP are also included in the VISAN application (there is even an ingestion dialog in the VISAN menu for it), but to integrate things in a workflow, it might be easier to just install HARP itself and use the command line tools of HARP.

Note that with the latest release of HARP we included a way to raster the S5P data.
The new operation in HARP for L3 gridding is called ‘bin_spatial’.

As an example, for the product S5P_NRTI_L2__O3_____20180712T064913_20180712T065413_03858_01_010000_20180712T072955.nc (available from s5phub) you can do the following:

harpconvert -a 'bin_spatial(241,24,0.1,361,74,0.1)' S5P_NRTI_L2__O3_____20180712T064913_20180712T065413_03858_01_010000_20180712T072955.nc l3.nc

This will grid all harp-supported quantities from the product to a 240x360 lat/lon grid with latitude edges [24, 24.1, 24.2, …, 47.9, 48] and longitude edges [74, 74.1, …, 109.9, 110].

The gridding algorithm uses an area weighted average and knows how to deal with satellite pixels that cross the dateline boundary or cover the poles.
You can adapt the resolution to whatever size you like. You can also combine the transformation with filtering and other operations. For example:

harpconvert -a 'O3_column_number_density_validity>50;keep(datetime_start,datetime_length,O3_column_number_density,latitude_bounds,longitude_bounds);derive(datetime_stop {time});exclude(datetime_length);bin_spatial(1201,24,0.02,1801,74,0.02);derive(latitude {latitude});derive(longitude {longitude});exclude(longitude_bounds,latitude_bounds)' S5P_NRTI_L2__O3_____20180712T064913_20180712T065413_03858_01_010000_20180712T072955.nc l3.nc

Furthermore, you can use harpmerge to combine data from multiple products into a single output. For example:

harpmerge -a 'latitude>9;latitude<21;longitude>19;longitude<31' -ap 'bin_spatial(101,10,0.1,101,20,0.1)' S5P_NRTI_L2__O3* l3.nc

More information about the HARP operations can be found at: HARP Operations

Note that we also plan to update the website of the Atmospheric Toolbox in the future and add examples of how to use the toolbox to access S5P data here as well. Unfortunately, due to unforeseen circumstances, it took a while for us to be able to start on that activity.

[svniemeijer]

If you want additional help on the atmospheric toolbox, I’d suggest to send an email to the beat email address that is mentioned on our toolbox page. BEAT is not part of SNAP or STEP and therefore not covered by this forum.

[luisgpaixao]

Thank you Sander. I’ll give it a go. Regards

[ABerni]

Hello svniemeijer,

I am trying to use harpconvert but I am always getting an error saying " syntax error, unexpected $undefined, expecting $end" but I can not figure out whats the problem. I wrote: harpconvert -a ‘bin_spatial(-5,279,0.1,111,451,0.1)’ E:\S5P_NRTI_L2__O3_____20181101T064128_20181101T064628_05447_01_010102_20181101T072721.nc l3.nc
(I tried different paths, with and without space between the values and so on)
I also tried to use your example file, but near real-time products are not in the S5phub. Do you have any suggestions?

Thank you very much

[svniemeijer]

Dear Anna,

First, I assume you are trying this using a Windows Command Prompt? If so, you should not use single quotes around the operation. Using proper quoting in a Windows Command Prompt is rather tricky. In this case though you can omit the quotes, or use a double-quote character.

Second, you seem to have swapped the arguments to the bin_spatial function. The offset and length values need to be swapped.

Finally, your longitude offset seems to high. You probably want to start at e.g. 86 deg for that specific product.

So your final command (on Windows) should look like:

harpconvert -a bin_spatial(279,-5,0.1,451,86,0.1) E:\S5P_NRTI_L2__O3_____20181101T064128_20181101T064628_05447_01_010102_20181101T072721.nc l3.nc

[ABerni]

Dear sniemeijer,

Thank you for your quick response. That was very helpful and the command is working now. Nevertheless, it seems like I have problems with the coordinates. I read the help from https://cdn.rawgit.com/stcorp/harp/1.3/doc/html/operations.html and understood it like this: to get the first value I have to read the rows and colums from the file +1 (means my file has an extent of 278 and 450 so I use 279 and 451). The second value is the coordinate in lat and long. I used R to read out the min coordinates from the metadata file (“PRODUCT/latidue” and “PRODUCT/longitude”). And the third value is the 0.1 is the step size I am counting the coordinates per pixel. Is it a problem that my sample file is at the equator? Even after projecting it to WGS84 the file is way to big and not at the right location.

Knd regards

[svniemeijer]

Dear Anna,

The pixels in the L2 product do not form a regular grid.
A ‘bin_spatial’ operation is what can turn the data into a regular lat/lon grid, but you should then determine yourself what grid resolution to use. This is usually application specific.

A resolution of 0.025x0.025 deg would be about the same magnitude as the S5P resolution (around the equator; near the poles you can use a much coarser grid).

Once you have chosen your resolution, can you determine the size of your grid based on the min/max lat/lon of the file you have. e.g. -5…18 deg at 0.025 -> 920 grid cells. 86…115 deg at 0.025 -> 1160 grid cells.

harpconvert -a bin_spatial(920,-5,0.025,1160,86,0.025) E:\S5P_NRTI_L2__O3_____20181101T064128_20181101T064628_05447_01_010102_20181101T072721.nc l3.nc