Georeference Sentinel-3 Data with irregular grid

Question

I need to georeference and subset a lot of Sentinel-3-Syn-2-Files. Every data product provides seperate NetCDF-Files for each band and the actual (irregular) coordinate grid called 'geolocation'. So the basic structure is like in the question here. And like the suggested answer I tried it with:

gdal_translate -of XYZ NETCDF:"geolocation.nc":lat latitude.csv
gdal_translate -of XYZ NETCDF:"geolocation.nc":lon longitude.csv
gdal_translate -of XYZ NETCDF:"Syn_Oa17_reflectance.nc":SDR_Oa17 sdr_oa17.csv

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

I don't know exactly what the term {print($3/1000000,$6/1000000,$9)} does but it worked apparently. Despite being the question the suggested answer doesn't provide the step how to actually write the merged data back to e.g. NetCDF. I found something similar here. But when I do:

gdal_translate merged.csv combined.nc

I get:

ERROR 1: Ungridded dataset: At line 4, too many stepY values

There is a workaround for this kind of error. But the suggested step:

tail -n +2 merged.csv| sort -n -t ' ' -k2 -k1 > merged_sorted.xyz

Also provides a dataset which causes the same error:

gdalinfo merged_sorted.xyz

ERROR 1: Ungridded dataset: At line 2, X spacing was -0.038000. Expected >0 value gdalinfo failed - unable to open 'merged_sorted.xyz'.

For the sake of completeness, when I do gdal_translate I'm getting a warning message:

Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute

But as it produced the desired csv-files I ingored it. If anyone is interested in the dataset it can be found here. Does anyone know how to do this?

Answer 1

We found that the best way is by using GDAL VRTs. It is still a pain, but... You basically extract the lat and lon rasters from geolocation.nc to e.g. VRTs:

gdal.Translate("lat.vrt", 'NETCDF:"geolocation.nc":latitude', 
               format="VRT")
gdal.Translate("lon.vrt", 'NETCDF:"geolocation.nc":longitude',
               format="VRT")

Then, you create VRTs around the data you want (e.g. the reflectance in your case), using the lat and lon rasters as geolocation metadata. This is basically a text file, which we'll call temp.vrt (or whatever)

<VRTDataset rasterXSize="4865" rasterYSize="4091">
  <metadata domain="GEOLOCATION">
    <mdi key="X_DATASET">longitude.vrt</mdi>
    <mdi key="X_BAND">1</mdi>
    <mdi key="Y_DATASET">latitude.vrt</mdi>
    <mdi key="Y_BAND">1</mdi>
    <mdi key="PIXEL_OFFSET">0</mdi>
    <mdi key="LINE_OFFSET">0</mdi>
    <mdi key="PIXEL_STEP">1</mdi>
    <mdi key="LINE_STEP">1</mdi>
  </metadata>
  <VRTRasterBand dataType="Float32" band="1">
    <SimpleSource>
      <SourceFilename relativeToVRT="0">NETCDF:Oa01_radiance.nc</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="4865" RasterYSize="4091" DataType="Float32" BlockXSize="4865" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="4865" ySize="4091" />
      <DstRect xOff="0" yOff="0" xSize="4865" ySize="4091" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

You now have 3 extra VRT files. We just "warp" test.vrt and convert to a sensible format:

gdal.Warp("super_dooper.tif", "test.vrt", geoloc=True, 
          dstSRS="EPSG:4326", format="GTiff")

We have some code that warps the S3 data to match MODIS tiles that implements this approach. It's not pretty, but maybe it helps clarify the above.