projecting INSAT3D L1B data from HDF5 to Geotiff

Question

I am attempting to use gdal to project a set of geostationary INSAT 3D L1B images.

I tried to follow the steps in Unable to warp HDF5 files

1. I created the .vrt files using

   gdal_translate -of VRT HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://Latitude lat.vrt 
   gdal_translate -of VRT HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://Longitude lon.vrt    
   gdal_translate -of VRT HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://IMG_MIR mir.vrt

2. I edited the .vrt files to remove the gcp lists and add the geolocation section to the mir.vrt file

3. I attempted to project the .vrt using

gdalwarp -geoloc -t_srs EPSG:4326 mir.vrt mir-out.tif

The resulting image is not what I expect

I re checked my input .vrt files for errors. I viewed the mir.vrt in QGIS and it does not seem to be correct

Could you advise me on how to correctly re-project this data.

Update 31 Aug 1700 HRS IST:

I looked further and thought that my problem could be sorted out by specifying the correct extents for my data. I searched and found, which I tried to follow

1) Georeferencing Himawari-8 in GDAL (or other)

2) Transforming geostationary satellite image to lon/lat

I used gdalinfo to find that the Observed_Altitude(km)=35782.063169 The INSAT TIR band has a spatial resolution of 4km

GDALINFO reports a size of 2816 x 2805 for the IMG_MIR dataset

So our extents from the centre are 2816 / 2 * 4000 = 5632000 and 2805 / 2 * 4000 = 5610000 in x and y respectively

I then ran

   gdal_translate -a_srs  "+proj=geos +a=6378137.0 +b=6356752.3 +lon_0=82 +h=35782063  +no_defs" -a_ullr -5632000 5610000 5632000 -5610000 HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://IMG_MIR temp.tif

   gdalwarp -overwrite -t_srs "+proj=latlong +ellps=WGS84 +pm=82" -te 0.8432 -81.041527 163.15671 81.041527 temp.tif insat.tif

   gdal_translate -projwin 68 38 98 6 insat.tif indinsat.tif

This worked and I could overload the India boundary successfully (it could be a slightly better fit esp on the west)

I reasoned that using the geolocation arrays should give me a better correction. Therefore I tried to create the .vrt files with what I think is the correct height and extents

   gdal_translate -of VRT -a_srs  "+proj=geos +a=6378137.0 +b=6356752.3 +lon_0=82 +h=35782063  +no_defs" -a_ullr -5632000 5610000 5632000 -5610000 HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://IMG_MIR mir.vrt

   gdal_translate -of VRT -a_srs  "+proj=geos +a=6378137.0 +b=6356752.3 +lon_0=82 +h=35782063  +no_defs" -a_ullr -5632000 5610000 5632000 -5610000 HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://Latitude lat.vrt

   gdal_translate -of VRT -a_srs  "+proj=geos +a=6378137.0 +b=6356752.3 +lon_0=82 +h=35782063  +no_defs" -a_ullr -5632000 5610000 5632000 -5610000 HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://Longitude lon.vrt

The mir.vrt now displays correctly in QGIS

And I ran

   gdalwarp -geoloc -overwrite -t_srs "+proj=latlong +ellps=WGS84 +pm=82" -te 0.8432 -81.041527 163.15671 81.041527 mir.vrt insat-mir.tif

for the final step. The result however is still wrong.

Update 01 Sep 1300 HRS IST:

The lat and long arrays are scaled.

The .vrt files show

   <Metadata>
  <MDI key="Latitude_add_offset">0 </MDI>
  <MDI key="Latitude_long_name">latitude</MDI>
  <MDI key="Latitude_scale_factor">0.0099999998 </MDI>
  <MDI key="Latitude_units">degrees_north</MDI>
  <MDI key="Latitude__FillValue">32767 </MDI>
</Metadata>

I attempted to correct this by passing the -unscale keyword to gdal_translate. I expected gdal_translate to use the offset and scale factor from the metadata shown above. This however did not work out

   gdallocationinfo -b 1 lat.vrt 2000 2000

showed the unscaled values.

I forced the unscaling by using the scale keyword as follows. (Is there a more elegant method ?)

   gdal_translate -of VRT -a_srs  "+proj=geos +a=6378137.0 +b=6356752.3 +lon_0=82 +h=35782063  +no_defs" -a_ullr -5632000 5610000 5632000 -5610000 HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://IMG_MIR mir.vrt

   gdal_translate -of VRT -scale -9000 9000 -90 90 -a_srs  "+proj=geos +a=6378137.0 +b=6356752.3 +lon_0=82 +h=35782063  +no_defs" -a_ullr -5632000 5610000 5632000 -5610000 HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://Latitude lat.vrt

   gdal_translate -of VRT -scale -18000 18000 -180 180 -a_srs "+proj=geos +a=6378137.0 +b=6356752.3 +lon_0=82 +h=35782063  +no_defs" -a_ullr -5632000 5610000 5632000 -5610000 HDF5:"3DIMG_27APR2016_0830_L1B_STD.h5"://Longitude lon.vrt

I ran (after editing the mrt.vrt to add the geolocation section (as done previously, I did not delete the gcps, since I was forcing the use of the geolocation arrays)

   gdalwarp -geoloc -overwrite -t_srs "+proj=latlong +ellps=WGS84 +pm=82" -te 0.8432 -81.041527 163.15671 81.041527 mir.vrt insat-mir1.tif

to obtain the correct result. I could overlay the India boundary correctly (the problem noted in the case without the use of geolocation arrays was corrected)

I note that not passing the extents results in an incorrect result

Answer 1

The geolocation arrays for the level 1b INSAT 3D product are scaled. (scale factor 0.0099999998)

The creation of the .vrt files should therefore use the -unscale keyword with gdal_translate. That solution however did not work.

I forced scaling for the latitude and longitude geolocation arrays using the scale keyword and passing ranges appropriately ( -9000 9000 -90 90 for latitude) in the creation of .vrt using gdal_translate.

Details are included in the second update to my question.