visad/src/main/java/ucar/nc2/iosp/grid/GridHorizCoordSys.java
private void makeMSGgeostationary() {
    double Lat0 = gds.getDouble(GridDefRecord.LAP);  // sub-satellite point lat
    double Lon0 = gds.getDouble(GridDefRecord.LOP);  // sub-satellite point lon

    //int nx = gds.getInt(GridDefRecord.NX);
    int ny = gds.getInt(GridDefRecord.NY);
    int x_off = gds.getInt(GridDefRecord.XP);  // sub-satellite point in grid lengths
    int y_off = gds.getInt(GridDefRecord.YP);

    double dx; // = gds.getDouble(GridDefRecord.DX);  // apparent diameter of earth in units of grid lengths
    double dy = gds.getDouble(GridDefRecord.DY);

    // per Simon Eliot 1/18/2010, there is a bug in Eumetsat grib files,
    // we need to "correct for ellipsoidal earth"
    // (Note we should check who the originating center is
    // "Originating_center" = "EUMETSAT Operation Centre" in the GRIB id (section 1))
    // although AFAIK, eumetsat is only one using this projection.
    if (dy < 2100) {
      dx = 1207;
      dy = 1203;
    } else {
      dx = 3622;
      dy = 3610;
    }

    // have to check both names because Grib1 and Grib2 used different names
    double major_axis = gds.getDouble(GridDefRecord.MAJOR_AXIS_EARTH);  // m
    if (Double.isNaN(major_axis))
      major_axis = gds.getDouble("major_axis_earth");

    double minor_axis = gds.getDouble(GridDefRecord.MINOR_AXIS_EARTH);  // m
    if (Double.isNaN(minor_axis))
      minor_axis = gds.getDouble("minor_axis_earth");
    // Nr = altitude of camera from center, in units of radius
    double nr = gds.getDouble(GridDefRecord.NR) * 1e-6; // altitude of the camera from the Earths centre, measured in units of the Earth (equatorial) radius

          // CFAC = 2^16 / {[2 * arcsine (10^6 / Nr)] / dx }
    double as = 2 * Math.asin(1.0/nr);
    double cfac = dx / as;
    double lfac = dy / as;

    // use km, so scale by the earth radius
    double scale_factor = (nr - 1) * major_axis / 1000; // this sets the units of the projection x,y coords in km
    double scale_x = scale_factor; // LOOK fake neg need scan value
    double scale_y = -scale_factor; // LOOK fake neg need scan value
    startx = scale_factor * (1 - x_off) / cfac;
    starty = scale_factor * (y_off - ny) / lfac;
    incrx = scale_factor/cfac;
    incry = scale_factor/lfac;

    attributes.add(new Attribute(GridCF.GRID_MAPPING_NAME, "MSGnavigation"));
    attributes.add(new Attribute(GridCF.LONGITUDE_OF_PROJECTION_ORIGIN, Lon0));
    attributes.add(new Attribute(GridCF.LATITUDE_OF_PROJECTION_ORIGIN, Lat0));
    //attributes.add(new Attribute("semi_major_axis", new Double(major_axis)));
    //attributes.add(new Attribute("semi_minor_axis", new Double(minor_axis)));
    attributes.add(new Attribute("height_from_earth_center", nr * major_axis));
    attributes.add(new Attribute("scale_x", scale_x));
    attributes.add(new Attribute("scale_y", scale_y));

    proj = new MSGnavigation(Lat0, Lon0, major_axis, minor_axis, nr * major_axis, scale_x, scale_y);

    if (GridServiceProvider.debugProj) {

      double Lo2 = gds.getDouble(GridDefRecord.LO2) + 360.0;
      double La2 = gds.getDouble(GridDefRecord.LA2);
      LatLonPointImpl endLL = new LatLonPointImpl(La2, Lo2);
      System.out.println("GridHorizCoordSys.makeMSGgeostationary end at latlon " + endLL);

      ProjectionPointImpl endPP =(ProjectionPointImpl) proj.latLonToProj(endLL);
      System.out.println("   end at proj coord " + endPP);

      double endx = 1 + getNx();
      double endy = 1 + getNy();
      System.out.println("   should be x=" + endx + " y=" + endy);
    }
  }
private void makeMSGgeostationary() {
    double Lat0 = gds.getDouble(GridDefRecord.LAP);  // sub-satellite point lat
    double Lon0 = gds.getDouble(GridDefRecord.LOP);  // sub-satellite point lon

    //int nx = gds.getInt(GridDefRecord.NX);
    int ny = gds.getInt(GridDefRecord.NY);
    int x_off = gds.getInt(GridDefRecord.XP);  // sub-satellite point in grid lengths
    int y_off = gds.getInt(GridDefRecord.YP);

    double dx; // = gds.getDouble(GridDefRecord.DX);  // apparent diameter of earth in units of grid lengths
    double dy = gds.getDouble(GridDefRecord.DY);

    // per Simon Eliot 1/18/2010, there is a bug in Eumetsat grib files,
    // we need to "correct for ellipsoidal earth"
    // (Note we should check who the originating center is
    // "Originating_center" = "EUMETSAT Operation Centre" in the GRIB id (section 1))
    // although AFAIK, eumetsat is only one using this projection.
    if (dy < 2100) {
      dx = 1207;
      dy = 1203;
    } else {
      dx = 3622;
      dy = 3610;
    }

    // have to check both names because Grib1 and Grib2 used different names
    double major_axis = gds.getDouble(GridDefRecord.MAJOR_AXIS_EARTH);  // m
    if (Double.isNaN(major_axis))
      major_axis = gds.getDouble("major_axis_earth");

    double minor_axis = gds.getDouble(GridDefRecord.MINOR_AXIS_EARTH);  // m
    if (Double.isNaN(minor_axis))
      minor_axis = gds.getDouble("minor_axis_earth");
    // Nr = altitude of camera from center, in units of radius
    double nr = gds.getDouble(GridDefRecord.NR) * 1e-6; // altitude of the camera from the Earths centre, measured in units of the Earth (equatorial) radius

          // CFAC = 2^16 / {[2 * arcsine (10^6 / Nr)] / dx }
    double as = 2 * Math.asin(1.0/nr);
    double cfac = dx / as;
    double lfac = dy / as;

    // use km, so scale by the earth radius
    double scale_factor = (nr - 1) * major_axis / 1000; // this sets the units of the projection x,y coords in km
    double scale_x = scale_factor; // LOOK fake neg need scan value
    double scale_y = -scale_factor; // LOOK fake neg need scan value
    startx = scale_factor * (1 - x_off) / cfac;
    starty = scale_factor * (y_off - ny) / lfac;
    incrx = scale_factor/cfac;
    incry = scale_factor/lfac;

    attributes.add(new Attribute(GridCF.GRID_MAPPING_NAME, "MSGnavigation"));
    attributes.add(new Attribute(GridCF.LONGITUDE_OF_PROJECTION_ORIGIN, Lon0));
    attributes.add(new Attribute(GridCF.LATITUDE_OF_PROJECTION_ORIGIN, Lat0));
    //attributes.add(new Attribute("semi_major_axis", new Double(major_axis)));
    //attributes.add(new Attribute("semi_minor_axis", new Double(minor_axis)));
    attributes.add(new Attribute("height_from_earth_center", nr * major_axis));
    attributes.add(new Attribute("scale_x", scale_x));
    attributes.add(new Attribute("scale_y", scale_y));

    proj = new MSGnavigation(Lat0, Lon0, major_axis, minor_axis, nr * major_axis, scale_x, scale_y);

    if (GridServiceProvider.debugProj) {

      double Lo2 = gds.getDouble(GridDefRecord.LO2) + 360.0;
      double La2 = gds.getDouble(GridDefRecord.LA2);
      LatLonPointImpl endLL = new LatLonPointImpl(La2, Lo2);
      System.out.println("GridHorizCoordSys.makeMSGgeostationary end at latlon " + endLL);

      ProjectionPointImpl endPP =(ProjectionPointImpl) proj.latLonToProj(endLL);
      System.out.println("   end at proj coord " + endPP);

      double endx = 1 + getNx();
      double endy = 1 + getNy();
      System.out.println("   should be x=" + endx + " y=" + endy);
    }
  }
Make a Eumetsat MSG "Normalized Geostationary Projection" projection.
Fake coordinates for now, then see if this can be generalized.

from  FM 92 GRIB-2 doc:

Grid Definition Template 3.90: Space view perspective or orthographic
Octet Number(s) Contents
15 Shape of the earth (see Code Table 3.2)
16 Scale factor of radius of spherical earth
17-20 Scaled value of radius of spherical earth
21 Scale factor of major axis of oblate spheroid earth
22-25 Scaled value of major axis of oblate spheroid earth
26 Scale factor of minor axis of oblate spheroid earth
27-30 Scaled value of minor axis of oblate spheroid earth
31-34 Nx - number of points along X-axis (columns)
35-38 Ny - number of points along Y-axis (rows or lines)
39-42 Lap - latitude of sub-satellite point
43-46 Lop - longitude of sub-satellite point
47 Resolution and component flags (see Flag Table 3.3)
48-51 dx - apparent diameter of Earth in grid lengths, in X-direction
52-55 dy - apparent diameter of Earth in grid lengths, in Y-direction
56-59 Xp - X-coordinate of sub-satellite point (in units of 10-3 grid length expressed as an integer)
60-63 Yp - Y-coordinate of sub-satellite point (in units of 10-3 grid length expressed as an integer)
64 Scanning mode (flags - see Flag Table 3.4)
65-68 Orientation of the grid; i.e., the angle between the increasing Y-axis and the meridian of the sub-satellite point in the direction of increasing latitude (see Note 3)
69-72 Nr - altitude of the camera from the Earths centre, measured in units of the Earth (equatorial) radius multiplied by a scale factor of 10 6 (see Notes 4 and 5)
73-76 Xo - X-coordinate of origin of sector image
77-80 Yo - Y-coordinate of origin of sector image

Notes:
(1) It is assumed that the satellite is at its nominal position, i.e., it is looking directly at its sub-satellite point.
(2) Octets 69-72 shall be set to all ones (missing) to indicate the orthographic view (from infinite distance)
(3) It is the angle between the increasing Y-axis and the meridian 180E if the sub-satellite point is the North Pole; or the meridian 0 if the sub-satellite point is the South Pole.
(4) The apparent angular size of the Earth will be given by 2 * Arcsin (10^6 )/Nr).
(5) For orthographic view from infinite distance, the value of Nr should be encoded as missing (all bits set to 1).
(6) The horizontal and vertical angular resolutions of the sensor (Rx and Ry), needed for navigation equation, can be calculated from the following:
Rx = 2 * Arcsin (106 )/Nr)/ dx
Ry = 2 * Arcsin (106 )/Nr)/ dy

=======

from  simon.elliott@eumetsat.int

For products on a single pixel resolution grid, the scan angle is 83.84333 E-6 rad.
So dx = 2 * arcsin(10e6/Nr) / 83.84333 E-6 = 3622.30, which encoded to the nearest integer is 3622.
This is correctly encoded in our products.

For products on a 3x3 pixel resolution grid, the scan angle is 3 * 83.84333 E-6 rad = 251.52999 E-6 rad.
So dx = 2 * arcsin(10e6/Nr) / 251.52999 E-6 = 1207.43, which encoded to the nearest integer is 1207.
This is correctly encoded in our products.

Due to the elliptical shape of the earth, the calculation is a bit different in the y direction (Nr is in multiples of
the equatorial radius, but the tangent point is much closer to the polar radius from the earth's centre.
Approximating that the tangent point is actually at the polar radius from the earth's centre:
The sine of the angle subtended by the Earths centre and the tangent point on the equator as seen from the spacecraft
= Rp / (( Nr * Re) / 10^6) = (Rp * 10^6) / (Re * Nr)

The angle subtended by the Earth equator as seen by the spacecraft is, by symmetry twice the inverse sine above,
= 2 * arcsine ((Rp * 10^6) / (Re * Nr))

For products on a single pixel resolution grid, the scan angle is 83.84333 E-6 rad.
So dy = 2 * arcsine ((Rp * 10^6) / (Re * Nr)) / 83.84333 E-6 = 3610.06, which encoded to the nearest integer is 3610.
This is currently encoded in our products as 3568.

For products on a 3x3 pixel resolution grid, the scan angle is 3 * 83.84333 E-6 rad = 251.52999 E-6 rad.
So dy = 2 * arcsine ((Rp * 10^6) / (Re * Nr)) / 251.52999 E-6 = 1203.35, which encoded to the nearest integer is 1203.
This is currently encoded in our products as 1189.

As you can see the dx and dy values we are using will lead to an error of around 1% in the y direction.
I will ensure that the values are corrected to those explained here (3610 and 1203) as soon as possible.