Converting Sentinel-3 LST data to geotiff

[OGMalan]

My apologies for the long post, but I have been stuck attempting to convert Sentinel-3 LST netCDF data to geotiffs for a while now as part of my master's thesis, and my supervisors are getting rather impatient.

Firstly, an introduction to the data. It is packaged as follow:

The data is packaged as several .nc files, as well as a xml file containing metadata. The metadata suggests that the data should be in EPSG:4346:
<metadataObject ID="measurementFrameSet" classification="DESCRIPTION" category="DMD"> <metadataWrap mimeType="text/xml" vocabularyName="Sentinel-SAFE" textInfo="Frame Set"> <xmlData> <sentinel-safe:frameSet> <sentinel-safe:footPrint srsName="http://www.opengis.net/def/crs/EPSG/0/4326"> <gml:posList>-29.1992 8.81452 -29.1684 9.00931 -29.0806 9.52497 -28.9909 10.035 -28.896 10.5474 -28.801 11.0545 -28.7054 11.5734 -28.605 12.0778 -28.5082 12.5866 -28.4088 13.0983 -28.2994 13.6035 -28.197 14.1087 -28.0927 14.6147 -27.9761 15.1199 -27.8675 15.6236 -27.7564 16.1214 -27.6423 16.627 -27.525 17.1239 -27.4036 17.6145 -27.2806 18.1169 -27.1611 18.6177 -27.0344 19.1076 -26.9124 19.605 -26.781 20.0911 -26.6548 20.5894 -26.5286 21.0864 -26.3885 21.5669 -26.2565 22.0561 -26.1192 22.5398 -25.9879 23.0344 -25.8475 23.5168 -28.389 24.4405 -30.9637 25.4388 -33.5273 26.5043 -36.0695 27.6424 -36.078 27.6463 -36.2329 27.1108 -36.4066 26.5843 -36.5511 26.0457 -36.7063 25.5069 -36.8573 24.9628 -37.0105 24.4346 -37.1568 23.8819 -37.2966 23.3343 -37.4393 22.7945 -37.5805 22.24 -37.7181 21.6969 -37.8487 21.1329 -37.9835 20.584 -38.1138 20.0259 -38.2315 19.464 -38.3634 18.9072 -38.4821 18.3398 -38.5973 17.7783 -38.7109 17.2074 -38.8288 16.6385 -38.9366 16.0638 -39.0387 15.4947 -39.148 14.9205 -39.2443 14.3444 -39.3456 13.7696 -39.4421 13.1931 -39.5308 12.6058 -39.6274 12.0296 -39.7142 11.4424 -39.7416 11.2224 -39.7345 11.2088 -37.1002 10.5838 -34.4547 9.9742 -31.8069 9.37884 -29.1992 8.81452</gml:posList> </sentinel-safe:footPrint> </sentinel-safe:frameSet> </xmlData> </metadataWrap> </metadataObject>

Now I attempted to load the data into xarray using this:
import rioxarray import xarray from pyproj import CRS root = r'G:\MSc\Code\JupyterLab\temp\S3B_SL_2_LST____20201214T214826_20201214T215126_20201216T025309_0179_046_385_5400_LN2_O_NT_004.SEN3\*.nc' xds = xarray.open_mfdataset(root,combine = 'by_coords', decode_times=False)

This did not work, and output the following error:
ValueError: arguments without labels along dimension 'columns' cannot be aligned because they have different dimension sizes: {130, 1500}

Now I believe the problem with this is that the _in files have a 1200 x 1500 grid, whereas the _tx files have a 500 x 1500 grid. Working around this problem, we can open up only the _in files (or the _tx files, but I am using _in for the example):

If I am correct in interpreting the information, it appears that the netCDF does not have any spatial dimensions defined, despite that I have attempted to set the CRS using xds.rio.set_crs(4326) and exporting that to a geotiff using xds.to_netcdf(r'G:\MSc\TestImage\testout\testout.tiff'), but that did not work and seemed to only map the row and column values as coordinates. This did not work.

Now the Sentinel-3 LST product does have some form of spatial referencing in the data that I assume I can work with - the longitude_in/_tx and latitude_in/_tx bands. These are gridded values containing the geographic longitude/latitude coordinates for every cell in the netCDF. I have used this to extract point values from the netCDF grid in the past and have attempted to solve this problem by turning every cell into a point feature, but this proved far too slow to be of much use.

Anyway, I have a feeling that this information could be used somehow to georeference the netCDFs, either by creating a new set of spatial dimensions, or failing that they should be able to act as GCPs, but I fear that this second option may also be too slow. I have also tried using SNAP, but unfortunately it proved too slow and unstable to work for this problem. Speed is a bit of an issue, as I am aiming to process a year's worth of information over an AOI.

I would greatly appreciate any assistance on this matter, be it some code that can already do this, or even just guidance on what to attempt next. Thanks a lot in advance!

[snowman2]

it appears that the netCDF does not have any spatial dimensions defined

Are you able to open it with rioxarray.open_rasterio and see if it finds the spatial information?

I believe the problem with this is that the _in files have a 1200 x 1500 grid, whereas the _tx files have a 500 x 1500 grid.

You will need to make sure they are on the same grid with the same projection/resolution. rio.reproject_match will likely be useful at some point when you figure out the projection/resolution.

xds.rio.set_crs(4326) and exporting that to a geotiff using xds.to_netcdf(r'G:\MSc\TestImage\testout\testout.tiff')

I would recommend xds.rio.write_crs(4326, inplace=True) and xds.rio.to_raster(r'G:\MSc\TestImage\testout\testout.tiff')

[OGMalan]

Thank you so much for your help!

On to your suggestions:

Are you able to open it with rioxarray.open_rasterio and see if it finds the spatial information?
I have two problems with this approach. Firstly, if I understand it correctly, I will have to load the nc files one by one using that. Secondly, it returns this warning:
NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix be returned.

So I am assuming it does not find the spatial information.

You will need to make sure they are on the same grid with the same projection/resolution. rio.reproject_match will likely be useful at some point when you figure out the projection/resolution.

I am assuming that this would not be of much use until I can solve the CRS issue above, but it is good to know that I can use that.

I would recommend xds.rio.write_crs(4326, inplace=True) and xds.rio.to_raster(r'G:\MSc\TestImage\testout\testout.tiff')

My apologies, this was my fault. I was in between modifying my code to output a .nc. When I use the above code I get the following error:
TypeError: ufunc 'isnan' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Replacing the nan values leads to the following:
DimensionError: y dimension not found. 'rio.set_spatial_dims()' or using 'rename()' to change the dimension name to 'y' can address this.

I have managed to address this by renaming rows and columns to x and y respectively. I am aware that I have not added any CRS info by doing this, as the rox and column index are not spatially reference in any way, but even if I was able to turn them into a set of spatial dimensions I fear I would be left with the following upon calling xds.rio.to_raster(r'G:\MSc\TestImage\testout\testout.tiff'):
Only 2D and 3D data arrays supported. Data variable: stack-1beda8259771c2b4ad8770ea09ae0653

I am assuming that the above means that I must try to process each band individually? If this is the case, would the following workflow make sense:

1. Load the geodetic_in/_tx.nc file (containing the lat and lon coordinates for the _in and _tx grids).
2. Somehow georeference the geodetic file, using their coordinate info.
3. Iterate the remaining fiiles, adding this georeferencing info to them.
4. Output and save the individual bands as geotiffs, resample them to a common resolution and CRS and then stack them into a final, multiband output image.

On this note, following point 2 - would I be able to use the identity matrix returned by the rasterio operation with rasterio's rasterio.transform.from_gcps method?

[snowman2]

Only 2D and 3D data arrays supported.

.squeeze() may help.

would I be able to use the identity matrix returned by the rasterio operation with rasterio's rasterio.transform.from_gcps method?

You could write the transform with from_gcps

transform = from_gcps(...)
cds.rio.write_transform(transform, inplace=True)

EDIT: See #339

[remi-braun]

Hello,

I am continuing this discussion as it seems I have a similar need on the same satellite: I want to geocode any netcdf array from Sentinel-3
The moot point is the same than above: applying lon/lat arrays to a non-geocoded array, and this seems waaay less trivial than hoped.

What I am doing is (ie. on a Sentinel-3 OLCI product):

1. Load the netcdf array of latitude/longitude:

>>> # LATITUDE
>>> lat = rioxarray.open_rasterio('netcdf:S3A_OL_1_EFR____20191215T105023_20191215T105323_20191216T153115_0179_052_322_2160_LN1_O_NT_002.SEN3\\geo_coordinates.nc:latitude')
<xarray.DataArray 'latitude' (band: 1, y: 4091, x: 4865)>
array([[[41970420, 41970130, ..., 39543947, 39543248],
        [41972986, 41972697, ..., 39546495, 39545797],
        ...,
        [52450015, 52449789, ..., 49880462, 49879650],
        [52452576, 52452349, ..., 49882966, 49882154]]])
Coordinates:
  * band         (band) int32 1
  * x            (x) float64 0.5 1.5 2.5 3.5 ... 4.862e+03 4.864e+03 4.864e+03
  * y            (y) float64 0.5 1.5 2.5 3.5 ... 4.088e+03 4.09e+03 4.09e+03
    spatial_ref  int32 0
Attributes:
    long_name:      DEM corrected latitude
    scale_factor:   1e-06
    standard_name:  latitude
    units:          degrees_north
    valid_max:      90000000
    valid_min:      -90000000
    _FillValue:     -2147483648.0
    add_offset:     0.0

>>> # LONGITUDE
>>> lon = rioxarray.open_rasterio('netcdf:S3A_OL_1_EFR____20191215T105023_20191215T105323_20191216T153115_0179_052_322_2160_LN1_O_NT_002.SEN3\\geo_coordinates.nc:longitude')
<xarray.DataArray 'longitude' (band: 1, y: 4091, x: 4865)>
array([[[-17401957, -17398723, ...,  -2160966,  -2157957],
        [-17401548, -17398314, ...,  -2159975,  -2156965],
        ...,
        [-15779766, -15775813, ...,   2597806,   2601348],
        [-15779381, -15775428, ...,   2599187,   2602729]]])
Coordinates:
  * band         (band) int32 1
  * x            (x) float64 0.5 1.5 2.5 3.5 ... 4.862e+03 4.864e+03 4.864e+03
  * y            (y) float64 0.5 1.5 2.5 3.5 ... 4.088e+03 4.09e+03 4.09e+03
    spatial_ref  int32 0
Attributes:
    long_name:      DEM corrected longitude
    scale_factor:   1e-06
    standard_name:  longitude
    units:          degrees_east
    valid_max:      180000000
    valid_min:      -180000000
    _FillValue:     -2147483648.0
    add_offset:     0.0

---

2. Compute the transform from those two arrays:

# Compute the transform
>>> tr = rasterio.transform.from_bounds(
    west=np.min(lon.data)*lon.scale_factor, 
    south=np.min(lat.data) * lat.scale_factor, 
    east=np.max(lon.data) * lon.scale_factor, 
    north=np.max(lat.data) * lat.scale_factor, 
    width=lat.x.size, 
    height=lat.y.size)

⚠️ However here I am not sure I am not losing a lot of information, does a simple GDAL transform suffice to describe the coordinates of a big Sentinel-3 product ?

---

3. Open a radiance band

>>> path = r"netcdf:S3A_OL_1_EFR____20191215T105023_20191215T105323_20191216T153115_0179_052_322_2160_LN1_O_NT_002.SEN3\Oa06_radiance.nc:Oa06_radiance"
>>> A = rioxarray.open_rasterio(path)
<xarray.DataArray 'Oa06_radiance' (band: 1, y: 4091, x: 4865)>
[19902715 values with dtype=uint16]
Coordinates:
  * band         (band) int32 1
  * x            (x) float64 0.5 1.5 2.5 3.5 ... 4.862e+03 4.864e+03 4.864e+03
  * y            (y) float64 0.5 1.5 2.5 3.5 ... 4.088e+03 4.09e+03 4.09e+03
    spatial_ref  int32 0
Attributes:
    add_offset:           0.0
    ancillary_variables:  Oa06_radiance_err
    coordinates:          time_stamp altitude latitude longitude
    long_name:            TOA radiance for OLCI acquisition band Oa06
    scale_factor:         0.012353800237178802
    standard_name:        toa_upwelling_spectral_radiance
    units:                mW.m-2.sr-1.nm-1
    valid_max:            65534
    valid_min:            0
    _FillValue:           65535.0

---

4. Apply CRS and transform to the opened band

>>> A.rio.write_crs("EPSG:4326", inplace=True)
<xarray.DataArray 'Oa06_radiance' (band: 1, y: 4091, x: 4865)>
[19902715 values with dtype=uint16]
Coordinates:
  * band         (band) int32 1
  * x            (x) float64 0.5 1.5 2.5 3.5 ... 4.862e+03 4.864e+03 4.864e+03
  * y            (y) float64 0.5 1.5 2.5 3.5 ... 4.088e+03 4.09e+03 4.09e+03
    spatial_ref  int32 0
Attributes:
    add_offset:           0.0
    ancillary_variables:  Oa06_radiance_err
    coordinates:          time_stamp altitude latitude longitude
    long_name:            TOA radiance for OLCI acquisition band Oa06
    scale_factor:         0.012353800237178802
    standard_name:        toa_upwelling_spectral_radiance
    units:                mW.m-2.sr-1.nm-1
    valid_max:            65534
    valid_min:            0
    _FillValue:           65535.0

>>> A.rio.write_transform(transform=tr, inplace=True)
<xarray.DataArray 'Oa06_radiance' (band: 1, y: 4091, x: 4865)>
[19902715 values with dtype=uint16]
Coordinates:
  * band         (band) int32 1
  * x            (x) float64 0.5 1.5 2.5 3.5 ... 4.862e+03 4.864e+03 4.864e+03
  * y            (y) float64 0.5 1.5 2.5 3.5 ... 4.088e+03 4.09e+03 4.09e+03
    spatial_ref  int32 0
Attributes:
    add_offset:           0.0
    ancillary_variables:  Oa06_radiance_err
    coordinates:          time_stamp altitude latitude longitude
    long_name:            TOA radiance for OLCI acquisition band Oa06
    scale_factor:         0.012353800237178802
    standard_name:        toa_upwelling_spectral_radiance
    units:                mW.m-2.sr-1.nm-1
    valid_max:            65534
    valid_min:            0
    _FillValue:           65535.0

⚠️ As you see here, the coordinates are not recomputed, and I don't know how to do that
⚠️ Moreover, only the cached transform is set, when the transform is recalculated, it goes back to the identity matrix

---

5. Write the raster

>>> A.rio.to_raster(r"test_green.tif", recalc_transform=False)

⚠️ As you can see, I HAVE TO set recalc_transform=False, as the recompted transform sets back the identity

---

Can you help me on those warnings ⚠️ ?

Context:

I am currently trying to remove SNAP from my lib (EOReader).
I have to main axes:

• remove from SAR processing (that will be difficult)
• remove from Sentinel-3 processing: this is about what I am asking here

[remi-braun]

Creating a transform with gcps seems to work better, but is still not exact. Maybe a reprojection is needed ? I am a bit lost here

Here is my code for creating the gcps and the corresponding transform:

from rasterio.control import GroundControlPoint
nof_gcp_x = np.arange(0, A.x.size, 100)
nof_gcp_y = np.arange(0, A.y.size, 100)
alt = rioxarray.open_rasterio(r"netcdf:S3A_OL_1_EFR____20191215T105023_20191215T105323_20191216T153115_0179_052_322_2160_LN1_O_NT_002.SEN3\geo_coordinates.nc:altitude")
gcps = []
id = 0
for x in nof_gcp_x:
    for y in nof_gcp_y:
        gcps.append(GroundControlPoint(
                                 row=y, 
                                 col=x, 
                                 x=lon.data[:, y, x]*lon.scale_factor, 
                                 y=lat.data[:, y, x] * lat.scale_factor, 
                                 z=alt.data[:, y, x] * alt.scale_factor, 
                                 id=id))
        id += 1
        
tr_gcp = rasterio.transform.from_gcps(gcps)

The obtained transform is a bit better (we got the rotation) but still cannot describe the raster:

Affine(0.0034268101089098874, 0.0007897957564530836, -18.0695103078887,
       -0.0005084449740915641, 0.0025636394854227575, 42.165547283051865)

See hereunder the problem:

[snowman2]

See: #209. It looks like you may have an irregular grid.

[snowman2]

Also see #339

[remi-braun]

I was trying to reproject my data with the GCPs, but I think I'm doing it wrong:

>>> A.rio.set_crs("EPSG:4326", inplace=True)
>>> A.rio.reproject(dst_crs="EPSG:4326", gcps=gcps)

The error I get is:

Traceback (most recent call last):
  File "<input>", line 2, in <module>
  File "C:\Users\rbraun\Anaconda3\envs\eoreader\lib\site-packages\rioxarray\raster_array.py", line 399, in reproject
    self._obj, self.crs, dst_crs, resolution, shape, **kwargs
  File "C:\Users\rbraun\Anaconda3\envs\eoreader\lib\site-packages\rioxarray\raster_array.py", line 134, in _make_dst_affine
    **kwargs,
  File "C:\Users\rbraun\Anaconda3\envs\eoreader\lib\site-packages\rasterio\env.py", line 387, in wrapper
    return f(*args, **kwds)
  File "C:\Users\rbraun\Anaconda3\envs\eoreader\lib\site-packages\rasterio\warp.py", line 509, in calculate_default_transform
    src_crs, dst_crs, width, height, left, bottom, right, top, gcps, rpcs, **kwargs
  File "rasterio\_warp.pyx", line 697, in rasterio._warp._calculate_default_transform
  File "rasterio\_warp.pyx", line 676, in rasterio._warp._calculate_default_transform
  File "rasterio\_err.pyx", line 215, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_AppDefinedError: Failed to compute GCP transform: Transform is not solvable

However, as shown hereover, I can compute a transform from my GCPs, so I am confused 😰

[remi-braun]

OK the problem was me creating my GCPs all along:
The corrected version is:

gcps = []
id = 0
for x in nof_gcp_x:
    for y in nof_gcp_y:
        
        gcps.append(GroundControlPoint(row=y, col=x, x=LON.data[0, y, x]*LON.scale_factor, y=LAT.data[0, y, x] * LAT.scale_factor, z=ALT.data[0, y, x] * ALT.scale_factor, id=id))
        id += 1

Aaaaaaand this almost works... However I got some errors in the corners, independant of the number of GCPs.

This is solved by reprojecting with thin plate spline option (as you specified):

Ap = A.rio.reproject(dst_crs="EPSG:4326",
                gcps=gcps, **{"SRC_METHOD": "GCP_TPS"})

[HafezAhmad]

any one tried with OLCI level 2 data? I am tired of searching and how to convert to geotif. Thanks