2

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
  1. I edited the .vrt files to remove the gcp lists and add the geolocation section to the mir.vrt file

  2. 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
Result_of_geolocation

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

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 Correctly_done (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

  • I assume the -te is in EPSG:4326? – AndreJ Aug 31 '17 at 12:09
  • Thanks. I intended it to be in EPSG:4326. Your comment got me thinking why I should pass extents when I am using the geoloc arrays. I figured that out and have edited my question. – user104715 Sep 1 '17 at 7:09
  • Trying your last update, I did delete the gcp point section, which sounds reasonable if you want to use the arrays instead. – AndreJ Jan 16 '18 at 14:34
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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