Best practice to convert and reproject Sentinel-3 radiances to reflectance

Optical Toolbox
#8

Which MODIS data do you use?

#9

I use the MOD09GA collection 6.

#10

Hi, jumping onto this thread - how would I implement the radiance to reflectance conversion with snappy? Thanks, Sam

#11

Actually this is already implemented. You just need to call

olci_refl_Product = GPF.createProduct(‘Rad2Refl’, parameters, olci_rad_Product)

Or you implement you own operator if you know a better way. You know how this can be down.

The equation is pretty simple.
it is implemented here:


Respectivly here:

#12

Great thanks Marco

#13

Good morning,

I seem to have a problem with implementing this conversion. I have written this sort of script:

Script

import numpy as np
import snappy
from snappy import GPF
import jpy

# 1. Open Product
s3prod = snappy.ProductIO.readProduct(product_path)

# 2. Set parameters
HashMap = jpy.get_type('java.util.HashMap')
params = HashMap()
params.put("sensor", "OLCI")
params.put("conversionMode", "RAD_TO_REFL")
params.put("copyNonSpectralBands", "false")

# 3. Run conversion
s3prod_refl = GPF.createProduct("Rad2Refl", params, s3prod)

Testing

To test this I took a test product and compared the results of my script and running the SNAP desktop implementation. I got the values as:

# Get pixel values of Oa01 band
h = s3prod_refl.getSceneRasterHeight()
w = s3prod_refl.getSceneRasterWidth()
pixel_values = np.zeros(w * h, np.float64)
obj = s3prod_refl.getBand("Oa01_reflectance")
obj.readPixels(0, 0, w, h, pixel_values)
pixel_values.shape = h, w

#Results

I seem to have some difference between the pixel values.

For example, for test product: S3A_OL_1_EFR____20161023T100950_20161023T101250_20161023T120602_0179_010_122_1979_SVL_O_NR_002.SEN3

I find things like:

  • values from my script (blocks of the same value):
    pixel_values[0, 126] = 0.4881261…
    pixel_values[0, 127] = 0.4881261…
    pixel_values[0, 128] = 0.4881261…
    pixel_values[0, 129] = 0.4881261…

  • same values from snap desktop:
    desktop values[0, 126] = 0.49819
    desktop values[0, 127] = 0.49814
    desktop values[0, 128] = 0.49809
    desktop values[0, 129] = 0.49804

Any ideas what could be the problem here?

Thanks,

Sam

#14

Can you test something for me, Sam?

Please copy the file
C:\Users\<User_Home>\AppData\Roaming\SNAP\modules\org-esa-s3tbx-s3tbx-rad2refl.jar
to <INSTALL_DIR>\snap\modules\ and replace the existing one. Maybe you create a backup of this file first.

If I’m right the values will be the same as in the desktop.
I fear that the updates are not correctly considered by snappy. And that it is not deterministic which version is used.

#15

Do you mean:
C:\Users\<User_Home>\AppData\Roaming\SNAP\modules\org-esa-s3tbx-s3tbx-rad2refl.jar
to <INSTALL_DIR>\s3tbx\modules where the existing copy is?

If I put it in s3tbx\modules now on the line:
obj.readPixels(0, 0, w, h, pixel_values)

get the error:
RuntimeError: java.lang.RuntimeException: Cannot construct DataBuffer.

EDIT: Correction this now seems to work!

Confusing (wrong?) band names for subset with snappy
#16

Ahm yeah, in the INSTALL_DIR it is in the s3tbx directory. You’re right.

That it now works is not good news. But I expected it. This means that modules updates are not considered by snappy. Or at least not correct. It might be random from run to run which module is used. Either the old 5.0 module or the latest module e.g. 5.0.7. This means snappy is not reliable.

I consider this as critical issue for the next release.

Resampling in SNAP Desktop and snappy
#17

Okay, thanks for the information Marco. Will await the future release

#18

FYI after testing this doesn’t seem to be a problem in the linux installation of snap and snappy. Best, A.

#19

Thanks for this info.
Did you try several times? because it can happen that it sometimes works and sometimes not.
Actually there should be no difference between Windows and Linux.

#20

Marco,

In the end I managed to create a graph to convert OLCI TOA radiances to TOA reflectances using a graph. As you suggested I used the per-pixel geocoding instead if the tie-points, therefore I’ve got the full size latitude and longitude bands. Now I need additional information at full resolution, for instance, the view and Sun angles, I need per-pixel data, as in the latitude and longitude bands.

I checked the full resolution latitude and longitude and the tie-point and the interpolation it doesn’t seem to be linear, is it?. My question is, how do you create the full resolution bands from the tie-point bands?

Thanks again for your help.

Gerardo

#21

Actually the interpolation is linear between the tie-points. Why do you think it’s not?

#22

I was checking the TP_latitude first two tie points in IPython:

In [1]: import osgeo.gdal as gdal

In [2]: d = gdal.Open( 'TP_latitude.img' )

In [3]: lat = d.ReadAsArray()

In [4]: lat.shape
Out[4]: (4091, 77)

In [5]: lat[0,0]
Out[5]: 52.672676

In [6]: lat[0,1]
Out[8]: 52.658318

In the metadata, the STEP_X is 64:

<Tie_Point_Grid_Info>
    <TIE_POINT_GRID_INDEX>0</TIE_POINT_GRID_INDEX>
    <TIE_POINT_DESCRIPTION>Latitude</TIE_POINT_DESCRIPTION>
    <PHYSICAL_UNIT>degrees_north</PHYSICAL_UNIT>
    <TIE_POINT_GRID_NAME>TP_latitude</TIE_POINT_GRID_NAME>
    <DATA_TYPE>float32</DATA_TYPE>
    <NCOLS>77</NCOLS>
    <NROWS>4091</NROWS>
    <OFFSET_X>0.0</OFFSET_X>
    <OFFSET_Y>0.0</OFFSET_Y>
    <STEP_X>64.0</STEP_X>
    <STEP_Y>1.0</STEP_Y>
    <CYCLIC>true</CYCLIC>
</Tie_Point_Grid_Info>

Therefore when I did the Rad2Ref using the per-pixel geocoding, I was expecting equidistant latitude values between the two tie points in the full resolution latitude raster:

In [1]: import osgeo.gdal as gdal

In [2]: d = gdal.Open( 'latitude.tif' )

In [3]: lat_full = d.ReadAsArray()

In [4]: lat_full.shape
Out[4]: (4091, 4865)

In [5]: lat_full[0, 0]
Out[5]: 52.672672

In [6]: lat_full[0,64]
Out[6]: 52.658318

In [7]: np.diff( lat_full[0,0:65] )
Out[7]: 
array([-0.00022507, -0.00022125, -0.00022507, -0.00022507, -0.00022507,
       -0.00022125, -0.00022507, -0.00022507, -0.00022507, -0.00022507,
       -0.00022507, -0.00022125, -0.00022507, -0.00022507, -0.00022507,
       -0.00022507, -0.00022125, -0.00022507, -0.00022507, -0.00022507,
       -0.00022507, -0.00022125, -0.00022507, -0.00022507, -0.00022507,
       -0.00022507, -0.00022125, -0.00022507, -0.00022507, -0.00022507,
       -0.00022507, -0.00022125, -0.00022507, -0.00022507, -0.00022507,
       -0.00022507, -0.00022507, -0.00022125, -0.00022507, -0.00022507,
       -0.00022507, -0.00022507, -0.00022125, -0.00022507, -0.00022507,
       -0.00022507, -0.00022507, -0.00022507, -0.00022125, -0.00022507,
       -0.00022507, -0.00022507, -0.00022125, -0.00022507, -0.00022507,
       -0.00022507, -0.00022507, -0.00022507, -0.00022125, -0.00022507,
       -0.00022507, -0.00022507, -0.00022507, -0.00022125], dtype=float32)

As you can see the distance is not always the same, any idea why is that?, if is a linear interpolation the distance between points should be equidistant.

Thanks for your help.
Gerardo

#23

The difference is only at the 6th digit. So it is probably caused by the floating point precision.

#24

Hey Marpet,

is there a way to set the pixel based geocoding without user interface?

thanks!

1 Like
#25

Dear marpet,

if there is not a way to set these options from command line maybe it would be nice to add them among the reprojection option for s3, what do you think about?

thanks for considering this point!

Carlo

#26

There is a way. You can either create a s3tbx.properties file in the etc folder of the SNAP installation directory and add the property s3tbx.reader.olci.pixelGeoCoding=true.
You can specify this also on the command line. Then add -Ds3tbx.reader.olci.pixelGeoCoding=true to the call.

Similar properties are:
s3tbx.reader.slstrl1b.pixelGeoCodings
s3tbx.reader.meris.pixelGeoCoding (for upcoming MERIS 4th reprocessing in SAFE format)

1 Like
Manually georeference offset imagery in SNAP
#27

Solved!
thanks!

C