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.
QGIS
Open the output file in QGIS. The georeferencing is wrong. Original *.nc files are not georeferenced either.
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.
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
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
