4

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.

6
  • Which GDAL version do you have ? I have no error messages with GDAL 2.2.2, released 2017/09/15. But it is painfully slow because GDAL seems to have problems reading the compressed netCDF. I had good results when I first copied the files with nccopy -d 0 to remove compression. With the results I'm using more or less the same solution as you do. I first converted the Band to GeoTiff for better performance and did not build VRTs of lat.tif and lon.tif, just linked the VRT to these.
    – pLumo
    Mar 5, 2018 at 9:52
  • @RoVo - To be honest, I'm still running on GDAL version 1.11.3, released 2015/09/16, and currently have no possibility of updating.
    – s6hebern
    Mar 5, 2018 at 9:59
  • Ok, that might be problem. Although it's far from perfect, there has been a lot of improvements in netCDF handling.
    – pLumo
    Mar 5, 2018 at 10:19
  • So you're basically saying that I have no chance of correctly georeferencing my image with Python/GDAL unless I can update my GDAL version?
    – s6hebern
    Mar 5, 2018 at 10:37
  • I don't say that, but I only know it works with 2.2.2
    – pLumo
    Mar 5, 2018 at 10:39

3 Answers 3

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>

1
  • This is not python solution. Apr 3, 2020 at 8:23
1

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

0
0

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).

2
  • Hi s6hebern where is the code which you have tried.
    – RRSC NGP
    May 15, 2022 at 2:11
  • @RRSC NGP, I updated the link in the answer, so it should be working again. But this script is already more than 3 years old and I have not tried it out since, so it might require some fixes. Additionally, it was still written in Py2. Maybe it can help you anyways.
    – s6hebern
    May 17, 2022 at 10:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.