Georeferencing Sentinel-3 LST product

Short answer:

1. use GDAL tools to extract lat and long data from an .nc file.
2. Merge it with the data from any band you want to see to make an x,y,z dataset.
3. Convert the irregular x,y,z data to a shapefile and apply a CRS.
4. Grid the shapefile at desired resolution.

The lat/long data are not regularly sampled, so the data can't be loaded as an XYZ raster directly, and need to be gridded.

This flow works on the one dataset that I tried.

   gdal_translate -of XYZ NETCDF:"S4_radiance_an.nc":S4_radiance_an S4_radiance.csv
   gdal_translate -of XYZ NETCDF:"geodetic_an.nc":latitude_an latitude.csv
   gdal_translate -of XYZ NETCDF:"geodetic_an.nc":longitude_an longitude.csv

   echo X,Y,Z >merged.csv
   paste longitude.csv latitude.csv S4_radiance.csv|gawk 'BEGIN{OFS=","}{print($3/1000000,$6/1000000,$9)}' >>merged.csv
   ogr2ogr -a_srs "EPSG:4326"  -mapFieldType All=Integer -oo X_POSSIBLE_NAMES=X* -oo Y_POSSIBLE_NAMES=Y* -oo KEEP_GEOM_COLUMNS=NO newS4_radiance.shp merged.csv

Long Answer

Use GDAL tools from the command line. For example, file S4_radiance_an.nc

gdalinfo S4_radiance_an.nc

At the bottom of the output, there are lines like:

Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"S4_radiance_an.nc":S4_exception_an
  SUBDATASET_1_DESC=[2400x3000] S4_exception_an (8-bit unsigned integer)
  SUBDATASET_2_NAME=NETCDF:"S4_radiance_an.nc":S4_exception_orphan_an
  SUBDATASET_2_DESC=[2400x374] toa_radiance_status_flag (8-bit unsigned integer)
  SUBDATASET_3_NAME=NETCDF:"S4_radiance_an.nc":S4_radiance_an
  SUBDATASET_3_DESC=[2400x3000] toa_radiance (16-bit integer)
  SUBDATASET_4_NAME=NETCDF:"S4_radiance_an.nc":S4_radiance_orphan_an
  SUBDATASET_4_DESC=[2400x374] toa_radiance (16-bit integer)

So resolution is 2400 x 3000 for bands that contain 'an' (nadir c.f. oblique).

The same information from geodetic_an.nc yields:

Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"geodetic_an.nc":elevation_an
  SUBDATASET_1_DESC=[2400x3000] surface_altitude (16-bit integer)
  SUBDATASET_2_NAME=NETCDF:"geodetic_an.nc":elevation_orphan_an
  SUBDATASET_2_DESC=[2400x374] surface_altitude (16-bit integer)
  SUBDATASET_3_NAME=NETCDF:"geodetic_an.nc":latitude_an
  SUBDATASET_3_DESC=[2400x3000] latitude (32-bit integer)
  SUBDATASET_4_NAME=NETCDF:"geodetic_an.nc":latitude_orphan_an
  SUBDATASET_4_DESC=[2400x374] latitude (32-bit integer)
  SUBDATASET_5_NAME=NETCDF:"geodetic_an.nc":longitude_an
  SUBDATASET_5_DESC=[2400x3000] longitude (32-bit integer)
  SUBDATASET_6_NAME=NETCDF:"geodetic_an.nc":longitude_orphan_an
  SUBDATASET_6_DESC=[2400x374] longitude (32-bit integer)

Resolution is again 2400 x 3000, so we assume they are related.

The NAME and DESCRIPTION of set 3 are:

  SUBDATASET_3_NAME=NETCDF:"geodetic_an.nc":latitude_an
  SUBDATASET_3_DESC=[2400x3000] latitude (32-bit integer)

We can use the NAME to load using GDAL tools. First, translate to XYZ format with:

gdal_translate -of XYZ NETCDF:"geodetic_an.nc":latitude_an latitude.csv

Similarly for the longitude file. In each case, records are like:

2990.5 2399.5 7575671
2991.5 2399.5 7574636
2992.5 2399.5 7573601

So the format is I, J, Latitude, and latitude is in (degrees * 1000000).

We can merge the XYZ files representing Latitude, Longitude and, say, Radiance (from S4_radiance_an.nc, for example, which is also 2400 x 3000 in the band "S4_radiance_an.nc":S4_radiance_an"). The first and last records of these three files are:

S4_radiance.csv
0.5 0.5 -32768
1.5 0.5 -32768
2.5 0.5 -32768
2997.5 2399.5 -32768
2998.5 2399.5 -32768
2999.5 2399.5 -32768
latitude.csv
0.5 0.5 -73571
1.5 0.5 -74567
2.5 0.5 -75563
2997.5 2399.5 7568425
2998.5 2399.5 7567390
2999.5 2399.5 7566355
longitude.csv
0.5 0.5 95229130
1.5 0.5 95233511
2.5 0.5 95237892
2997.5 2399.5 110768320
2998.5 2399.5 110772730
2999.5 2399.5 110777141

So they have the same limits and are therefore related. Merge them using any method you like. They have the same number of rows (wc <S4_radiance.csv =>> 7200000 21600000 127768857), so on Unix paste and AWK will work:

paste longitude.csv latitude.csv S4_radiance.csv|gawk '{print($3/1000000,$6/1000000,$9)}' >merged.csv

This gives a file like:

X,Y,Z
95.2291,-0.073571,-32768
95.2335,-0.074567,-32768
.
.
110.773,7.56739,-32768
110.777,7.56635,-32768

The Z values are -32768, which means nodata, and is normal at the edge of a raster. X and Y values are in degrees These can be loaded as point data and gridded as required.