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);
}
}
