# 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.

<figure><img src="/files/oiwhT5cyQ1c2OoskzHSS" alt=""><figcaption><p>SNAP window for data expoort</p></figcaption></figure>

## QGIS

Open the output file in QGIS. The georeferencing is wrong. Original \*.nc files are not georeferenced either.&#x20;

<figure><img src="/files/ZPPWocLR0409spGcS4CK" alt=""><figcaption><p>SNAP output in QGIS</p></figcaption></figure>

## Python

The full notebook download is below. To run it, [set up conda environment](/working-with-scientific-data-formats-in-gis/intro-notes.md#configuring-environment-for-jupyter-notebook), savethe notebook to your working directory, and run jupyter notebook from this directory.&#x20;

{% file src="/files/1YIxwk2yVb1lKSV6i35f" %}

### Read image data

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

<pre class="language-python"><code class="lang-python">import numpy as np
<strong>from osgeo import gdal, osr
</strong>import xarray as xr
import matplotlib.pyplot as plt
from glob import glob
</code></pre>

```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();
```

<img src="/files/fDkxZ0jIVTYR7zRN7FOr" alt="" data-size="original">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.&#x20;

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

<figure><img src="/files/e7JyWDRyIdymuJ1fQJdD" alt=""><figcaption></figcaption></figure>

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.

<figure><img src="/files/uj17yBh5rX0BZl4j0NTf" alt=""><figcaption></figcaption></figure>

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

<figure><img src="/files/Ur16uM4KYOy5hNtDI5jn" alt=""><figcaption><p>The output is georeferenced correctly</p></figcaption></figure>


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://guides.geospatial.bas.ac.uk/working-with-scientific-data-formats-in-gis/reading-sentinel-3-data.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.