# Reading Sentinel-3 data
## DATA source
Sentinel-3 L1 EFR (OL\_1\_EFR) downloaded from [Copernicus Open Access Hub](https://scihub.copernicus.eu/dhus/#/home).
## 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](https://guides.geospatial.bas.ac.uk/working-with-scientific-data-formats-in-gis/intro-notes#configuring-environment-for-jupyter-notebook), savethe notebook to your working directory, and run jupyter notebook from this directory.
{% file src="" %}
### 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
```python
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
```
```python
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.
```python
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.
```python
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 ](https://gdal.org/tutorials/geotransforms_tut.html)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
```python
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
```python
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
```
{% code overflow="wrap" %}
```python
save_file_multiband(rgb_array,'S3_OLCI_2023-03-16_RGB.tif', gcps, epsg=4326, no_data=-32768)
```
{% endcode %}
### Reproject
{% code overflow="wrap" %}
```python
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')
```
{% endcode %}
### Check the output in QGIS
Every file is read as 2D array with no geotransformation.
```python
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.
```python
geo = xr.open_dataset('geo_coordinates.nc')
lon = geo.longitude.values
lat = geo.latitude.values
```
