5

I'm trying to convert Sentinel-3 data to GeoTiff file. The Sentinel-3 data comes as separate nc files (netcdf): a file for each spectral band and a nc file for qualityFlags, geo_coordinates etc' and an xml file. The official why to read Sentinel-3 data is to use the SNAP software and loading to SNAP the xml file that directs SNAP to the different files. I would like to do it using a python script, but I can't figure how. I tried GDAL translate but it does not make the connection between the spectral band file and the coordinate file. I am able to read each nc file separately and to read the coordinates nc file as text/array.

How can I combine the spectral band data with the coordinate data to create a Geotiff?

My guess is that it is possible doing that by using GDAL but I can't find the way. Is anyone familiar with Sentinel-3 data? I'm using python 2.7 on win10, I have GDAL 2.1

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
user88484
  • 1,779
  • 4
  • 26
  • 39

2 Answers2

3

it's a bit old but here I am. Unfortunately GDAL does not fully support the NetCDF standard and Sentinel 3 netcdf are one of those unfortunate case.

The following is a simple script written in kotlin unsing the java bindings:

    println(" * Converting $prodName...")
    val wgs84 = SpatialReference()
    wgs84.ImportFromEPSG(4326)

    val lst = gdal.Open("NETCDF:$prodName/LST_in.nc:LST")
    val map = mapOf(
            "LINE_OFFSET" to "1", "LINE_STEP" to "1",
            "PIXEL_OFFSET" to "1", "PIXEL_STEP" to "1",
            "X_BAND" to "1", "X_DATASET" to "NETCDF:$prodName/geodetic_in.nc:longitude_in",
            "Y_BAND" to "1", "Y_DATASET" to "NETCDF:$prodName/geodetic_in.nc:latitude_in"
    )

    lst.SetMetadata(Hashtable(map), "GEOLOCATION")

    val warp = gdal.AutoCreateWarpedVRT(lst, wgs84.ExportToWkt())

    gdal.Translate("$prodName/lst.tif", warp, TranslateOptions( gdal.ParseCommandLine("-of gtiff -oo COMPRESS=LZW ") ) )

    println(" * Complete")

this perfectly works with Sentinel 1 OCN products, but with Sentinel 3 GDAL gives some strange errors while parsing the file:

Warning 1: NetCDF driver detected file type=5, but libnetcdf detected type=3
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Warning 1: NetCDF driver detected file type=5, but libnetcdf detected type=3
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Warning 1: NetCDF driver detected file type=5, but libnetcdf detected type=3
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
ERROR 1: Too many points (441 out of 441) failed to transform, unable to compute output bounds.
ERROR 1: IReadBlock failed at X offset 2, Y offset 0
ERROR 1: GetBlockRef failed at X block offset 2, Y block offset 0

After some throubleshooting it seems that it does not automatically apply the scale factor and offset of each variable, and so he have to deal with raw lat-lon coordinates and he does not know what to do with it.

I hope this will get fixed soon

fer.marino
  • 131
  • 3
  • This still seems to be the case with gdalinfo --version -> GDAL 2.3.2, released 2018/09/21. The same error message isn't exactly the same with I have (S-5P), but the type 5/3 mismatch is. One way of testing gdalinfo NETCDF:"<filename>.nc". – Veksi Oct 29 '18 at 15:28
  • I'll take the previous back, some valuable information at https://trac.osgeo.org/gdal/ticket/6856. Specifically the command gdalinfo --format netcdf to check which formats are supported. At least mine doesn't currently list support for the latest version. Something's amiss here. – Veksi Oct 29 '18 at 15:59
  • It appears http://www.gisinternals.com/ packages come with support to NC and NC2 whereas https://trac.osgeo.org/osgeo4w/ lists NC, NC2, NC4 and NC4C. – Veksi Oct 29 '18 at 16:19
  • 3
    This is not GDAL's fault. The product format for Sentinel 3 products is horrible. Yes, the data values are stored in netCDF, but the coordinate axes are in separate files and it is all just a bunch of files and metadata. Use ESA's SNAP software to reproject and convert to GeoTIFF for sanity. – j08lue Dec 21 '18 at 14:52
  • 2
    Oh, for sustained sanity, remember the per-pixel geocoding when you reproject Sentinel 3 products with SNAP. – j08lue Dec 21 '18 at 14:55
-1

You could use high level libraries such as EOReader to geocode directly your data:

from eoreader.reader import Reader
from eoreader.bands import SWIR_2 # which is S6 band

Open the Sentinel-3 product

s3_path = "S3B_SL_1_RBT____20191115T233722_20191115T234022_20191117T031722_0179_032_144_3420_LN2_O_NT_003.SEN3"

prod = Reader().open(s3_path)

Save SWIR band directly on disk

swir = prod.stack(SWIR_2, stack_path="your/path.tif")

More information on Sentinel-3 data are available in this notebook.

Disclaimer: I am the maintainer of EOReader.

Kadir Şahbaz
  • 76,800
  • 56
  • 247
  • 389
remi.braun
  • 423
  • 2
  • 6
  • 2
    You have posted this answer several times. In case this is your development or you are associated to it, please add a disclaimer. Otherwise this may be detected as spam. – MrXsquared Dec 22 '22 at 13:53
  • Oh yes I am the main maintainer ;) How should I add a disclaimer ? – remi.braun Dec 22 '22 at 13:56
  • Just a short sentence: "Disclaimer: I am the maintainer of EOReader" or similar. You can add it by editing your answer. – MrXsquared Dec 22 '22 at 13:57