Georeference calculation for msg datasets is incorrect

[jdroenner]

I'm working with MSG data and i noticed that while VIS_IR and HRV have the same corner coordinates the HRV channel appears to be shifted. When comparing VIS_IR or HRV with other data e.g. SRTM both appear to be shifted.

I think the problem is the calculation of the georeference information in msgdataset.cpp ~line 217. The existing code assumes that the MSG data is symmetric when calculating rMinX and rMaxY:

double rPixelSizeX;
double rPixelSizeY;
double rMinX;
double rMaxY;

if (command.channel[11] != 0)
{
  rPixelSizeX = 1000 * pp.idr()->ReferenceGridHRV->ColumnDirGridStep;
  rPixelSizeY = 1000 * pp.idr()->ReferenceGridHRV->LineDirGridStep;
  rMinX = -rPixelSizeX * (pp.idr()->ReferenceGridHRV->NumberOfColumns / 2.0); // assumption: (0,0) falls in centre
  rMaxY = rPixelSizeY * (pp.idr()->ReferenceGridHRV->NumberOfLines / 2.0);
}
else
{
  rPixelSizeX = 1000 * pp.idr()->ReferenceGridVIS_IR->ColumnDirGridStep;
  rPixelSizeY = 1000 * pp.idr()->ReferenceGridVIS_IR->LineDirGridStep;
  rMinX = -rPixelSizeX * (pp.idr()->ReferenceGridVIS_IR->NumberOfColumns / 2.0); // assumption: (0,0) falls in centre
  rMaxY = rPixelSizeY * (pp.idr()->ReferenceGridVIS_IR->NumberOfLines / 2.0);
}

The ​MSG Level 1.5 Image Data Format Description defines the coordinate system used for MSG data on page 22f. For the reference grid of the VIS_IR channels (3712x3712 pixels) the center or sub satellite point is defined as the center of the pixel at (1856,1856) where the origin is the SE corner. For the reference grid of the HRV channel (11136x11136 pixels) the center is defined as the center of the pixel at (5566,5566) also using the SE corner as origin.

A simple fix would be to introduce the central pixels as constants:

const int MSGDataset::iCenterPixelVIS_IR = 1856;
const int MSGDataset::iCenterPixelHRV = 5566;

And use the following code to calculate rMinX and rMaxY:

    if (command.channel[11] != 0)
    {
      rPixelSizeX = 1000 * pp.idr()->ReferenceGridHRV->ColumnDirGridStep;
      rPixelSizeY = 1000 * pp.idr()->ReferenceGridHRV->LineDirGridStep;
      rMinX = -rPixelSizeX * (pp.idr()->ReferenceGridHRV->NumberOfColumns - iCenterPixelHRV + 0.5); // The MSG Level 1.5 Image Data Format Description defines the center for HRV as the center of pixel 5566 from SE
      rMaxY = rPixelSizeY * (pp.idr()->ReferenceGridHRV->NumberOfLines - iCenterPixelHRV + 0.5);
    }
    else
    {
      rPixelSizeX = 1000 * pp.idr()->ReferenceGridVIS_IR->ColumnDirGridStep;
      rPixelSizeY = 1000 * pp.idr()->ReferenceGridVIS_IR->LineDirGridStep;
      rMinX = -rPixelSizeX * (pp.idr()->ReferenceGridVIS_IR->NumberOfColumns - iCenterPixelVIS_IR + 0.5); // The MSG Level 1.5 Image Data Format Description defines the center as the center of pixel 1856 from SE
      rMaxY = rPixelSizeY * (pp.idr()->ReferenceGridVIS_IR->NumberOfLines - iCenterPixelVIS_IR + 0.5);
    }

This solution works for the operational data from MSG. Gdalinfo using the adapted code produces:

VIS_IR
Origin = (-5570248.477339744567871,5570248.477339744567871)
Pixel Size = (3000.403165817260742,-3000.403165817260742)
Corner Coordinates:
Upper Left  (-5570248.477, 5570248.477) 
Lower Left  (-5570248.477,-5567248.074) 
Upper Right ( 5567248.074, 5570248.477) 
Lower Right ( 5567248.074,-5567248.074) 
Center      (   -1500.202,    1500.202) 

HRV
Origin = (-5571248.390376567840576,5571248.390376567840576)
Pixel Size = (1000.134348869323730,-1000.134348869323730)
Corner Coordinates:
Upper Left  (-5571248.390, 5571248.390) 
Lower Left  (-5571248.390,-5566247.719) 
Upper Right ( 5566247.719, 5571248.390) 
Lower Right ( 5566247.719,-5566247.719) 
Center      (   -2500.336,    2500.336) 

Gdalinfo using the old code produces:

VIS_IR
Origin = (-5568748.275756835937500,5568748.275756835937500)
Pixel Size = (3000.403165817260742,-3000.403165817260742)
Corner Coordinates:
Upper Left  (-5568748.276, 5568748.276)
Lower Left  (-5568748.276,-5568748.276)
Upper Right ( 5568748.276, 5568748.276)
Lower Right ( 5568748.276,-5568748.276)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)

HRV
Origin = (-5568748.054504394531250,5568748.054504394531250)
Pixel Size = (1000.134348869323730,-1000.134348869323730)
Corner Coordinates:
Upper Left  (-5568748.055, 5568748.055)
Lower Left  (-5568748.055,-5568748.055)
Upper Right ( 5568748.055, 5568748.055)
Lower Right ( 5568748.055,-5568748.055)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)

While this works with the operational MSG data, i'm not really sure if this is the correct approach. The scan direction of operational MSG satellites is always east->west (​MSG Level 1.5 Image Data Format Description page 97) but a reversed scan direction is possible and i don't know if the central pixel depends on the scan direction...

The attached patch ignores the scan direction.

[Even Rouault]

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to ​https://github.com/OSGeo/gdal/issues if it is not already reported there.