ucar.unidata.geoloc.projection.RotatedPole Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of cdm Show documentation
The NetCDF-Java Library is a Java interface to NetCDF files, as well as to many other types of scientific data formats.
The newest version!
/*
 * Copyright 1998-2009 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.
 */
/*
 *	This software is provided by the National Aeronatics and Space
 *	Administration Goddard Institute for Space Studies (NASA GISS) for full,
 *	free and open release. It is understood by the recipient/user that NASA
 *	assumes no liability for any errors contained in the code.  Although this
 *	software is released without conditions or restrictions in its use, it is
 *	expected that appropriate credit be given to its author and to NASA GISS
 *	should the software be included by the recipient in acknowledgments or
 *	credits in other product development.
 */
package ucar.unidata.geoloc.projection;

import ucar.nc2.constants.CF;
import ucar.unidata.geoloc.*;

/**
 * Rotated-pole longitude-latitude grid.
 *
 * This is probably the same as rotated lat lon, using matrix to do rotation.
 * Follows CF convention with "north_pole lat/lon", whereas RotatedLatLon uses
 * south pole. LOOK
 *
 * @author Robert Schmunk
 * @author jcaron
 */
public class RotatedPole extends ProjectionImpl {
  private static final double RAD_PER_DEG = Math.PI / 180.;
  private static final double DEG_PER_RAD = 1. / RAD_PER_DEG;
  private static boolean show = false;

  /*	Coordinates of north pole for rotated pole. */
  private final ProjectionPointImpl northPole;

  /*	Y-axis rotation matrix. */
  private final double[][] rotY = new double[3][3];

  /*	Z-axis rotation matrix. */
  private final double[][] rotZ = new double[3][3];

  /**
   * Default Constructor, needed for beans.
   */
  public RotatedPole() {
    this(0.0, 0.0);
  }

  /**
   * Constructor.
   * @param northPoleLat
   * @param northPoleLon
   */
  public RotatedPole(double northPoleLat, double northPoleLon) {
    super("RotatedPole", false);

    northPole = new ProjectionPointImpl(northPoleLon, northPoleLat);
    buildRotationMatrices();

    addParameter(CF.GRID_MAPPING_NAME, CF.ROTATED_LATITUDE_LONGITUDE);
    addParameter(CF.GRID_NORTH_POLE_LATITUDE, northPoleLat);
    addParameter(CF.GRID_NORTH_POLE_LONGITUDE, northPoleLon);
  }

  private void buildRotationMatrices() {

    double betaRad = 0.;
    double gammaRad = 0.;

    if (northPole.getY() == 90.)
    {
        betaRad  = 0.;
        gammaRad = northPole.getX() * RAD_PER_DEG;
    }
    else
    {
        betaRad  = -(90. - northPole.getY()) * RAD_PER_DEG;
        gammaRad = (northPole.getX() + 180.) * RAD_PER_DEG;
    }

    double cosBeta = Math.cos(betaRad);
    double sinBeta = Math.sin(betaRad);

    double cosGamma = Math.cos(gammaRad);
    double sinGamma = Math.sin(gammaRad);

    rotY[0][0] = cosBeta;
    rotY[0][1] = 0.;
    rotY[0][2] = -sinBeta;
    rotY[1][0] = 0.;
    rotY[1][1] = 1.;
    rotY[1][2] = 0.;
    rotY[2][0] = sinBeta;
    rotY[2][1] = 0.;
    rotY[2][2] = cosBeta;

    rotZ[0][0] = cosGamma;
    rotZ[0][1] = sinGamma;
    rotZ[0][2] = 0.;
    rotZ[1][0] = -sinGamma;
    rotZ[1][1] = cosGamma;
    rotZ[1][2] = 0.;
    rotZ[2][0] = 0.;
    rotZ[2][1] = 0.;
    rotZ[2][2] = 1.;
  }

  public ProjectionPointImpl getNorthPole() {
    return new ProjectionPointImpl(northPole);
  }

  @Override
  public ProjectionImpl constructCopy() {
    ProjectionImpl result =  new RotatedPole(northPole.getY(), northPole.getX());
    result.setDefaultMapArea(defaultMapArea);
    result.setName(name);
    return result;
  }

  public String paramsToString() {
    return " northPole= " + northPole;
  }

  /**
   * Transform a "real" longitude and latitude into the rotated longitude (X) and
   * rotated latitude (Y).
   */
  public ProjectionPoint latLonToProj(LatLonPoint latlon, ProjectionPointImpl destPoint) {
    double lat = latlon.getLatitude();
    double lon = latlon.getLongitude();

    //	Lon-lat pair to xyz coordinates on sphere with radius 1
    double[] p0 = new double[]
            {Math.cos(lat * RAD_PER_DEG) * Math.cos(lon * RAD_PER_DEG),
                    Math.cos(lat * RAD_PER_DEG) * Math.sin(lon * RAD_PER_DEG),
                    Math.sin(lat * RAD_PER_DEG)};

    //	Rotate around Z-axis
    double[] p1 = new double[]
            {rotZ[0][0] * p0[0] + rotZ[0][1] * p0[1] + rotZ[0][2] * p0[2],
                    rotZ[1][0] * p0[0] + rotZ[1][1] * p0[1] + rotZ[1][2] * p0[2],
                    rotZ[2][0] * p0[0] + rotZ[2][1] * p0[1] + rotZ[2][2] * p0[2]};

    //	Rotate around Y-axis
    double[] p2 = new double[]
            {rotY[0][0] * p1[0] + rotY[0][1] * p1[1] + rotY[0][2] * p1[2],
                    rotY[1][0] * p1[0] + rotY[1][1] * p1[1] + rotY[1][2] * p1[2],
                    rotY[2][0] * p1[0] + rotY[2][1] * p1[1] + rotY[2][2] * p1[2]};

    final double lonR = LatLonPointImpl.range180( Math.atan2(p2[1], p2[0]) * DEG_PER_RAD);
    //final double lonR = Math.atan2(p2[1], p2[0]) * DEG_PER_RAD;
    final double latR = Math.asin(p2[2]) * DEG_PER_RAD;

    if (destPoint == null)
      destPoint =  new ProjectionPointImpl(lonR, latR);
    else
      destPoint.setLocation(lonR, latR);

     if (show)
        System.out.println("LatLon= " + latlon+" proj= " + destPoint);

    return destPoint;
  }

  /**
   * Transform a rotated longitude (X) and rotated latitude (Y) into a "real"
   * longitude-latitude pair.
   */
  public LatLonPoint projToLatLon(ProjectionPoint ppt, LatLonPointImpl destPoint) {
    //	"x" and "y" input for rotated pole coords are actually a lon-lat pair
    final double lonR = LatLonPointImpl.range180( ppt.getX()); // LOOK guessing
    final double latR = ppt.getY();

    //	Lon-lat pair to xyz coordinates on sphere with radius 1
    double[] p0 = new double[]
            {Math.cos(latR * RAD_PER_DEG) * Math.cos(lonR * RAD_PER_DEG),
                    Math.cos(latR * RAD_PER_DEG) * Math.sin(lonR * RAD_PER_DEG),
                    Math.sin(latR * RAD_PER_DEG)};

    //	Inverse rotate around Y-axis (using transpose of Y matrix)
    double[] p1 = new double[]
            {rotY[0][0] * p0[0] + rotY[1][0] * p0[1] + rotY[2][0] * p0[2],
                    rotY[0][1] * p0[0] + rotY[1][1] * p0[1] + rotY[2][1] * p0[2],
                    rotY[0][2] * p0[0] + rotY[1][2] * p0[1] + rotY[2][2] * p0[2]};

    //	Inverse rotate around Z-axis (using transpose of Z matrix)
    double[] p2 = new double[]
            {rotZ[0][0] * p1[0] + rotZ[1][0] * p1[1] + rotZ[2][0] * p1[2],
                    rotZ[0][1] * p1[0] + rotZ[1][1] * p1[1] + rotZ[2][1] * p1[2],
                    rotZ[0][2] * p1[0] + rotZ[1][2] * p1[1] + rotZ[2][2] * p1[2]};

    final double lon = Math.atan2(p2[1], p2[0]) * DEG_PER_RAD;
    final double lat = Math.asin(p2[2]) * DEG_PER_RAD;

    if (destPoint == null)
      destPoint = new LatLonPointImpl(lat, lon);
    else
      destPoint.set(lat, lon);

    if (show)
       System.out.println("Proj= " + ppt+" latlon= " + destPoint);

    return destPoint;
  }

  public boolean crossSeam(ProjectionPoint pt1, ProjectionPoint pt2) {
     return Math.abs(pt1.getX() - pt2.getX()) > 270.0;
  }

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

  @Override
  public int hashCode() {
    return northPole.hashCode();
  }

}