NetCDF with separate lat/lon bands to GeoTiff with Python/GDAL

Question

I know this has been asked before, but none of the suggested solutions works for me, and Converting Sentinel-3 data (netcdf) to GeoTIFF is for Java. I tried the accepted answer from Map project a raster having separate latitude and longitude raster bands, but got the error message

ERROR 1: Unable to compute a transformation between pixel/line
    and georeferenced coordinates for Oa01_radiance.tif.
    There is no affine transformation and no GCPs.

although I made sure that my latitude.vrt and longitude.vrt have the same dimensions as Oa01_radiance.vrt and also point to the correct files.

---

So I have some Sentinel-3 data, which come as a folder of a number of NetCDF's, one per radiance band, together with another NetCDF holding tie points for every x-th pixel (every 64th, in that case). I managed to interpolate the lats and lons and create single tifs with the same dimensions as the radiance bands, so basically the lat/lon information for each pixel (I checked them, they are correct). From the other questions around the web, I tried this:

# create vrt for data, lat and lon
gdal_translate -of VRT lat.tif latitude.vrt
gdal_translate -of VRT lon.tif longitude.vrt
gdal_translate -of VRT NETCDF:"Oa01_radiance.nc" Oa01_radiance.vrt

Then I edited data.vrt to only contain the following:

<VRTDataset rasterXSize="4865" rasterYSize="4091">
  <metadata domain="GEOLOCATION">
    <mdi key="X_DATASET">MY_FILE_PATH\longitude.vrt</mdi>
    <mdi key="X_BAND">1</mdi>
    <mdi key="Y_DATASET">MY_FILE_PATH\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>

and then tried

gdalwarp -geoloc -t_srs epsg:4326 Oa01_radiance.vrt Oa01_radiance.tif

which resulted in the above mentioned error. So

1. is there something I'm missing and
2. is there some Python code that can convert the data from the given structure to GeoTiff?

---

One more comment: using min/max of the lat/lon bands (even the interpolated ones) and calculating the GeoTransformObject from them results in unacceptable offsets and distortions of several hundred km.

Answer 1

Here is a solution that works for the Sentinel 3 Land Surface Temperature (LST) product, and it should work for the other products as well (you might need to account for different pixel steps in the lat/lon grid, though). It has been implemented in R, but it should work in any other language that can interface to GDAL.

  library("rgdal")
  library("gdalUtils")
  library("raster")

  #Sen 3 LST data have lat/lon coordinates stored in seperate nc files:
  geolayer <- "geodetic_in.nc"
  #while this nc contains our target parameter (LST):
  lstlayer <- "LST_in.nc"
 #You need to build a VRT from lat, lon and target band nc's:

  lat_nc.vrt <- "LAT.vrt"
  lon_nc.vrt <- "LON.vrt"
  lst_nc.vrt <- "LST.vrt"

  gdalbuildvrt(gdalfile=geolayer, sd = 3, output.vrt=lat_nc.vrt, separate=TRUE,verbose=TRUE)
  gdalbuildvrt(gdalfile=geolayer, sd = 5, output.vrt=lon_nc.vrt, separate=TRUE,verbose=TRUE)
  gdalbuildvrt(gdalfile=lstlayer, sd = 1, output.vrt=lst_nc.vrt, separate=TRUE,verbose=TRUE)

You need to edit the VRT XML in order to correctly treat scale and offset of the nc's (which GDAL does not, neither in building the VRT nor in gdalwarp). If you do not
treat this, your data will be both wrongly scaled and wrongly located. To do the scaling correctly, it is necessary to output to temporary files, e.g. as a tif:

  writeRaster(raster(lat_nc.vrt), "lat_in.tif", overwrite = T)
  writeRaster(raster(lon_nc.vrt), "lon_in.tif", overwrite = T)
  writeRaster(raster(lst_nc.vrt), "lst_in.tif", overwrite = T)

Having done so, you need to build a VRT for the output temporary tif of your target parameter (in this case LST) and again edit the VRT tif file.

  gdalbuildvrt(gdalfile="lst_in.tif", sd = 1, output.vrt="LST_tif.vrt", separate=TRUE,verbose=TRUE)


  res <- gdalwarp(srcfile = "LST_tif.vrt", geoloc = TRUE, t_srs = "EPSG:4326", dstfile = outTifName)

The editing of the VRT files can be done using XML functions of your programming language of choice (I used the R XML package). It needs to follow the guidelines outlined here: Unable to warp HDF5 files

As said, before converting the files to tif you need to correctly treat scaling and offset (this is found in the nc metadata for each band) - here is an example for the latitude band:

<VRTDataset rasterXSize="1500" rasterYSize="1200">
  <VRTRasterBand dataType="Int32" band="1">
    <NoDataValue>-2147483648</NoDataValue>
    <Scale>1e-06</Scale>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">NETCDF:"D:/Data/Sentinel3/in/S3A_SL_2_LST____20190205T093623_20190205T093923_20190205T112255_0179_041_093_2880_LN2_O_NR_003.SEN3/geodetic_in.nc":latitude_in</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1500" RasterYSize="1200" DataType="Int32" BlockXSize="1500" BlockYSize="1"/>
      <SrcRect xOff="0" yOff="0" xSize="1500" ySize="1200"/>
      <DstRect xOff="0" yOff="0" xSize="1500" ySize="1200"/>
      <NODATA>-2147483648</NODATA>
    </ComplexSource>
    <Offset>0.0</Offset>
    <Scale>1.0E-6</Scale>
  </VRTRasterBand>
</VRTDataset>

For the editing of the LST VRT to be georeferenced, you need to add the metadata section for gdalwarp:

<VRTDataset rasterXSize="1500" rasterYSize="1200">
  <GeoTransform>  0.0000000000000000e+00,  1.0000000000000000e+00,  0.0000000000000000e+00,  1.2000000000000000e+03,  0.0000000000000000e+00, -1.0000000000000000e+00</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>-3.4e+38</NoDataValue>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">lst_in.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1500" RasterYSize="1200" DataType="Float32" BlockXSize="1500" BlockYSize="1"/>
      <SrcRect xOff="0" yOff="0" xSize="1500" ySize="1200"/>
      <DstRect xOff="0" yOff="0" xSize="1500" ySize="1200"/>
      <NODATA>-3.4e+38</NODATA>
    </ComplexSource>
  </VRTRasterBand>
  <metadata domain="GEOLOCATION">
    <mdi key="X_DATASET">D:/Data/Sentinel3/temp/lon_in.tif</mdi>
    <mdi key="X_BAND">1</mdi>
    <mdi key="Y_DATASET">D:/Data/Sentinel3/temp/lat_in.tif</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>
</VRTDataset>

Answer 2

I added this pre-processing data from netCDF to tiff with gdal binary, in case, someone needs:

    Example OLCI level 2 data: S3A_OL_2_LFR____20180131T225040_20180131T225340_20180202T040253_0180_027_215_2520_LN1_O_NT_002.SEN3

+ First, create the lat.vrt file:

gdal_translate -of VRT NETCDF:"/home/rasdaman/RASDAMAN_PROJECTS/PARSEC_MUNDI/Sentinel-3/S3A_OL_2_LFR____20180131T225040_20180131T225340_20180202T040253_0180_027_215_2520_LN1_O_NT_002.SEN3/geo_coordinates.nc":latitude lat.vrt

+ Then, create the lon.vrt file:

gdal_translate -of VRT NETCDF:"/home/rasdaman/RASDAMAN_PROJECTS/PARSEC_MUNDI/Sentinel-3/S3A_OL_2_LFR____20180131T225040_20180131T225340_20180202T040253_0180_027_215_2520_LN1_O_NT_002.SEN3/geo_coordinates.nc":longitude lon.vrt

+ Then, create the band.vrt file manually with this XML below, for example: iwv.nc (IWV):
change sizes of X and Y and paths to lat.vrt and lon.vrt

<VRTDataset rasterXSize="4865" rasterYSize="4090">
 <metadata domain="GEOLOCATION">
 <mdi key="X_DATASET">/home/rasdaman/RASDAMAN_PROJECTS/PARSEC_MUNDI/Sentinel-3/S3A_OL_2_LFR____20180131T225040_20180131T225340_20180202T040253_0180_027_215_2520_LN1_O_NT_002.SEN3/lon.vrt</mdi>
 <mdi key="X_BAND">1</mdi>
 <mdi key="Y_DATASET">/home/rasdaman/RASDAMAN_PROJECTS/PARSEC_MUNDI/Sentinel-3/S3A_OL_2_LFR____20180131T225040_20180131T225340_20180202T040253_0180_027_215_2520_LN1_O_NT_002.SEN3/lat.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="Byte" band="1">

    <SimpleSource>
      <SourceFilename relativeToVRT="1">NETCDF:"iwv.nc":IWV</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="4865" RasterYSize="4090" dataType="Byte"/>
      <SrcRect xOff="0" yOff="0" xSize="4865" ySize="4090" />
      <DstRect xOff="0" yOff="0" xSize="4865" ySize="4090" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

+ Then, run this gdalwarp command to project netCDF to tiff.
 export GDAL_NETCDF_BOTTOMUP=NO
 gdalwarp -geoloc -t_srs EPSG:4326 band.vrt test.tif -srcnodata 255 -overwrite

Answer 3

Thanks to GRuC, I was able to convert the data. It was a bit of a pain, but I translated GRuC's suggestion into a Python script, which I posted into my GitHub because of its length. There, I made use of the GDAL installation, calling all necessary modules via subprocess. Feel free to comment, use and improve (I'm sure this can be written much cleaner).