Gpt outputs geotiff with corrupted bands after BandSelect

Dear SNAP team,

I get strange output from gpt with the following graph test.xml (4.7 KB) and gpt call:

gpt test.xml -PSLSTRsource=“S3A_SL_1_RBT____20180315T103348_20180315T103648_20181006T082919_0179_029_051______LR1_R_NT_003.SEN3” -PtargetFolder="." -Ds3tbx.reader.olci.pixelGeoCoding=true -Ds3tbx.reader.slstrl1b.pixelGeoCodings=true -Dsnap.log.level=ERROR -e

which runs without problem:

INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: GDAL 2.4.1 found on system. JNI driver will be used.
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: Installed GDAL 2.4.1 set to be used by SNAP.
INFO: org.esa.snap.core.util.EngineVersionCheckActivator: Please check regularly for new updates for the best SNAP experience.
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: Installed GDAL 2.4.1 set to be used by SNAP.
Executing processing graph
INFO: org.hsqldb.persist.Logger: dataFileCache open start
…10%…20%…30%…40%…50%…60%…70%…80%…90% done.

The geotiff contains 67 bands. Many of them are empty and after inspection, I am guessing that the processed SLSTR data is in the last band.

Band 67:

Description of the geotiff in QGIS:

Name r_TOA_S1
Path [/application/pi/Tests/r_TOA_S1.tif](file:///application/pi/Tests/r_TOA_S1.tif)
CRS EPSG:3413 - WGS 84 / NSIDC Sea Ice Polar Stereographic North - Projected
Extent 360060.8351618153974414,-1507438.5418066726997495 : 2242688.6933637140318751,293635.9603251670487225
Unit meters
Width 3001
Height 2871
Data type Float32 - Thirty two bit floating point
GDAL Driver Description GTiff
GDAL Driver Metadata GeoTIFF
Dataset Description /application/pi/Tests/r_TOA_S1.tif
Compression
Band 1 * STATISTICS_MAXIMUM=128
* STATISTICS_MEAN=2.0650142046
* STATISTICS_MINIMUM=0
* STATISTICS_STDDEV=15.17489063947
Band 2 * STATISTICS_MAXIMUM=128
* STATISTICS_MEAN=3.7547184724564
* STATISTICS_MINIMUM=0
* STATISTICS_STDDEV=20.927269505439
Band 3 * STATISTICS_MAXIMUM=128
* STATISTICS_MEAN=2.0643456709138
* STATISTICS_MINIMUM=0
* STATISTICS_STDDEV=15.174256142999
Band 4 * STATISTICS_MAXIMUM=128
* STATISTICS_MEAN=3.7542504988758
* STATISTICS_MINIMUM=0
* STATISTICS_STDDEV=20.926979342654
Band 5
Band 6
Band 7
Band 8
Band 9
Band 10
Band 11
Band 12
Band 13
Band 14
Band 15
Band 16
Band 17
Band 18
Band 19
Band 20
Band 21
Band 22
Band 23
Band 24
Band 25
Band 26
Band 27
Band 28
Band 29
Band 30
Band 31
Band 32
Band 33
Band 34
Band 35
Band 36
Band 37
Band 38
Band 39
Band 40
Band 41
Band 42
Band 43
Band 44
Band 45
Band 46
Band 47
Band 48
Band 49
Band 50
Band 51
Band 52
Band 53
Band 54
Band 55
Band 56
Band 57
Band 58
Band 59
Band 60
Band 61
Band 62
Band 63
Band 64
Band 65
Band 66
Band 67 * STATISTICS_MAXIMUM=1.6927000284195
* STATISTICS_MEAN=0.40414608065538
* STATISTICS_MINIMUM=0.058100000023842
* STATISTICS_STDDEV=0.17680351524881
More information * AREA_OR_POINT=Area
* TIFFTAG_IMAGEDESCRIPTION=r_TOA_S1
* TIFFTAG_RESOLUTIONUNIT=1 (unitless)
* TIFFTAG_XRESOLUTION=1
* TIFFTAG_YRESOLUTION=1
Dimensions X: 3001 Y: 2871 Bands: 67
Origin 360061,293636
Pixel Size 627.3335082312224813,-627.3335082312224813

This issue occurred with SNAP7, remained when updating to SNAP8 and I may have seen similar outputs when processing OLCI scenes after similar geotiff generation.

Thanks for your help!

You shouldn’t have to guess at the contents of a band. BEAM Dimap or NETCDF-CF should preserve more of the metadata than GeoTIFF’s, and can be imported to many (current) 3rd party applications.

It could be helpful to know if the problem is limited to GeoTIFF’s.

Thank you for your quick response.
The content of band 67 in the output geotiff is indeed slightly different (different shape of histograms) than what I get when calculating S1 reflectances in SNAP.

Band 67 of r_TOA_S1:

S1_reflectance_an computed in SNAP:
image

But before we go into how these rasters differ, I’d be interested to know if it is normal that I get 67 bands as output and if the problem is reproducible by others.

Link to the SLSTR scene in the example: https://scihub.copernicus.eu/dhus/odata/v1/Products(‘53516775-7347-46cb-a505-be070182be1d’)/$value

After replicating the issue on different machines, I can confirm that this problem appears with the installation of SNAP8. It occurs with both SLSTR and OLCI data when using Rad2Refl -> Reproject -> BandSelect -> Write as geotiff.

The additional bands seem to be tie-point grids of auxiliary data.
I can for example reduce the number of bands outputed by changing

  <node id="Rad2Refl_OLCI">
    <operator>Rad2Refl</operator>
    <sources><source>Read_OLCI</source></sources>
    <parameters>
      <sensor>OLCI</sensor>
      <copyTiePointGrids>true</copyTiePointGrids>
      <copyFlagBandsAndMasks>true</copyFlagBandsAndMasks>
      <copyNonSpectralBands>true</copyNonSpectralBands>

to

  <node id="Rad2Refl_OLCI">
    <operator>Rad2Refl</operator>
    <sources><source>Read_OLCI</source></sources>
    <parameters>
      <sensor>OLCI</sensor>
      <copyTiePointGrids>true</copyTiePointGrids>
      <copyFlagBandsAndMasks>false</copyFlagBandsAndMasks>
      <copyNonSpectralBands>false</copyNonSpectralBands>

However that makes me unable to filter the output of Rad2Refl with quality flags or extract auxiliary information (altitude, ozone… etc) from the same node.

It seems that there has been a change, in SNAP8, in the way these auxiliary data are dealt with when written as geotiffs. I hope an easy fix is possible.

We remain stalled for processing 2020 Sentinel-3 data.

In recent days, Adrien and Baptiste independently tested SNAP 8, Adrien got “Error: SPI not found for operator 'BandSelect’”

On a different system Baptiste had an unexpected extraction of auxiliary data as real data and writes them as additional bands in geotiffs.

We’re hoping others can illuminate the situation!

best,

Jason

Hi Baptiste, Hi Jason.

The additional bands in the geotiff are probably caused by one fix:
[SNAP-1232] Masks not correctly copied by ProductUtils.copyProductNodes method - JIRA (atlassian.net)
The masks band were just not included and should have been copied in earlier version.
I first thought you see the results of the following fix:
[SNAP-1152] Subsetting in graph does not consider tie-point grids correctly - JIRA (atlassian.net)
But it is more likely that the first solved issue.
To solve this, you can add a subset operator to your graph which leaves out the flag bands.

Regarding the "Error: SPI not found for operator 'BandSelect’ problem.
Was the S1TBX installed when this occurred?
Because this operator is provided by this toolbox.

Thanks for the reply Marco.

Can you give a bit more details? I thought Subset was to clip a product for a given area.
Which parameter should I use to keep the full scene? and to leave the masks out?

Since I need some of the masks for BandMath before I write each band to geotiff, can I do something like:
Read > Reproject > BandMath (where I use the masks to process a given band) > Subset (to drop the masks) > Write

Subset can do both regional subset and band subset, too.
You can use sourceBands to define the bands to include and tiePointGridNames for tie-point grids.
Leave region and geoRegion empty.
See also help on gpt:

Sounds reasonable. Should work.

1 Like

Thanks!
I manage to extract the desired band without the masks. But it seems that some of the masks are applied to the band automatically. I think it is the flag called “quality_flags_duplicated”.
How can I prevent this behaviour from happening?

<graph id="S3">
  <version>1.0</version>

  <node id="Read_OLCI">
    <operator>Read</operator>
    <parameters>
      <file>${OLCIsource}/xfdumanifest.xml</file>
      <formatName>Sen3</formatName>
    </parameters>
  </node>

  <node id="Rad2Refl_OLCI">
    <operator>Rad2Refl</operator>
    <sources><source>Read_OLCI</source></sources>
    <parameters>
      <sensor>OLCI</sensor>
      <copyTiePointGrids>true</copyTiePointGrids>
      <copyFlagBandsAndMasks>true</copyFlagBandsAndMasks>
      <copyNonSpectralBands>true</copyNonSpectralBands>
    </parameters>
  </node>

  <node id="Reproject_OLCI">
    <operator>Reproject</operator>
    <sources><source>Rad2Refl_OLCI</source></sources>
    <parameters>
      <crs>EPSG:3413</crs>
      <resampling>Nearest</resampling>
      <noDataValue>NaN</noDataValue>
      <includeTiePointGrids>true</includeTiePointGrids>
    </parameters>
  </node>
    
  <node id="Subset_OZA">
    <operator>Subset</operator>
    <sources><source>Reproject_OLCI</source></sources>
    <parameters>
	<sourceBands>OZA</sourceBands>
	<copyMetadata>false</copyMetadata>
    </parameters>
  </node>
  
  <node id="Write_OZA">
    <operator>Write</operator>
    <sources><sourceProduct>Subset_OZA</sourceProduct></sources>
    <parameters><file>${targetFolder}/OZA_x.tif</file><formatName>GeoTIFF</formatName></parameters>
  </node>

</graph>

I think your image does not show the OZA band. It seems there is a coastline visible. Which should not be the case for the observation zenith angle.
I’ve used your graph and got the result on the right. No NaN values visisble.

Hi Marco,

Sorry, my illustration was not very useful. Here is a better one:

Have you run the graph through gpt or using SNAP’s graph builder?

When I run the previous graph using GPT, I get the output above, with NaN.

When I load the graph into SNAP’s graph builder, I have an error related to the EPSG 3413 CRS: “This “AxisDirection” object is too complex for WKT syntax”

When I choose another CRS (e.g. EPSG 4326) then SNAP produces a gap-free geotiff:

But using GPT and either i) Reproject with EPSG 4326; or ii) no Reproject (feeding directly Rad2Refl to Subset), still replaces many pixels by NaNs as illustrated above.

Thanks for your help

I simply used gpt from the command line.
When I only use the Reprojection operator, it also works for me. Even if the scene is not located north.

But I was able to reproduce the AxisDirection issue in the GraphBuilder
So it seems, the GrqaphBuilder has an issue here.

@lveci can you check this?

I found the issue!
I was using gpt with the option " -Ds3tbx.reader.olci.pixelGeoCoding=true"
If I remove it, the output is normal: no NaN pixels.

Any idea why this command is not valid any more?
Thanks a lot for your help!

Ah, that’s good.
Actually, it is still valid, but there are currently issues at high latitudes. We are already working on it. That’s why it wasn’t visible in my scene. It wasn’t located so far in the north.

The location accuracy might be a bit lower now in your case. But I guess still good enough.
ow the tie-points are used.

To summarize the solutions:

Multiple bands in outputs when using BandSelect
This is the expected behaviour of BandSelect. It outputs masks and flags as extra bands.
Use Subset instead.
Example:

  <node id="Subset_OZA">
    <operator>Subset</operator>
    <sources><source>Reproject_OLCI</source></sources>
    <parameters>
	<sourceBands>OZA</sourceBands>
	<copyMetadata>false</copyMetadata>
    </parameters>
  </node>

NaN in outputs after subset
Using pixel geocoding replaces some pixels by NaNs. This appears only when using the option:
" -Ds3tbx.reader.olci.pixelGeoCoding=true" in gpt

Error when reprojecting in EPSG 3413 in SNAP’s graphbuilder

Under investigation.

1 Like