Extracting pixel values Sentinel-2a L2A with Snap or Snappy

development python
#1

Hey guys, I’m working with Sentinel-2a with level of preprocessing of level L2A and resampled all bands to 10m. I need to extract the bands values for some lat and lon coordinates that I have. I try this export-data-to-excel-from-snap, but I was wondering if I can just give some lat and lon coordinates and extract the values just for these coordinates instead of extract a whole area. By the way, I’m working with Snappy so I would like to do it not just in Snap Desktop, but also in Python with snappy. I really would appreciated any help, since I’m still pretty new using Snap.

Thanks in advance.

#2

Take a look at the PixEx tool (Raster / Export / Extract Pixel Values). You can let it run on a batch products from SNAP Desktop, and because it is an operator you can invoke it also from the command line and from Python.

#3

Thanks!!

It works in Snap Desktop, I tried with Snappy, but the file is never created. I already verified that the coordinates are inside the image. Do you know what could be happening? I modified the code from one this post extracting-single-pixel-information

import pandas as pd
import numpy as np
import traceback

from snappy import Product
from snappy import ProductIO
from snappy import ProductUtils
from snappy import WKTReader
from snappy import HashMap
from snappy import GPF
from snappy import jpy
from snappy import GeoPos
from snappy import PixelPos

path_to_sentinel_data = "S2A_MSIL2A_20160110T152632_N0201_R025_T18NXM_20160110T153105_subsetoutput2.dim"
product_subset = ProductIO.readProduct(path_to_sentinel_data)
lat, lon = (5.9147, 73.5136)


my_coordinates=jpy.array('org.esa.snap.pixex.Coordinate', 1)
Coord = jpy.get_type('org.esa.snap.pixex.Coordinate', 1)
my_coordinates[0] = Coord('bin1',lat, lon, None)
print(my_coordinates)
parameters = HashMap()
parameters.put('PexportBands', 1)
parameters.put('PexportExpressionResult', 0)
parameters.put('PexportMasks', 0)
parameters.put('PexportTiePoints', 0)
parameters.put('PoutputDir', 'C:\Users\camil\Desktop\output')

parameters.put('coordinates', my_coordinates)
c=GPF.createProduct('PixEx', parameters, product_subset)
c
[Lorg.esa.snap.pixex.Coordinate;@5c8ab9de

Out[45]:

org.esa.snap.core.datamodel.Product(objectRef=0x000000003CB9B328)
#4

Apart from the marpet propose solution, I come up with this solution, I suppose it could be of use for those who are trying to get the bands data into a pandas dataframe. :smiley:

import pandas as pd
import numpy as np
import traceback

from snappy import ProductIO
from snappy import GeoPos
from snappy import PixelPos

path_to_sentinel_data = "S2A_MSIL2A_20160110T152632_N0201_R025_T18NXM_20160110T153105_subsetoutput_ndvi.dim"
product_subset = ProductIO.readProduct(path_to_sentinel_data)

bands_names = list(product_subset.getBandNames())[1:13]
gc = product_subset.getSceneGeoCoding()
lat, lon = (5.914702982, -73.51362747)
print(bands_names)

cols = ['lat', 'lon', 'X', 'Y']
cols.extend(bands_names)
print(cols)

df_bands = pd.DataFrame(columns=cols)

pixel_pos = gc.getPixelPos(GeoPos(lat, lon), None)
data = list()
print("(lat, lon) -→ (X, Y) : (%s, %s) -→ (%s, %s)" % (lat, lon, int(pixel_pos.x), int(pixel_pos.y)))

data.append(lat)
data.append(lon)
data.append(int(pixel_pos.x))
data.append(int(pixel_pos.y))

for i, band_name in enumerate(bands_names):
    temp_band = product_subset.getBand(band_name)
    width, height = temp_band.getRasterWidth(), temp_band.getRasterHeight()
    try:
        tmp = np.zeros(1)
        temp_band.readPixels(int(pixel_pos.x), int(pixel_pos.y), 1, 1, tmp)
        data.append(tmp[0])
    except Exception as e:
        print(band_name)
        print(width, height)
        print(int(pixel_pos.x), int(pixel_pos.y))
        print(e)
        traceback.print_exc()

df_bands.loc[df_bands.shape[0]] = data
df_bands

Out[]

['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B9', 'B11', 'B12']
['lat', 'lon', 'X', 'Y', 'B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B9', 'B11', 'B12']
(lat, lon) -→ (X, Y) : (5.914702982, -73.51362747) -→ (2071, 3162)

image
image780×62 2.37 KB

1 Like
#5

One issue with the code is that the parameters should not start with ‘P’ and I’m not sure if the one and zero are properly converted to True and False.

1 Like
#6

Jep, your solution is working too. :slight_smile:

1 Like
#7

I already made the changes, but it still doesn’t generate any file. :confused:

my_coordinates=jpy.array('org.esa.snap.pixex.Coordinate', 1)
Coord = jpy.get_type('org.esa.snap.pixex.Coordinate', 1)
my_coordinates[0] = Coord('bin1',lat, lon, None)
print(my_coordinates)
parameters = HashMap()
parameters.put('exportBands', True)
parameters.put('exportExpressionResult', False)
parameters.put('exportMasks', False)
parameters.put('exportTiePoints', False)
parameters.put('outputDir', 'C:\Users\camil\Desktop\output')

parameters.put('coordinates', my_coordinates)
c=GPF.createProduct('PixEx', parameters, product_subset)
c

However, it’s pretty weird since I’m using the same coordinates that with the other code.

#8

Yes, PixEx behaves sometimes like a diva…

It could be related to this issue: https://senbox.atlassian.net/browse/SNAP-484

#9

hahaha I see. Anyway, thanks!!