[netcdf-java] 2 fix on GRIB1 reader

Hello,

I work with martin desruisseaux on the geotoolkit.org project.

While trying to access datas from : http://www.globalmarinenet.com/grib_downloads.php I found a few bugs in the grib1 to netcdf metamodel.
Those are GRIB 1 files with LatLon projections.

grib/src/main/java/ucar/nc2/grib/grib1/Grib1Gds.java$LatLon

There was 2 bugs in this projection creation :

1 - a copy/paste error when calculating the latitude delta.

if (!Misc.closeEnough(deltaLat, calcDelta)) {
  log.debug("deltaLat != calcDeltaLat");
  deltaLat = calcDelta;  ----> was deltaLon = calcDelta;
}

2 - wrong value order when datas start by the lower values.

BEFORE
  if (deltaLat != GribNumbers.UNDEFINED) {
        deltaLat *= scale3; // undefined for thin grids
        if (la2 < la1) deltaLat *= -1.0;
  } else deltaLat = calcDelta;

AFTER
    if (deltaLat != GribNumbers.UNDEFINED) {
        deltaLat *= scale3; // undefined for thin grids
        if (la2 < la1) {
            //flip declaration order
            float latemp = la1;
            la1 = la2;
            la2 = latemp;
            calcDelta *= -1.0;

//we must also consider the cell corner, since we flipped the order //we should specify that the true value is at the BOTTOM-LEFT corner //but we can't show this information so we make a one cell displacement
            //to move the value on a TOP-LEFT corner.
            la1 -= calcDelta;
            la2 -= calcDelta;
        }
      } else {
          deltaLat = calcDelta;
      }



I attached a snapshot of the visualized datas before and after fix over the indian ocean and the patched java class.


Johann Sorel
Geomatys



Attachment: after.jpg
Description: JPEG image

Attachment: before.jpg
Description: JPEG image

/*
 * Copyright (c) 1998 - 2011. University Corporation for Atmospheric 
Research/Unidata
 * Portions of this software were developed by the Unidata Program at the
 * University Corporation for Atmospheric Research.
 *
 * Access and use of this software shall impose the following obligations
 * and understandings on the user. The user is granted the right, without
 * any fee or cost, to use, copy, modify, alter, enhance and distribute
 * this software, and any derivative works thereof, and its supporting
 * documentation for any purpose whatsoever, provided that this entire
 * notice appears in all copies of the software, derivative works and
 * supporting documentation.  Further, UCAR requests that the user credit
 * UCAR/Unidata in any publications that result from the use of this
 * software or in any product that includes this software. The names UCAR
 * and/or Unidata, however, may not be used in any advertising or publicity
 * to endorse or promote any products or commercial entity unless specific
 * written permission is obtained from UCAR/Unidata. The user also
 * understands that UCAR/Unidata is not obligated to provide the user with
 * any support, consulting, training or assistance of any kind with regard
 * to the use, operation and performance of this software nor to provide
 * the user with any updates, revisions, new versions or "bug fixes."
 *
 * THIS SOFTWARE IS PROVIDED BY UCAR/UNIDATA "AS IS" AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL UCAR/UNIDATA BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING
 * FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
 * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
 * WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE.
 */

package ucar.nc2.grib.grib1;

import ucar.nc2.grib.GribNumbers;
import ucar.nc2.grib.GdsHorizCoordSys;
import ucar.nc2.grib.QuasiRegular;
import ucar.nc2.util.Misc;
import ucar.unidata.geoloc.*;
import ucar.unidata.geoloc.projection.LatLonProjection;
import ucar.unidata.geoloc.projection.Stereographic;

import java.util.Formatter;

/**
 * GRIB1 GDS - subclass for each type
 * <p/>
 * <ul>
 * <li>0 Latitude/longitude grid â?? equidistant cylindrical or Plate Carrée 
projection
 * <li>1 Mercator projection
 * <li>2 Gnomonic projection
 * <li>3 Lambert conformal, secant or tangent, conic or bi-polar, projection
 * <li>4 Gaussian latitude/longitude grid
 * <li>5 Polar stereographic projection
 * <li>6 Universal Transverse Mercator (UTM) projection
 * <li>7 Simple polyconic projection
 * <li>8 Albers equal-area, secant or tangent, conic or bi-polar, projection
 * <li>9 Millerâ??s cylindrical projection
 * <li>10 Rotated latitude/longitude grid
 * <li>13 Oblique Lambert conformal, secant or tangent, conic or bi-polar, 
projection
 * <li>14 Rotated Gaussian latitude/longitude grid
 * <li>20 Stretched latitude/longitude grid
 * <li>24 Stretched Gaussian latitude/longitude grid
 * <li>30 Stretched and rotated latitude/longitude grids
 * <li>34 Stretched and rotated Gaussian latitude/longitude grids
 * <li>50 Spherical harmonic coefficients
 * <li>60 Rotated spherical harmonic coefficients
 * <li>70 Stretched spherical harmonics
 * <li>80 Stretched and rotated spherical harmonic coefficients
 * <li>90 Space view, perspective or orthographic
 * </ul>
 *
 * @author John
 * @since 9/3/11
 */
public abstract class Grib1Gds {
  static private final org.slf4j.Logger log = 
org.slf4j.LoggerFactory.getLogger(Grib1Gds.class);

  public static Grib1Gds factory(int template, byte[] data) {
    switch (template) {
      case 0:
        return new LatLon(data, 0);
      case 1:
        return new Mercator(data, 1);
      case 3:
        return new LambertConformal(data, 3);
      case 4:
        return new GaussianLatLon(data, 4);
      case 5:
        return new PolarStereographic(data, 5);
      case 10:
        return new RotatedLatLon(data, 10);
      default:
        throw new UnsupportedOperationException("Unsupported GDS type = " + 
template);
    }
  }

  ///////////////////////////////////////////////////
  private static final float scale3 = (float) 1.0e-3;
  private static final float scale6 = (float) 1.0e-6;

  protected final byte[] data;
  protected int[] nptsInLine; // thin grids, else null

  public int template;
  public float earthRadius, majorAxis, minorAxis;
  public int earthShape;
  protected int nx, ny;
  public int scanMode, resolution;
  protected int lastOctet;

  protected Grib1Gds(int template) {
    this.template = template;
    this.data = null;
  }

  public Grib1Gds(byte[] data, int template) {
    this.data = data;
    this.template = template;

    nx = getOctet2(7);
    ny = getOctet2(9);
  }

  public byte[] getRawBytes() {
    return data;
  }

  public int[] getNptsInLine() {
    return nptsInLine;
  }

  public void setNptsInLine(int[] nptsInLine) {
    this.nptsInLine = nptsInLine;
  }

  protected int getOctet(int index) {
    if (index > data.length) return GribNumbers.UNDEFINED;
    return data[index - 1] & 0xff;
  }

  protected int getOctet2(int start) {
    return GribNumbers.int2(getOctet(start), getOctet(start + 1));
  }

  protected int getOctet3(int start) {
    return GribNumbers.int3(getOctet(start), getOctet(start + 1), 
getOctet(start + 2));
  }

  protected int getOctet4(int start) {
    return GribNumbers.int4(getOctet(start), getOctet(start + 1), 
getOctet(start + 2), getOctet(start + 3));
  }

  /* protected float getScaledValue(int start) {
    int scaleFactor = getOctet(start);
    int scaleValue = getOctet4(start + 1);
    if (scaleFactor != 0)
      return (float) (scaleValue / Math.pow(10, scaleFactor));
    else
      return (float) scaleValue;
  } */

  public boolean isLatLon() {
    return false;
  }

  /*
  Flag/Code table 8 â?? Scanning mode
 Bit Value Meaning
  1   0   Points scan in +i direction
      1   Points scan in â??i direction
  2   0   Points scan in â??j direction
      1   Points scan in +j direction
  3   0   Adjacent points in i direction are consecutive
      1   Adjacent points in j direction are consecutive
   */
  public int getScanMode() {
    return scanMode;
  }

  public int getResolution() {
    return resolution;
  }

  public int getNxRaw() {
    return nx;
  }

  public int getNyRaw() {
    return ny;
  }

  ///////// thin grid crapola
  /*
    For data on a quasi-regular grid, in which all the rows or columns do not 
necessarily have the same number of grid
    points, either Ni (octets 7â??8) or Nj (octets 9â??10) and the 
corresponding Di (octets 24â??25) or Dj (octets 26â??27) shall be
    coded with all bits set to 1 (missing); the actual number of points along 
each parallel or meridian shall be coded.
  */

  // number of points along nx, adjusted for thin grid
  public int getNx() {
    if (nptsInLine == null || nx > 0) return nx;
    return QuasiRegular.getMax(nptsInLine);
  }

  // number of points along ny, adjusted for thin grid
  public int getNy() {
    if (nptsInLine == null || ny > 0) return ny;
    return QuasiRegular.getMax(nptsInLine);
  }

  // delta in x direction
  public float getDx() {
    return getDxRaw();
  }

  // delta in y direction
  public float getDy() {
    return getDyRaw();
  }

  public abstract float getDxRaw();

  public abstract float getDyRaw();

  public abstract GdsHorizCoordSys makeHorizCoordSys();

  public abstract void testHorizCoordSys(Formatter f);

  public String getNameShort() {
    String className = getClass().getName();
    int pos = className.lastIndexOf("$");
    return className.substring(pos + 1);
  }

  @Override
  public boolean equals(Object o) {
    if (this == o) return true;
    if (o == null || getClass() != o.getClass()) return false;

    Grib1Gds Grib1Gds = (Grib1Gds) o;
    if (nx != Grib1Gds.nx) return false;
    if (ny != Grib1Gds.ny) return false;
    if (template != Grib1Gds.template) return false;

    return true;
  }

  @Override
  public int hashCode() {
    int result = template;
    result = 31 * result + nx;
    result = 31 * result + ny;
    return result;
  }

  protected int hashCode = 0;

  /*
  Code Table Code table 3.2 - Shape of the Earth (3.2)
      0: Earth assumed spherical with radius = 6 367 470.0 m
      1: Earth assumed spherical with radius specified (in m) by data producer
      2: Earth assumed oblate spheroid with size as determined by IAU in 1965 
(major axis = 6 378 160.0 m, minor axis = 6 356 775.0 m, f = 1/297.0)
      3: Earth assumed oblate spheroid with major and minor axes specified (in 
km) by data producer
      4: Earth assumed oblate spheroid as defined in IAG-GRS80 model (major 
axis = 6 378 137.0 m, minor axis = 6 356 752.314 m, f = 1/298.257 222 101)
      5: Earth assumed represented by WGS84 (as used by ICAO since 1998)
      6: Earth assumed spherical with radius of 6 371 229.0 m
      7: Earth assumed oblate spheroid with major or minor axes specified (in 
m) by data producer
      8: Earth model assumed spherical with radius of 6 371 200 m, but the 
horizontal datum of the resulting
         latitude/longitude field is the WGS84 reference frame
  */
  protected Earth getEarth() {
    switch (earthShape) {
      case 0:
        return new Earth(6367470.0);
      case 1:
        return new Earth(earthRadius);
      case 2:
        return EarthEllipsoid.IAU;
      case 3:
        return new EarthEllipsoid("Grib2 Type 3", -1, majorAxis, minorAxis, 0);
      case 4:
        return EarthEllipsoid.IAG_GRS80;
      case 5:
        return EarthEllipsoid.WGS84;
      case 6:
        return new Earth(6371229.0);
      case 7:
        return new EarthEllipsoid("Grib2 Type 37", -1, majorAxis * 1000, 
minorAxis * 1000, 0);
      case 8:
        return new Earth(6371200.0);
      default:
        return new Earth();
    }
  }

  /*
    Grid definition â??   latitude/longitude grid (or equidistant cylindrical, 
or Plate Carrée)
    Octet No. Contents
    7â??8 Ni â?? number of points along a parallel
    9â??10 Nj â?? number of points along a meridian
    11â??13 La1 â?? latitude of first grid point
    14â??16 Lo1 â?? longitude of first grid point
    17 Resolution and component flags (see Code table 7)
    18â??20 La2 â?? latitude of last grid point
    21â??23 Lo2 â?? longitude of last grid point
    24â??25 Di â?? i direction increment
    26â??27 Dj â?? j direction increment
    28 Scanning mode (flags â?? see Flag/Code table 8)
    29â??32 Set to zero (reserved)
    33â??35 Latitude of the southern pole in millidegrees (integer) / Latitude 
of pole of stretching in millidegrees (integer)
    36â??38 Longitude of the southern pole in millidegrees (integer) / 
Longitude of pole of stretching in millidegrees (integer)
    39â??42 Angle of rotation (represented in the same way as the reference 
value) / Stretching factor (representation as for the reference value)
    43â??45 Latitude of pole of stretching in millidegrees (integer)
    46â??48 Longitude of pole of stretching in millidegrees (integer)
    49â??52 Stretching factor (representation as for the reference value)
*/
  public static class LatLon extends Grib1Gds {
    public float la1, lo1, la2, lo2, deltaLon, deltaLat;

    public LatLon(int template) {
      super(template);
    }

    LatLon(byte[] data, int template) {
      super(data, template);

      la1 = getOctet3(11) * scale3;
      lo1 = getOctet3(14) * scale3;
      resolution = getOctet(17);
      la2 = getOctet3(18) * scale3;
      lo2 = getOctet3(21) * scale3;

      if (lo2 < lo1) lo2 += 360.0F;

      deltaLon = getOctet2(24);
      float calcDelta = (lo2 - lo1) / (nx-1); // more accurate - deltaLon may 
have roundoff
      if (deltaLon != GribNumbers.UNDEFINED) {
          deltaLon *= scale3; // undefined for thin grids
      } else {
          deltaLon = calcDelta;
      }
      
      if (!Misc.closeEnough(deltaLon, calcDelta)) {
        log.debug("deltaLon != calcDeltaLon");
        deltaLon = calcDelta;
      }

      deltaLat = getOctet2(26);
      calcDelta = (la2 - la1) / (ny-1); // more accurate - deltaLon may have 
roundoff
      if (deltaLat != GribNumbers.UNDEFINED) {
        deltaLat *= scale3; // undefined for thin grids
        if (la2 < la1) {
            //flip declaration order
            float latemp = la1;
            la1 = la2;
            la2 = latemp;
            calcDelta *= -1.0;
            
            //we must also consider the cell corner, since we flipped the order
            //we should specify that the true value is at the BOTTOM-LEFT corner
            //but we can't show this information so we make a one cell 
displacement
            //to move the value on a TOP-LEFT CORNER.
            la1 -= calcDelta;
            la2 -= calcDelta;
        }
      } else {
          deltaLat = calcDelta;
      }
      
      if (!Misc.closeEnough(deltaLat, calcDelta)) {
        log.debug("deltaLat != calcDeltaLat");
        deltaLat = calcDelta;
      }

      scanMode = (byte) getOctet(28);

      lastOctet = 28;
    }

    @Override
    public void setNptsInLine(int[] nptsInLine) {
      super.setNptsInLine(nptsInLine);
      // now that we know this, we may have to recalc some stuff
      int n = QuasiRegular.getMax(nptsInLine);
      if (nx < 0) {
        deltaLon = (lo2 - lo1) / (n-1);
      } if (ny < 0) {
        deltaLat = (la2 - la1) / (n-1);
      }
    }


    @Override
    public boolean isLatLon() {
      return true;
    }

    @Override
    public float getDx() {
      if (nptsInLine == null || deltaLon != GribNumbers.UNDEFINED) return 
deltaLon;
      return (float) (lo2 - lo1) / (getNx() - 1);
  }

    @Override
    public float getDy() {
      if (nptsInLine == null || deltaLat != GribNumbers.UNDEFINED) return 
deltaLat;
      return (float) (la2 - la1) / (getNy() - 1);
  }

    @Override
    public float getDxRaw() {
       return deltaLon;
    }

    @Override
    public float getDyRaw() {
      return deltaLat;
    }

    @Override
    public boolean equals(Object o) {
      if (this == o) return true;
      if (o == null || getClass() != o.getClass()) return false;
      if (!super.equals(o)) return false;

      LatLon other = (LatLon) o;
      if (!Misc.closeEnough(la1, other.la1)) return false;
      if (!Misc.closeEnough(lo1, other.lo1)) return false;
      if (!Misc.closeEnough(deltaLat, other.deltaLat)) return false;
      if (!Misc.closeEnough(deltaLon, other.deltaLon)) return false;
      return true;
    }

    @Override
    public int hashCode() {
      if (hashCode == 0) {
        int result = super.hashCode();
        result = 31 * result + (la1 != +0.0f ? Float.floatToIntBits(la1) : 0); 
// LOOK this is an exact comparision
        result = 31 * result + (lo1 != +0.0f ? Float.floatToIntBits(lo1) : 0);
        result = 31 * result + (deltaLon != +0.0f ? 
Float.floatToIntBits(deltaLon) : 0);
        result = 31 * result + (deltaLat != +0.0f ? 
Float.floatToIntBits(deltaLat) : 0);
        hashCode = result;
      }
      return hashCode;
    }

    @Override
    public String toString() {
      return "LatLon{" +
              "la1=" + la1 +
              ", lo1=" + lo1 +
              ", la2=" + la2 +
              ", lo2=" + lo2 +
              ", deltaLon=" + deltaLon +
              ", deltaLat=" + deltaLat +
              '}';
    }

    public GdsHorizCoordSys makeHorizCoordSys() {
      LatLonProjection proj = new LatLonProjection();
      ProjectionPoint startP = proj.latLonToProj(new LatLonPointImpl(la1, lo1));
      double startx = startP.getX();
      double starty = startP.getY();
      return new GdsHorizCoordSys(getNameShort(), template, 0, scanMode, proj, 
startx, getDx(), starty, getDy(), getNx(), getNy());
    }

    public void testHorizCoordSys(Formatter f) {
      GdsHorizCoordSys cs = makeHorizCoordSys();
      double Lo2 = lo2;
      if (Lo2 < lo1) Lo2 += 360;
      LatLonPointImpl startLL = new LatLonPointImpl(la1, lo1);
      LatLonPointImpl endLL = new LatLonPointImpl(la2, lo2);

      f.format("%s testProjection%n", getClass().getName());
      f.format("  start at latlon= %s%n", startLL);
      f.format("    end at latlon= %s%n", endLL);

      ProjectionPointImpl endPP = (ProjectionPointImpl) 
cs.proj.latLonToProj(endLL, new ProjectionPointImpl());
      f.format("   start at proj coord= %s%n", new 
ProjectionPointImpl(cs.startx, cs.starty));
      f.format("     end at proj coord= %s%n", endPP);

      double endx = cs.startx + (getNx() - 1) * cs.dx;
      double endy = cs.starty + (getNy() - 1) * cs.dy;
      f.format("   should end at x= (%f,%f)%n", endx, endy);
    }
  }

  /*
  Grid definition â??   Gaussian latitude/longitude grid (including rotated, 
stretched or stretched and rotated)
 Octet    Contents
    7â??8     Ni â?? number of points along a parallel
    9â??10    Nj â?? number of points along a meridian
    11â??13   La1 â?? latitude of first grid point
    14â??16   Lo1 â?? longitude of first grid point
    17      Resolution and component flags (see Code table 7)
    18â??20   La2 â?? latitude of last grid point
    21â??23   Lo2 â?? longitude of last grid point
    24â??25   Di â?? i direction increment
    26â??27   N â?? number of parallels between a pole and the equator
    28      Scanning mode (flags â?? see Flag/Code table 8)
    29â??32   Set to zero (reserved)
    33â??35   Latitude of the southern pole in millidegrees (integer)
            Latitude of pole of stretching in millidegrees (integer)
    36â??38   Longitude of the southern pole in millidegrees (integer)
            Longitude of pole of stretching in millidegrees (integer)
    39â??42   Angle of rotation (represented in the same way as the reference 
value)
            Stretching factor (representation as for the reference value)
    43â??45   Latitude of pole of stretching in millidegrees (integer)
    46â??48   Longitude of pole of stretching in millidegrees (integer)
    49â??52   Stretching factor (representation as for the reference value)
   */
  public static class GaussianLatLon extends LatLon {
    int nparellels;
    public float latSouthPole, lonSouthPole, rotAngle, latPole, lonPole, 
stretchFactor;

    GaussianLatLon(byte[] data, int template) {
      super(data, template);
      nparellels = getOctet2(26);

      if (data.length > 32) {
        latSouthPole = getOctet3(33) * scale3;
        lonSouthPole = getOctet3(36) * scale3;
        rotAngle = getOctet4(39) * scale3;
        latPole = getOctet3(43) * scale3;
        lonPole = getOctet3(46) * scale3;
        stretchFactor = getOctet4(49) * scale3;
      }

      lastOctet = 52;
    }

    @Override
    public GdsHorizCoordSys makeHorizCoordSys() {
      GdsHorizCoordSys hc = super.makeHorizCoordSys();
      hc.setGaussianLats(nparellels, la1, la2);
      return hc;
    }

    @Override
    public boolean equals(Object o) {
      if (this == o) return true;
      if (o == null || getClass() != o.getClass()) return false;
      if (!super.equals(o)) return false;

      GaussianLatLon that = (GaussianLatLon) o;

      if (nparellels != that.nparellels) return false;

      return true;
    }

    @Override
    public int hashCode() {
      int result = super.hashCode();
      result = 31 * result + nparellels;
      return result;
    }

    @Override
    public String toString() {
      final StringBuilder sb = new StringBuilder();
      sb.append(super.toString());
      sb.append("\nGaussianLatLon");
      sb.append("{nparellels=").append(nparellels);
      sb.append(", latSouthPole=").append(latSouthPole);
      sb.append(", lonSouthPole=").append(lonSouthPole);
      sb.append(", rotAngle=").append(rotAngle);
      sb.append(", latPole=").append(latPole);
      sb.append(", lonPole=").append(lonPole);
      sb.append(", stretchFactor=").append(stretchFactor);
      sb.append('}');
      return sb.toString();
    }
  }

  /*
Grid definition â??   polar stereographic
 Octet No. Contents
 7â??8    Nx â?? number of points along x-axis
 9â??10   Ny â?? number of points along y-axis
 11â??13  La1 â?? latitude of first grid point
 14â??16  Lo1 â?? longitude of first grid point
 17     Resolution and component flags (see Code table 7)
 18â??20  LoV â?? orientation of the grid; i.e. the longitude value of the 
meridian which is parallel to the y-axis (or columns
              of the grid) along which latitude increases as the Y-coordinate 
increases (the orientation longitude may or may not appear on a particular grid)
 21â??23  Dx â?? X-direction grid length (see Note 2)
 24â??26  Dy â?? Y-direction grid length (see Note 2)
 27     Projection centre flag (see Note 5)
 28     Scanning mode (flags â?? see Flag/Code table 8)
 29â??32  Set to zero (reserved)
  */
  public static class PolarStereographic extends Grib1Gds {
    public float la1, lo1, lov, dX, dY;
    public int projCenterFlag;
    private float lad = (float) 60.0; // LOOK

    protected PolarStereographic(int template) {
      super(template);
    }

    PolarStereographic(byte[] data, int template) {
      super(data, template);

      la1 = getOctet3(11) * scale3;
      lo1 = getOctet3(14) * scale3;
      resolution =  getOctet(17);
      lov = getOctet3(18) * scale3;

      dX = getOctet3(21) * scale3;
      dY = getOctet3(24) * scale3;

      projCenterFlag =  getOctet(27);
      scanMode =  getOctet(28);

      lastOctet = 28;
    }

    @Override
    public boolean equals(Object o) {
      if (this == o) return true;
      if (o == null || getClass() != o.getClass()) return false;
      if (!super.equals(o)) return false;

      PolarStereographic that = (PolarStereographic) o;

      if (Float.compare(that.dX, dX) != 0) return false;
      if (Float.compare(that.dY, dY) != 0) return false;
      if (Float.compare(that.la1, la1) != 0) return false;
      if (Float.compare(that.lo1, lo1) != 0) return false;
      if (Float.compare(that.lov, lov) != 0) return false;
      if (projCenterFlag != that.projCenterFlag) return false;

      return true;
    }

    @Override
    public int hashCode() {
      if (hashCode == 0) {
        int result = super.hashCode();
        result = 31 * result + (la1 != +0.0f ? Float.floatToIntBits(la1) : 0);
        result = 31 * result + (lo1 != +0.0f ? Float.floatToIntBits(lo1) : 0);
        result = 31 * result + (lov != +0.0f ? Float.floatToIntBits(lov) : 0);
        result = 31 * result + (dX != +0.0f ? Float.floatToIntBits(dX) : 0);
        result = 31 * result + (dY != +0.0f ? Float.floatToIntBits(dY) : 0);
        result = 31 * result + (int) projCenterFlag;
        hashCode = result;
      }
      return hashCode;
    }

    @Override
    public float getDxRaw() {
      return dX;
    }

    @Override
    public float getDyRaw() {
      return dY;
    }

    public GdsHorizCoordSys makeHorizCoordSys() {
      boolean northPole = (projCenterFlag & 128) == 0;
      double latOrigin = northPole ? 90.0 : -90.0;

      // Why the scale factor?. according to GRIB docs:
      // "Grid lengths are in units of meters, at the 60 degree latitude circle 
nearest to the pole"
      // since the scale factor at 60 degrees = k = 2*k0/(1+sin(60))  
[Snyder,Working Manual p157]
      // then to make scale = 1 at 60 degrees, k0 = (1+sin(60))/2 = .933
      double scale;
      if (Double.isNaN(lad)) { // LOOK ??
        scale = 0.9330127018922193;
      } else {
        scale = (1.0 + Math.sin(Math.toRadians( Math.abs(lad)))) / 2;
      }

      ProjectionImpl proj = null;

      Earth earth = getEarth();
      if (earth.isSpherical()) {
        proj = new Stereographic(latOrigin, lov, scale);
      } else {
        proj = new 
ucar.unidata.geoloc.projection.proj4.StereographicAzimuthalProjection(
                latOrigin, lov, scale, lad, 0.0, 0.0, earth);
      }

      ProjectionPointImpl start = (ProjectionPointImpl) proj.latLonToProj(new 
LatLonPointImpl(la1, lo1));
      return new GdsHorizCoordSys(getNameShort(), template, 0, scanMode, proj, 
start.getX(), getDx(), start.getY(), getDy(), getNx(), getNy());
    }

    public void testHorizCoordSys(Formatter f) {
      GdsHorizCoordSys cs = makeHorizCoordSys();
      f.format("%s testProjection %s%n", getClass().getName(), 
cs.proj.getClass().getName());

      double endx = cs.startx + (getNx() - 1) * cs.dx;
      double endy = cs.starty + (getNy() - 1) * cs.dy;
      ProjectionPointImpl endPP = new ProjectionPointImpl(endx, endy);
      f.format("   start at proj coord= %s%n", new 
ProjectionPointImpl(cs.startx, cs.starty));
      f.format("     end at proj coord= %s%n", endPP);

      LatLonPointImpl startLL = new LatLonPointImpl(la1, lo1);
      LatLonPoint endLL = cs.proj.projToLatLon(endPP, new LatLonPointImpl());

      f.format("  start at latlon= %s%n", startLL);
      f.format("    end at latlon= %s%n", endLL);
    }

  }

  /*
  Grid definition â??   Lambert conformal, secant or tangent, conic or bi-polar 
(normal or oblique), or
  Albers  equal-area,  secant  or  tangent,  conic  or  bi-polar  (normal  or  
oblique), projection
  Octet No. Contents
  7â??8   Nx â?? number of points along x-axis
  9â??10  Ny â?? number of points along y-axis
  11â??13 La1 â?? latitude of first grid point
  14â??16 Lo1â?? longitude of first grid point
  17    Resolution and component flags (see Code table 7)
  18â??20 LoV â??   orientation of the grid; i.e. the east longitude value of 
the meridian which is parallel to the y-axis
    (or columns of the grid) along which latitude increases as the y-coordinate 
increases (the orientation longitude may
    or may not appear on a particular grid)
  21â??23 Dx â?? x-direction grid length (see Note 2)
  24â??26 Dy â?? y-direction grid length (see Note 2)
  27    Projection centre flag (see Note 5)
  28    Scanning mode (flags â?? see Flag/Code table 8)
  29â??31 Latin 1 â?? first latitude from the pole at which the secant cone 
cuts the sphere
  32â??34 Latin 2 â?? second latitude from the pole at which the secant cone 
cuts the sphere
  35â??37 Latitude of the southern pole in millidegrees (integer)
  38â??40 Longitude of the southern pole in millidegrees (integer)
  41â??42 Set to zero (reserved)
   */
  public static class LambertConformal extends Grib1Gds {
    public float la1, lo1, lov, lad, dX, dY, latin1, latin2, latSouthPole, 
lonSouthPole;
    public int projCenterFlag;

    // private int hla1, hlo1, hlov, hlad, hdX, hdY, hlatin1, hlatin2; // 
hasheesh

    LambertConformal(byte[] data, int template) {
      super(data, template);

      la1 = getOctet3(11) * scale3;
      lo1 = getOctet3(14) * scale3;
      resolution =  getOctet(17);
      lov = getOctet3(18) * scale3;

      dX = getOctet3(21) * scale3;
      dY = getOctet3(24) * scale3;

      projCenterFlag =  getOctet(27);
      scanMode = getOctet(28);

      latin1 = getOctet3(29) * scale3;
      latin2 = getOctet3(32) * scale3;
      latSouthPole = getOctet3(35) * scale3;
      lonSouthPole = getOctet3(38) * scale3;
    }

    @Override
    public float getDxRaw() {
      return dX;
    }

    @Override
    public float getDyRaw() {
      return dY;
    }

    @Override
    public boolean equals(Object o) {
      if (this == o) return true;
      if (o == null || getClass() != o.getClass()) return false;
      if (!super.equals(o)) return false;

      LambertConformal that = (LambertConformal) o;

      if (Float.compare(that.dX, dX) != 0) return false;
      if (Float.compare(that.dY, dY) != 0) return false;
      if (Float.compare(that.la1, la1) != 0) return false;
      if (Float.compare(that.lad, lad) != 0) return false;
      if (Float.compare(that.latSouthPole, latSouthPole) != 0) return false;
      if (Float.compare(that.latin1, latin1) != 0) return false;
      if (Float.compare(that.latin2, latin2) != 0) return false;
      if (Float.compare(that.lo1, lo1) != 0) return false;
      if (Float.compare(that.lonSouthPole, lonSouthPole) != 0) return false;
      if (Float.compare(that.lov, lov) != 0) return false;
      if (projCenterFlag != that.projCenterFlag) return false;

      return true;
    }

    @Override
    public int hashCode() {
      int result = super.hashCode();
      result = 31 * result + (la1 != +0.0f ? Float.floatToIntBits(la1) : 0);
      result = 31 * result + (lo1 != +0.0f ? Float.floatToIntBits(lo1) : 0);
      result = 31 * result + (lov != +0.0f ? Float.floatToIntBits(lov) : 0);
      result = 31 * result + (lad != +0.0f ? Float.floatToIntBits(lad) : 0);
      result = 31 * result + (dX != +0.0f ? Float.floatToIntBits(dX) : 0);
      result = 31 * result + (dY != +0.0f ? Float.floatToIntBits(dY) : 0);
      result = 31 * result + (latin1 != +0.0f ? Float.floatToIntBits(latin1) : 
0);
      result = 31 * result + (latin2 != +0.0f ? Float.floatToIntBits(latin2) : 
0);
      result = 31 * result + (latSouthPole != +0.0f ? 
Float.floatToIntBits(latSouthPole) : 0);
      result = 31 * result + (lonSouthPole != +0.0f ? 
Float.floatToIntBits(lonSouthPole) : 0);
      result = 31 * result + projCenterFlag;
      return result;
    }

    /*
  @Override
  public boolean equals(Object o) {
    if (this == o) return true;
    if (o == null || getClass() != o.getClass()) return false;
    if (!super.equals(o)) return false;

    LambertConformal that = (LambertConformal) o;

    if (hdX != that.hdX) return false;
    if (hdY != that.hdY) return false;
    if (hla1 != that.hla1) return false;
    if (hlad != that.hlad) return false;
    if (hlatin1 != that.hlatin1) return false;
    if (hlatin2 != that.hlatin2) return false;
    if (hlo1 != that.hlo1) return false;
    if (hlov != that.hlov) return false;

    return true;
  }

  @Override
  public int hashCode() {
    if (hashCode == 0) {
      int result = super.hashCode();
      result = 31 * result + hla1;
      result = 31 * result + hlo1;
      result = 31 * result + hlov;
      result = 31 * result + hlad;
      result = 31 * result + hdX;
      result = 31 * result + hdY;
      result = 31 * result + hlatin1;
      result = 31 * result + hlatin2;
      hashCode = result;
    }
    return hashCode;
  }

  private static int round(int a) { // NCEP rounding (!)
    return (a + 5) / 10;
  }  */

    public GdsHorizCoordSys makeHorizCoordSys() {
      ProjectionImpl proj = null;

      Earth earth = getEarth();
      if (earth.isSpherical()) {
        proj = new ucar.unidata.geoloc.projection.LambertConformal(latin1, lov, 
latin1, latin2, 0.0, 0.0, earth.getEquatorRadius() * .001);
      } else {
        proj = new 
ucar.unidata.geoloc.projection.proj4.LambertConformalConicEllipse(
                latin1, lov, latin1, latin2, 0.0, 0.0, earth);
      }

      LatLonPointImpl startLL = new LatLonPointImpl(la1, lo1);
      ProjectionPointImpl start = (ProjectionPointImpl) 
proj.latLonToProj(startLL);
      return new GdsHorizCoordSys(getNameShort(), template, 0, scanMode, proj, 
start.getX(), getDx(), start.getY(), getDy(), getNx(), getNy());
    }

    public void testHorizCoordSys(Formatter f) {
      GdsHorizCoordSys cs = makeHorizCoordSys();
      f.format("%s testProjection %s%n", getClass().getName(), 
cs.proj.getClass().getName());

      double endx = cs.startx + (getNx() - 1) * cs.dx;
      double endy = cs.starty + (getNy() - 1) * cs.dy;
      ProjectionPointImpl endPP = new ProjectionPointImpl(endx, endy);
      f.format("   start at proj coord= %s%n", new 
ProjectionPointImpl(cs.startx, cs.starty));
      f.format("     end at proj coord= %s%n", endPP);

      LatLonPointImpl startLL = new LatLonPointImpl(la1, lo1);
      LatLonPoint endLL = cs.proj.projToLatLon(endPP, new LatLonPointImpl());

      f.format("  start at latlon= %s%n", startLL);
      f.format("    end at latlon= %s%n", endLL);
    }

  }

  /*
  Grid definition â?? Mercator
  Octet  Contents
    7â??8   Ni â?? number of points along a parallel
    9â??10  Nj â?? number of points along a meridian
    11â??13 La1 â?? latitude of first grid point
    14â??16 Lo1 â?? longitude of first grid point
    17    Resolution and component flags (see Code table 7)
    18â??20 La2 â?? latitude of last grid point
    21â??23 Lo2 â?? longitude of last grid point
    24â??26 Latin â?? latitude(s) at which the Mercator projection cylinder 
intersects the Earth
    27    Set to zero (reserved)
    28    Scanning mode (flags â?? see Flag/Code table 8)
    29â??31 Di â?? longitudinal direction grid length (see Note 2)
    32â??34 Dj â?? latitudinal direction grid length (see Note 2)
    35â??42 Set to zero (reserved)
   */

  public static class Mercator extends Grib1Gds {
    public float la1, lo1, la2, lo2, latin, dX, dY;
    protected int lastOctet;

    Mercator(byte[] data, int template) {
      super(data, template);

      la1 = getOctet3(11) * scale3;
      lo1 = getOctet3(14) * scale3;
      resolution = getOctet(47);
      la2 = getOctet3(18) * scale3;
      lo2 = getOctet3(21) * scale3;
      latin = getOctet3(24) * scale3;

      scanMode = getOctet(28);

      dX = getOctet3(29) * scale3;  // km
      dY = getOctet3(32) * scale3;  // km

      lastOctet = 42;
    }

    @Override
    public float getDxRaw() {
      return dX;
    }

    @Override
    public float getDyRaw() {
      return dY;
    }

    @Override
    public boolean equals(Object o) {
      if (this == o) return true;
      if (o == null || getClass() != o.getClass()) return false;
      if (!super.equals(o)) return false;

      Mercator mercator = (Mercator) o;

      if (Float.compare(mercator.dX, dX) != 0) return false;
      if (Float.compare(mercator.dY, dY) != 0) return false;
      if (Float.compare(mercator.la1, la1) != 0) return false;
      if (Float.compare(mercator.latin, latin) != 0) return false;
      if (Float.compare(mercator.lo1, lo1) != 0) return false;

      return true;
    }

    @Override
    public int hashCode() {
      if (hashCode == 0) {
        int result = super.hashCode();
        result = 31 * result + (la1 != +0.0f ? Float.floatToIntBits(la1) : 0);
        result = 31 * result + (lo1 != +0.0f ? Float.floatToIntBits(lo1) : 0);
        result = 31 * result + (latin != +0.0f ? Float.floatToIntBits(latin) : 
0);
        result = 31 * result + (dX != +0.0f ? Float.floatToIntBits(dX) : 0);
        result = 31 * result + (dY != +0.0f ? Float.floatToIntBits(dY) : 0);
        hashCode = result;
      }
      return hashCode;
    }

    public GdsHorizCoordSys makeHorizCoordSys() {
      // put longitude origin at first point - doesnt actually matter
      // param par standard parallel (degrees). cylinder cuts earth at this 
latitude.
      Earth earth = getEarth();
      ucar.unidata.geoloc.projection.Mercator proj = new 
ucar.unidata.geoloc.projection.Mercator(lo1, latin, 0, 0, 
earth.getEquatorRadius() * .001);

      // find out where things start
      ProjectionPoint startP = proj.latLonToProj(new LatLonPointImpl(la1, lo1));
      double startx = startP.getX();
      double starty = startP.getY();

      return new GdsHorizCoordSys(getNameShort(), template, 0, scanMode, proj, 
startx, getDx(), starty, getDy(), getNx(), getNy());
    }

    public void testHorizCoordSys(Formatter f) {
      GdsHorizCoordSys cs = makeHorizCoordSys();
      double Lo2 = lo2;
      if (Lo2 < lo1) Lo2 += 360;
      LatLonPointImpl startLL = new LatLonPointImpl(la1, lo1);
      LatLonPointImpl endLL = new LatLonPointImpl(la2, lo2);

      f.format("%s testProjection%n", getClass().getName());
      f.format("  start at latlon= %s%n", startLL);
      f.format("    end at latlon= %s%n", endLL);

      ProjectionPointImpl endPP = (ProjectionPointImpl) 
cs.proj.latLonToProj(endLL, new ProjectionPointImpl());
      f.format("   start at proj coord= %s%n", new 
ProjectionPointImpl(cs.startx, cs.starty));
      f.format("     end at proj coord= %s%n", endPP);

      double endx = cs.startx + (getNx() - 1) * cs.dx;
      double endy = cs.starty + (getNy() - 1) * cs.dy;
      f.format("   should end at x= (%f,%f)%n", endx, endy);
    }

  }

  public static class RotatedLatLon extends LatLon {
    public float angleRotation; // Angle of rotation (represented in the same 
way as the reference value)
    public float latSouthPole; // Latitude of pole of stretching in 
millidegrees (integer)
    public float lonSouthPole;  // Longitude of pole of stretching in 
millidegrees (integer)
    protected int lastOctet;

    RotatedLatLon(byte[] data, int template) {
      super(data, template);

      latSouthPole = getOctet3(33) * scale3;
      lonSouthPole = getOctet3(36) * scale3;
      angleRotation = getOctet4(39) * scale3;

      lastOctet = 43;
    }

    public String getName() {
      return "Rotated latitude/longitude";
    }

    @Override
    public String toString() {
      final StringBuilder sb = new StringBuilder();
      sb.append(super.toString());
      sb.append("\nRotLatLon");
      sb.append("{angleRotation=").append(angleRotation);
      sb.append(", latSouthPole=").append(latSouthPole);
      sb.append(", lonSouthPole=").append(lonSouthPole);
      sb.append(", lastOctet=").append(lastOctet);
      sb.append('}');
      return sb.toString();
    }

    @Override
    public boolean equals(Object o) {
      if (this == o) return true;
      if (o == null || getClass() != o.getClass()) return false;
      if (!super.equals(o)) return false;

      RotatedLatLon other = (RotatedLatLon) o;
      if (!Misc.closeEnough(angleRotation, other.angleRotation)) return false;
      return true;
    }

    @Override
    public int hashCode() {
      if (hashCode == 0) {
        int result = super.hashCode();
        result = 31 * result + (angleRotation != +0.0f ? 
Float.floatToIntBits(angleRotation) : 0);
        hashCode = result;
      }
      return hashCode;
    }

    public GdsHorizCoordSys makeHorizCoordSys() {
      ucar.unidata.geoloc.projection.RotatedLatLon proj =
              new ucar.unidata.geoloc.projection.RotatedLatLon(latSouthPole, 
lonSouthPole, angleRotation);
      // LOOK dont transform - works for grib1 
Q:/cdmUnitTest/transforms/HIRLAMhybrid.grib
      // LatLonPoint startLL = proj.projToLatLon(new ProjectionPointImpl(lo1, 
la1));
      //double startx = startLL.getLongitude();
      //double starty = startLL.getLatitude();
      return new GdsHorizCoordSys(getNameShort(), template, 0, scanMode, proj, 
lo1, deltaLon, la1, deltaLat, nx, ny);
    }

    public void testHorizCoordSys(Formatter f) {
      GdsHorizCoordSys cs = makeHorizCoordSys();
      LatLonPoint startLL = cs.proj.projToLatLon(new ProjectionPointImpl(lo1, 
la1));
      LatLonPoint endLL = cs.proj.projToLatLon(new ProjectionPointImpl(lo2, 
la2));

      f.format("%s testProjection%n", getClass().getName());
      f.format("  start at latlon= %s%n", startLL);
      f.format("    end at latlon= %s%n", endLL);

      ProjectionPointImpl endPP = (ProjectionPointImpl) 
cs.proj.latLonToProj(endLL, new ProjectionPointImpl());
      f.format("   start at proj coord= %s%n", new 
ProjectionPointImpl(cs.startx, cs.starty));
      f.format("     end at proj coord= %s%n", endPP);

      double endx = cs.startx + (nx - 1) * cs.dx;
      double endy = cs.starty + (ny - 1) * cs.dy;
      f.format("   should end at x= (%f,%f)%n", endx, endy);
    }

  }

}
  • 2012 messages navigation, sorted by:
    1. Thread
    2. Subject
    3. Author
    4. Date
    5. ↑ Table Of Contents
  • Search the netcdf-java archives: