☕Reading Sentinel-3 data

Sentinel-3 data from Copernicus Open Access Hub could be easily read in SNAP.

DATA source

Sentinel-3 L1 EFR (OL_1_EFR) downloaded from Copernicus Open Access Hub.

SNAP

To open data: File > Import > Optical Sensors > Sentinel-3 >Sentinel-3. Select *.xml file

To visualize data: Window > New RGB window > Select one of the profiles

To export data: File > Export > GeoTIFF / BigTIFF. Click Subset... button, change the spatial subset and select at least 3 bands (e.g. 8-6-4 for RGB). Export file.

SNAP window for data expoort

QGIS

Open the output file in QGIS. The georeferencing is wrong. Original *.nc files are not georeferenced either.

SNAP output in QGIS

Python

The full notebook download is below. To run it, set up conda environment, savethe notebook to your working directory, and run jupyter notebook from this directory.

341KB

Sentinel3_NETCDF_exp.ipynb

Read image data

GDAL can read NetCDF files as arrays. We read all the bands and store them to dict with band names.

import numpy as np
from osgeo import gdal, osr
import xarray as xr
import matplotlib.pyplot as plt
from glob import glob

band_files = sorted(glob('Oa*radiance.nc'))
bands_dict = {}
for n, file in enumerate(band_files):
    band_ds = xr.open_dataset(file)
    band = band_ds[file.replace('.nc', '')].values
    bands_dict[str(n+1).zfill(2)] = band

plt.imshow(bands_dict['04'],
           cmap='gray',
           vmin=np.nanpercentile(bands_dict['04'], 2),
           vmax=np.nanpercentile(bands_dict['04'], 98))
plt.axis('off')
plt.colorbar();

Every file is read as 2D array with no geotransformation.

rgb_array = np.stack([
    bands_dict['08'],
    bands_dict['06'],
    bands_dict['04']
], axis=2)

Sentinel-3 georeference info is stored in the geo_coordinates.nc file, which contains 2D latitude and longitude grids.

geo = xr.open_dataset('geo_coordinates.nc')
lon = geo.longitude.values
lat = geo.latitude.values

You can see that the coordinates gradient is not parallel to the image sides. This means, that it is not possible to produce a GDAL geotransform line. For this case, we can create a GDAL GCP object to warp the image with it.

GCP object consists of list of geocoordinates with corresponding pixel coordinates

gcps = [] #empty list
num = 20 #number of points along one of the axes
for row in np.linspace(0, rgb_array.shape[0]-num, num):
    for col in np.linspace(0, rgb_array.shape[1]-num, num):
        row = int(row)
        col = int(col)
        gcp = gdal.GCP(lon[row, col], lat[row, col], 0, col, row)
        gcps.append(gcp)

Custom function for saving array to GeoTIFF

def save_file_multiband(input_array, out_file, gcps, epsg=4326 , no_data=-32768):
    
    """Save numpy array to multiband GTiff file"""
    
    driver = gdal.GetDriverByName('GTiff')
    
    sr = osr.SpatialReference()
    sr.ImportFromEPSG(epsg) #create projection reference from EPSG code
    
    dataset = driver.Create(out_file,
                        input_array.shape[1],
                        input_array.shape[0],
                        input_array.shape[2],
                        eType=gdal.GDT_Float32)
    
    for band in range(input_array.shape[2]):
        outband = dataset.GetRasterBand(band + 1)
        outband.SetNoDataValue(no_data)
        outband.WriteArray(input_array[:, :, band])
        
    dataset.SetGCPs(gcps, sr.ExportToWkt())        
    dataset = None

save_file_multiband(rgb_array,'S3_OLCI_2023-03-16_RGB.tif', gcps, epsg=4326, no_data=-32768)

Reproject

gdal.Warp('/mnt/c/DATA/DropIn/S3/S3_OLCI_{}_RGB_3031.tif'.format(date),
          '/mnt/c/DATA/DropIn/S3/S3_OLCI_{}_RGB.tif'.format(date),
          dstSRS='EPSG:3031')

Check the output in QGIS

The output is georeferenced correctly