Pixel extraction from many SENTINEL-3 products in SNAP


I am trying to extract pixel values (Raster->Export->Extract Pixel Values) from many different Sentinel OLCI .xml files (many source paths). I add my coordinates and I select my output directory and then I select extract. Then a successful message appears (The pixel extraction tool has run successfully and written the result files to C:…etc), but there is no an output file (.txt *file /prefix: pixEx ) in the output directory.

Can you help me?

Thank you a lot,

1 Like

This often happens because some contraints are not met. Either the allowed time distance too high or the coordinates are not located in the source products or the source products are rejected for some other reason.
Unfortunately the error reporting within the desktop application is not optimal. When you invoke it from the command line you should see the reasons why no pixels are extracted.

I’ve attached a template file where you can insert the parameters from the GUI.
PixExTemplate.xml (409 Bytes)

On the command line you can invoke it like this:

gpt PixExTemplate.xml

1 Like

Hello there @marpet, first of all thanks for making yourself available and dedicating your time helping us all over the forum.

I am having a similar issue, so I thought it would be ok to post in the same thread.

I have a region of interest in ESRI Shapefile format, I can subset my image using it, ok. But how do I extract all the values from the pixels within my shapefile? I don’t see any options related to it in the GUI…

I have managed to programmatically extract a single pixel by a given coordinate using python, but I am having a hard time manupulating the Sentinel-3 NetCDF files. Any directions would be of a great help on this subject.

1 Like

Hi, and thanks for the nice words.

The question is do you want the extracted data in an ASCII formator or do you still want an image format? It is also different if you want to do it on a batch of products or just a few in the desktop.

1 Like

And also thanks for the quick reply, you see? You deserve more than nice words for helping us with such complex tasks.

Both will work, so long I can get rid of all the data outside the shapefile/ROI or get all the data inside it.

My final goal is to have the value of every individual pixel within the green area for all available bands of a given Sentinel-3 WFR or EFR product, even if that data is invalid or under a cloud pixel.

My biggest struggle so far is that the individual Oa**_reflectance.nc files inside the image folder seems not to have geographical information attached to them, and that information is inside the geo_coordinates.nc

I have tried converting the WFR NetCDF images into a GeoTIFF/BigTIFF to make them easier to manipulate but the resulting stack was like 20Gb+ and I needed only to operate in the pixels within this green shapefile above…

(just for the sake of curiosity, this is part of the Manacapuru river in Brazil btw)

When you say ‘different’ do you mean the method? I’ve managed to do it on 209 WFR images on a desktop, but I was extracting all the band values for a single pixel.

You are correct, I don’t really know how computationally expensive it will be when I do it for all the pixels within the shapefile instead of only one pixel…

This is my python3 code so far, but it is not very efficient because it loads up in memory all the NetCDF data from all the S3 bands of a given image:

bands_dictionary is a list of numpy N-D-Arrays containing the NetCDF data for ALL the bands
lon and lat are the longitude and latitude dimensions inside de file geo_coordinates.nc
target_lon and target_lat are the point of interest to be extracted like: -60.014493, -3.158980

def get_point_data_in_bands(self, bands_dictionary, lon=None, lat=None, target_lon=None, target_lat=None):
    lat = lat[:, :, np.newaxis]
    lon = lon[:, :, np.newaxis]
    grid = np.concatenate([lat, lon], axis=2)
    vector = np.array([target_lat, target_lon]).reshape(1, 1, -1)
    subtraction = vector - grid
    dist = np.linalg.norm(subtraction, axis=2)
    result = np.where(dist == dist.min())
    target_x_y = result[0][0], result[1][0]

    rad_in_bands = []
    for b in bands_dictionary:
        rad = bands_dictionary[b][target_x_y]

    return target_x_y, rad_in_bands

The output target_x_y is just in case I want to re-use it later, and represents a tuple containing the absolute X,Y position ex: (3825, 726) in the image matrix (like SNAP in the Pin Manager tool):


And rad_in_bands is what I actually want, a list with the radiance of the pixel in each band e.g.:


I feel like I just did way too much code gymnastics to do something that probably has a better way of doing. And before I started writing more spaghetti code to do the same above for all the pixels inside a given ESRI shapefile:


I thought it would be smarter to ask you @marpet how to do it the correct way :sweat_smile:

Ok, sorry for flooding the thread but I have just found this answer from you:

And it is exactly what I need, now, how do I do it from the command line?

Good that you have found this way. That’s what I wanted to explain for the Desktop use case.

When you use python you can use snappy. This helps reading the data.


However, the simplest way to do on a batch of products is probably to use gpt. There are two wiki pages available which explain how to do it. They are linked here:

For your use case, you need 3 operators chained in a graph.
I’ve attached a template. shapefileExtraction.xml (1.5 KB)

In the first step, the data is tailored to the bounding box. Here I just added a random rectangle.
In the second step, the shapefile is imported and in the last step, the data outside of the shapefile is masked.

Unfortunately, this approach does not export only the pixels covered by your shapefile but always the bounding box. Especially when exporting to CSV/ASCII format the amount of data is huge.

1 Like

Thanks so much for all the tutorials and also for the xml template @marpet

through some trial and error I’ve managed to get snappy up and running system-wide alongside GDAL:

Ps. I’ve also had this SEVERE message in the console but after some searching I saw this guy saying it was ok:

Is there a way of making the gpt/gpf or snappy to simulate this behavior of the GUI?:

I did it in SNAP Desktop GUI and got this file: S3A_OL_2_WFR____20190309T141223_20190309T141523_20190310T211622_0179_042_167_3060_MAR_O_NT_002_manacapuru_Mask.txt (3.1 MB)

And I didn’t tested it to see where the lon/lat pixel data will fall on a map, was planning on doing it later today, but… My real question is: does that mean that the pixels that I extracted this way in the above file are actually everything inside the bounding box rectangle of the shapefile or only the ones actually inside the shapefile?

Sorry to keep bothering you guys with this, you already did such a great job by making all these great tools and data freely available to us.

Thanks again and keep safe! :raised_hands: :microbe:

gotta love this SNAP man… The points extracted by the GUI fall right inside the ROI/SHP file:

Now… so you can get rid of me @marpet, how do I replicate this for my 209 images?
Cheers and a huge thanks,

Unfortunately, this is not possible to do it on a stack of data.
Only what I already wrote regarding gpt is possible.
We have this improvement already on the todo list (SNAP-828).

What you could do and if it is applicable for your use case, you could create a list of coordinates and use the Pixel Extraction tool (usable on a stack of products even in the GUI).
This way you wouldn’t get all pixels covered by your region, but maybe enough?
The data exported with the mask pixel extraction could be a starting point for this list.