Wrong projection- Creating GeoTIFF from Sentinel 3 Land Surface Temperature (LST)

Question

I have been working on creating a Land Surface Temperature (LST) map from the Sentinel-3 LST product. After selecting a polygon region in the Mediterranean I have downloaded the data from the Copernicus Data Space website. The date of the satellite image is 2024-11-17. The data I have downloaded is in NetCDF format.

I wrote a script to convert the NetCDF file of LST into GeoTIFF format, but I seem to be facing an issue with the projection. I have mostly used xarray and rasterio libraries. After applying the projection and exporting the TIFF image, the output is not correctly located in QGIS (it is moved slightly to the left as can be seen on the map). The TIFF image is not in the right place. I also tried using EPSG:32634 when writing the TIFF image, as it is commonly used for this region in many other Sentinel products, but it did not yield the desired results. The code is attached below:

!pip install xarray netCDF4 rasterio

import xarray as xr
import rasterio
from rasterio.transform import from_bounds

# File paths
lst_file = "LST_in.nc"
geo_file = "geodetic_in.nc"

# Load the datasets
lst_data = xr.open_dataset(lst_file)
geo_data = xr.open_dataset(geo_file)

# Extract the LST variable and convert from Kelvin to Celsius
lst = lst_data['LST'] - 273.15  # Convert Kelvin to Celsius

# Print metadata to check the units
print(lst_data['LST'].attrs)


# Extract latitude and longitude from the geo dataset (assuming they are stored there)
longitude = geo_data['longitude_in']  # Replace with the actual variable name
latitude = geo_data['latitude_in']    # Replace with the actual variable name

# Define the transform using bounds and resolution
transform = from_bounds(
    west=longitude.min().values,  # Minimum longitude
    east=longitude.max().values,  # Maximum longitude
    south=latitude.min().values,  # Minimum latitude
    north=latitude.max().values,  # Maximum latitude
    width=lst.shape[1],           # Number of columns
    height=lst.shape[0],          # Number of rows
)

# Output file path
output_tiff = "LST_output_v7.tif"

# Write the data to a GeoTIFF
with rasterio.open(
    output_tiff,
    'w',
    driver='GTiff',
    height=lst.shape[0],
    width=lst.shape[1],
    count=1,  # Number of bands
    dtype=lst.dtype,
    crs="EPSG:4326", 
    transform=transform,
) as dst:
    dst.write(lst.values, 1)  # Write the first (and only) band

print(f"GeoTIFF saved to {output_tiff}")

What might be missing or what steps I need to follow to ensure the exported GeoTIFF is georeferenced correctly?

Answer 1

Sentinel-3 products don't have rectangular grids so you cannot use a classic transform.

Here is a thread with an example on how to geocode Sentinel-3 data with rioxarray (it's a bit complicated).

Maybe it's best to use directly specialized libraries such as pyresample.

1. Create a SwathDefinition from lat/lon arrays
2. Create your area definition with create_area_def
3. Use the XArrayBilinearResampler or XArrayResamplerNN to resample your data

-> Here is an example of such workflow.