Re: [thredds] Curious effect from attributes semi_major/minor_axis on southern hemisphere images created by ncWMS

  • To: Egil Støren <egil.storen@xxxxxx>
  • Subject: Re: [thredds] Curious effect from attributes semi_major/minor_axis on southern hemisphere images created by ncWMS
  • From: Heiko Klein <Heiko.Klein@xxxxxx>
  • Date: Mon, 28 Jan 2013 19:55:17 +0100
Hi,

I moved this request to the netcdf-java list, since I can see the problems also in the latest toolsUI. In newer versions (4.3.14) it looks even worth with a different scale (km -> m):

in https://github.com/Unidata/thredds/blob/master/cdm/src/main/java/ucar/unidata/geoloc/projection/proj4/StereographicAzimuthalProjection.java



Please find attached a replacement for above file, fixing both south-pole x-flipping and scale.

Heiko






On 2013-01-25 16:28, Egil Støren wrote:
Dear all,

We have a netCDF file containing ice concentrations around antarctica
that is available at
http://thredds.met.no/thredds/catalog/egiltest/ES/data/catalog.html?dataset=egiltest/ES/data/iceconc.nc.
Using the Godiva2 viewer, it is obvious that the data is displayed
inverted along a vertical line centred on the map. This is most easily
seen using the south polar stereographic projection.

The netCDF file has a variable containing projection info:

     int Polar_Stereographic_Grid ;
         Polar_Stereographic_Grid:grid_mapping_name =
"polar_stereographic" ;
         Polar_Stereographic_Grid:straight_vertical_longitude_from_pole
= 0.f ;
         Polar_Stereographic_Grid:latitude_of_projection_origin = -90.f ;
         Polar_Stereographic_Grid:standard_parallel = -70.f ;
         Polar_Stereographic_Grid:false_easting = 0.f ;
         Polar_Stereographic_Grid:false_northing = 0.f ;
         Polar_Stereographic_Grid:semi_major_axis = 6378273.f ;
         Polar_Stereographic_Grid:semi_minor_axis = 6356890.f ;
         Polar_Stereographic_Grid:proj4_string = "+proj=stere +a=6378273
+b=6356889.44891 +lat_0=-90 +lat_ts=-70 +lon_0=0" ;

When I removed the following attributes from this variable, the data
suddenly displayed correctly in Godiva2:

         Polar_Stereographic_Grid:semi_major_axis = 6378273.f ;
         Polar_Stereographic_Grid:semi_minor_axis = 6356890.f ;

The modified version can be seen at
http://thredds.met.no/thredds/catalog/egiltest/ES/data/catalog.html?dataset=egiltest/ES/data/iceconc_1.nc.


Since these attributes only gives a more exact description of the shape
of the earth (compared to not having them), I think there must be a bug
somewhere. Have anybody experienced this behaviour? Is this a known bug?

We are using thredds version 4.2.9.

Best regards,

   Egil Støren
   Norwegian Meteorological Institute

--
Dr. Heiko Klein                              Tel. + 47 22 96 32 58
Development Section / IT Department          Fax. + 47 22 69 63 55
Norwegian Meteorological Institute           http://www.met.no
P.O. Box 43 Blindern  0313 Oslo NORWAY
/*
Copyright 2006 Jerry Huxtable

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

/*
 * This file was semi-automatically converted from the public-domain USGS PROJ 
source.
 */
package ucar.unidata.geoloc.projection.proj4;

import java.awt.geom.Point2D;
import java.util.Formatter;

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

/**
 * taken from the USGS PROJ package.
 *
 * @author Heiko.Klein@xxxxxx
 */
public class StereographicAzimuthalProjection extends ProjectionImpl {
  // projection parameters
  double projectionLatitude, projectionLongitude; // origin in radian
  double n; // Math.sin(projectionLatitude)
  double scaleFactor, trueScaleLatitude;  // scale or trueScale in radian
  double falseEasting, falseNorthing; // km

  // earth shape
  private Earth earth;
  private double e;   // earth.getEccentricity
  private double totalScale; // scale to convert cartesian coords in km

  private final static int NORTH_POLE = 1;
  private final static int SOUTH_POLE = 2;
  private final static int EQUATOR = 3;
  private final static int OBLIQUE = 4;

  private final static double TOL = 1.e-8;

  private double akm1, sinphi0, cosphi0;
  private int mode;

  public StereographicAzimuthalProjection() {  // polar strerographic with true 
longitude at 60 deg
    this(90.0, 0.0, 0.9330127018922193, 60., 0, 0, new Earth());
  }

  /**
   * Construct a Stereographic Projection.
   *
   * @param latt         tangent point of projection, also origin of
   *                     projection coord system, in degree
   * @param lont         tangent point of projection, also origin of
   *                     projection coord system, in degree
   * @param trueScaleLat latitude in degree where scale is scale
   * @param scale        scale factor at tangent point, "normally 1.0 but may 
be reduced"
   */
  public StereographicAzimuthalProjection(double latt, double lont, double 
scale, double trueScaleLat, double false_easting, double false_northing, Earth 
earth) {
    super("StereographicAzimuthalProjection", false);

    projectionLatitude = Math.toRadians(latt);
    n = Math.abs(Math.sin(projectionLatitude));
    projectionLongitude = Math.toRadians(lont);
    trueScaleLatitude = Math.toRadians(trueScaleLat);
    scaleFactor = Math.abs(scale);
    falseEasting = false_easting;
    falseNorthing = false_northing;

    // earth figure
    this.earth = earth;
    this.e = earth.getEccentricity();
    this.totalScale = earth.getMajor() * 0.001; // scale factor for cartesion 
coords in km.
    initialize();

    // parameters
    addParameter(CF.GRID_MAPPING_NAME, CF.STEREOGRAPHIC);
    addParameter(CF.LONGITUDE_OF_PROJECTION_ORIGIN, lont);
    addParameter(CF.LATITUDE_OF_PROJECTION_ORIGIN, latt);
    addParameter(CF.SCALE_FACTOR_AT_PROJECTION_ORIGIN, scale);
    if ((false_easting != 0.0) || (false_northing != 0.0)) {
      addParameter(CF.FALSE_EASTING, false_easting);
      addParameter(CF.FALSE_NORTHING, false_northing);
      addParameter(CDM.UNITS, "km");
    }
    addParameter(CF.SEMI_MAJOR_AXIS, earth.getMajor());
    addParameter(CF.INVERSE_FLATTENING, 1.0 / earth.getFlattening());

    //System.err.println(paramsToString());

  }

  private void initialize() {
    double t;

    if (Math.abs((t = Math.abs(projectionLatitude)) - MapMath.HALFPI) < 
MapMath.EPS10)
      mode = projectionLatitude < 0. ? SOUTH_POLE : NORTH_POLE;
    else
      mode = t > MapMath.EPS10 ? OBLIQUE : EQUATOR;
    trueScaleLatitude = Math.abs(trueScaleLatitude);
    if (earth.isSpherical()) { // sphere
      switch (mode) {
        case OBLIQUE:
          sinphi0 = Math.sin(projectionLatitude);
          cosphi0 = Math.cos(projectionLatitude);
        case EQUATOR:
          akm1 = 2. * scaleFactor;
          break;
        case SOUTH_POLE:
        case NORTH_POLE:
          akm1 = Math.abs(trueScaleLatitude - MapMath.HALFPI) >= MapMath.EPS10 ?
                  Math.cos(trueScaleLatitude) / Math.tan(MapMath.QUARTERPI - .5 
* trueScaleLatitude) :
                  2. * scaleFactor;
          break;
      }
    } else { // ellipsoid
      double X;

      switch (mode) {
        case NORTH_POLE:
        case SOUTH_POLE:
          if (Math.abs(trueScaleLatitude - MapMath.HALFPI) < MapMath.EPS10)
            akm1 = 2. * scaleFactor /
                    Math.sqrt(Math.pow(1 + e, 1 + e) * Math.pow(1 - e, 1 - e));
          else {
            akm1 = Math.cos(trueScaleLatitude) /
                    MapMath.tsfn(trueScaleLatitude, t = 
Math.sin(trueScaleLatitude), e);
            t *= e;
            akm1 /= Math.sqrt(1. - t * t);
          }
          break;
        case EQUATOR:
          akm1 = 2. * scaleFactor;
          break;
        case OBLIQUE:
          t = Math.sin(projectionLatitude);
          X = 2. * Math.atan(ssfn(projectionLatitude, t, e)) - MapMath.HALFPI;
          t *= e;
          akm1 = 2. * scaleFactor * Math.cos(projectionLatitude) / Math.sqrt(1. 
- t * t);
          sinphi0 = Math.sin(X);
          cosphi0 = Math.cos(X);
          break;
      }
    }
  }

  public Point2D.Double project(double lam, double phi, Point2D.Double xy) {
    double coslam = Math.cos(lam);
    double sinlam = Math.sin(lam);
    double sinphi = Math.sin(phi);

    if (earth.isSpherical()) { // sphere
      double cosphi = Math.cos(phi);

      switch (mode) {
        case EQUATOR:
          xy.y = 1. + cosphi * coslam;
          if (xy.y <= MapMath.EPS10)
            throw new RuntimeException("I");
          xy.x = (xy.y = akm1 / xy.y) * cosphi * sinlam;
          xy.y *= sinphi;
          break;
        case OBLIQUE:
          xy.y = 1. + sinphi0 * sinphi + cosphi0 * cosphi * coslam;
          if (xy.y <= MapMath.EPS10)
            throw new RuntimeException("I");
          xy.x = (xy.y = akm1 / xy.y) * cosphi * sinlam;
          xy.y *= cosphi0 * sinphi - sinphi0 * cosphi * coslam;
          break;
        case NORTH_POLE:
          coslam = -coslam;
          phi = -phi;
        case SOUTH_POLE:
          if (Math.abs(phi - MapMath.HALFPI) < TOL)
            throw new RuntimeException("I");
          xy.x = sinlam * (xy.y = akm1 * Math.tan(MapMath.QUARTERPI + .5 * 
phi));
          xy.y *= coslam;
          break;
      }
    } else { // ellipsoid
      double sinX = 0, cosX = 0, X, A;

      if (mode == OBLIQUE || mode == EQUATOR) {
        sinX = Math.sin(X = 2. * Math.atan(ssfn(phi, sinphi, e)) - 
MapMath.HALFPI);
        cosX = Math.cos(X);
      }
      switch (mode) {
        case OBLIQUE:
          A = akm1 / (cosphi0 * (1. + sinphi0 * sinX + cosphi0 * cosX * 
coslam));
          xy.y = A * (cosphi0 * sinX - sinphi0 * cosX * coslam);
          xy.x = A * cosX;
          break;
        case EQUATOR:
          A = 2. * akm1 / (1. + cosX * coslam);
          xy.y = A * sinX;
          xy.x = A * cosX;
          break;
        case SOUTH_POLE:
          phi = -phi;
          coslam = -coslam;
          sinphi = -sinphi;
        case NORTH_POLE:
          xy.x = akm1 * MapMath.tsfn(phi, sinphi, e);
          xy.y = -xy.x * coslam;
          break;
      }
      xy.x = xy.x * sinlam;
    }
    return xy;
  }

  public Point2D.Double projectInverse(double x, double y, Point2D.Double lp) {
    if (earth.isSpherical()) {
      double c, rh, sinc, cosc;

      sinc = Math.sin(c = 2. * Math.atan((rh = MapMath.distance(x, y)) / akm1));
      cosc = Math.cos(c);
      lp.x = 0.;
      switch (mode) {
        case EQUATOR:
          if (Math.abs(rh) <= MapMath.EPS10)
            lp.y = 0.;
          else
            lp.y = Math.asin(y * sinc / rh);
          if (cosc != 0. || x != 0.)
            lp.x = Math.atan2(x * sinc, cosc * rh);
          break;
        case OBLIQUE:
          if (Math.abs(rh) <= MapMath.EPS10)
            lp.y = projectionLatitude;
          else
            lp.y = Math.asin(cosc * sinphi0 + y * sinc * cosphi0 / rh);
          if ((c = cosc - sinphi0 * Math.sin(lp.y)) != 0. || x != 0.)
            lp.x = Math.atan2(x * sinc * cosphi0, c * rh);
          break;
        case NORTH_POLE:
          y = -y;
        case SOUTH_POLE:
          if (Math.abs(rh) <= MapMath.EPS10)
            lp.y = projectionLatitude;
          else
            lp.y = Math.asin(mode == SOUTH_POLE ? -cosc : cosc);
          lp.x = (x == 0. && y == 0.) ? 0. : Math.atan2(x, y);
          break;
      }
    } else {
      double cosphi, sinphi, tp, phi_l, rho, halfe, halfpi;

      rho = MapMath.distance(x, y);
      switch (mode) {
        case NORTH_POLE:
          y = -y;
        case SOUTH_POLE:
          phi_l = MapMath.HALFPI - 2. * Math.atan(tp = -rho / akm1);
          halfpi = -MapMath.HALFPI;
          halfe = -.5 * e;
          break;
        case OBLIQUE:
        case EQUATOR:
        default:
          cosphi = Math.cos(tp = 2. * Math.atan2(rho * cosphi0, akm1));
          sinphi = Math.sin(tp);
          phi_l = Math.asin(cosphi * sinphi0 + (y * sinphi * cosphi0 / rho));
          tp = Math.tan(.5 * (MapMath.HALFPI + phi_l));
          x *= sinphi;
          y = rho * cosphi0 * cosphi - y * sinphi0 * sinphi;
          halfpi = MapMath.HALFPI;
          halfe = .5 * e;
          break;
      }
      for (int i = 8; i-- != 0; phi_l = lp.y) {
        sinphi = e * Math.sin(phi_l);
        lp.y = 2. * Math.atan(tp * Math.pow((1. + sinphi) / (1. - sinphi), 
halfe)) - halfpi;
        if (Math.abs(phi_l - lp.y) < MapMath.EPS10) {
          if (mode == SOUTH_POLE)
            lp.y = -lp.y;
          lp.x = (x == 0. && y == 0.) ? 0. : Math.atan2(x, y);
          return lp;
        }
      }
      throw new RuntimeException("Iteration didn't converge");
    }
    return lp;
  }

  private double ssfn(double phit, double sinphi, double eccen) {
    sinphi *= eccen;
    return Math.tan(.5 * (MapMath.HALFPI + phit)) *
            Math.pow((1. - sinphi) / (1. + sinphi), .5 * eccen);
  }

  @Override
  public String getProjectionTypeLabel() {
    return "Stereographic Azimuthal Ellipsoidal Earth";
  }

  @Override
  public ProjectionImpl constructCopy() {
    ProjectionImpl result =  new 
StereographicAzimuthalProjection(Math.toDegrees(projectionLatitude), 
Math.toDegrees(projectionLongitude),
            scaleFactor, Math.toDegrees(trueScaleLatitude), falseEasting, 
falseNorthing, earth);
    result.setDefaultMapArea(defaultMapArea);
    return result;
  }

  @Override
  public String paramsToString() {
    Formatter f = new Formatter();
    f.format("origin lat,lon=%f,%f scale,trueScaleLat=%f,%f earth=%s", 
Math.toDegrees(projectionLatitude),
            Math.toDegrees(projectionLongitude), scaleFactor, 
Math.toDegrees(trueScaleLatitude), earth);
    return f.toString();
  }

  @Override
  public ProjectionPoint latLonToProj(LatLonPoint latLon, ProjectionPointImpl 
destPoint) {
    double fromLat = Math.toRadians(latLon.getLatitude());
    double theta = computeTheta(latLon.getLongitude());

    //System.err.println(Math.toDegrees(theta) + " " + Math.toDegrees(fromLat));
    Point2D.Double res = project(theta, fromLat, new Point2D.Double());

    destPoint.setLocation(totalScale * res.x + falseEasting, totalScale * res.y 
+ falseNorthing);
    return destPoint;
  }

  @Override
  public LatLonPoint projToLatLon(ProjectionPoint world, LatLonPointImpl 
result) {
    double fromX = (world.getX() - falseEasting) / totalScale; // assumes 
cartesian coords in km
    double fromY = (world.getY() - falseNorthing) / totalScale;

    Point2D.Double dst = projectInverse(fromX, fromY, new Point2D.Double());
    if (dst.x < -Math.PI)
      dst.x = -Math.PI;
    else if (dst.x > Math.PI)
      dst.x = Math.PI;
    if (projectionLongitude != 0)
      dst.x = MapMath.normalizeLongitude(dst.x + projectionLongitude);

    result.setLongitude(Math.toDegrees(dst.x));
    result.setLatitude(Math.toDegrees(dst.y));
    return result;
  }

  @Override
  public boolean crossSeam(ProjectionPoint pt1, ProjectionPoint pt2) {
    // TODO: not sure what this is, HK
    //       just taken from ucar.unidata.geoloc.projection.Stereographic
    return false;
  }

  @Override
  public boolean equals(Object proj) {
    if (!(proj instanceof StereographicAzimuthalProjection)) {
      return false;
    }
    StereographicAzimuthalProjection oo = (StereographicAzimuthalProjection) 
proj;
    if ((this.getDefaultMapArea() == null) != (oo.defaultMapArea == null)) 
return false; // common case is that these are null
    if (this.getDefaultMapArea() != null && 
!this.defaultMapArea.equals(oo.defaultMapArea)) return false;

    return ((this.projectionLatitude == oo.projectionLatitude)
          && (this.projectionLongitude == oo.projectionLongitude)
          && (this.scaleFactor == oo.scaleFactor)
          && (this.trueScaleLatitude == oo.trueScaleLatitude)
          && (this.falseEasting == oo.falseEasting)
          && (this.falseNorthing == oo.falseNorthing)
          && this.earth.equals(oo.earth));
  }

  @Override
  public int hashCode() {
    int hash = 3;
    hash = 67 * hash + (int) (Double.doubleToLongBits(this.projectionLatitude) 
^ (Double.doubleToLongBits(this.projectionLatitude) >>> 32));
    hash = 67 * hash + (int) (Double.doubleToLongBits(this.projectionLongitude) 
^ (Double.doubleToLongBits(this.projectionLongitude) >>> 32));
    hash = 67 * hash + (int) (Double.doubleToLongBits(this.scaleFactor) ^ 
(Double.doubleToLongBits(this.scaleFactor) >>> 32));
    hash = 67 * hash + (int) (Double.doubleToLongBits(this.trueScaleLatitude) ^ 
(Double.doubleToLongBits(this.trueScaleLatitude) >>> 32));
    hash = 67 * hash + (int) (Double.doubleToLongBits(this.falseEasting) ^ 
(Double.doubleToLongBits(this.falseEasting) >>> 32));
    hash = 67 * hash + (int) (Double.doubleToLongBits(this.falseNorthing) ^ 
(Double.doubleToLongBits(this.falseNorthing) >>> 32));
    hash = 67 * hash + (this.earth != null ? this.earth.hashCode() : 0);
    return hash;
  }

  private double computeTheta(double lon) {
    double dlon = LatLonPointImpl.lonNormal(lon - 
Math.toDegrees(projectionLongitude));
    return n * Math.toRadians(dlon);
  }

  static private void test(ProjectionImpl proj, double[] lat, double[] lon) {
    double[] x = new double[lat.length];
    double[] y = new double[lat.length];
    for (int i = 0; i < lat.length; ++i) {
      LatLonPoint lp = new LatLonPointImpl(lat[i], lon[i]);
      ProjectionPointImpl p = (ProjectionPointImpl) proj.latLonToProj(lp, new 
ProjectionPointImpl());
      x[i] = p.x;
      y[i] = p.y;
    }
    for (int i = 0; i < lat.length; ++i) {
      ProjectionPointImpl p = new ProjectionPointImpl(x[i], y[i]);
      LatLonPointImpl lp = (LatLonPointImpl) proj.projToLatLon(p);
      if ((Math.abs(lp.getLatitude() - lat[i]) > 1e-5)
              || (Math.abs(lp.getLongitude() - lon[i]) > 1e-5)) {
        if (Math.abs(lp.getLatitude()) > 89.99 &&
                (Math.abs(lp.getLatitude() - lat[i]) < 1e-5)) {
          // ignore longitude singularities at poles
        } else {
          System.err.print("ERROR:");
        }
      }
      System.out.println("reverse:" + p.x + ", " + p.y + ": " + 
lp.getLatitude() + ", " + lp.getLongitude());

    }

  }

  static public void main(String[] args) {
    // test-code
    Earth e = new Earth(6378137., 0, 298.257224);
    StereographicAzimuthalProjection proj = new 
StereographicAzimuthalProjection(90., 0., 0.93306907, 90., 0., 0., e);

    double[] lat = {60., 90., 60.};
    double[] lon = {0., 0., 10.};
    test(proj, lat, lon);

    proj = new StereographicAzimuthalProjection(90., -45., 0.96985819, 90., 0., 
0., e);
    test(proj, lat, lon);

    // southpole
    proj = new StereographicAzimuthalProjection(-90., 0., -1, -70., 0., 0., e);

    double[] latS = {-60., -90., -60.};
    double[] lonS = {0., 0., 10.};
    test(proj, latS, lonS);


  }
}

  • 2013 messages navigation, sorted by:
    1. Thread
    2. Subject
    3. Author
    4. Date
    5. ↑ Table Of Contents
  • Search the thredds archives: