/*
* Copyright (c) 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.
*/
package ucar.unidata.geoloc.projection.proj4;
import ucar.nc2.constants.CDM;
import ucar.nc2.constants.CF;
import ucar.unidata.geoloc.*;
import ucar.unidata.util.Parameter;
import java.util.Formatter;
/**
* Adapted from com.jhlabs.map.proj.LambertConformalConicProjection
*
* @see "http://www.jhlabs.com/java/maps/proj/index.html"
* @see "http://trac.osgeo.org/proj/"
*
* @author caron
* @since Oct 20, 2009
*/
public class LambertConformalConicEllipse extends ProjectionImpl {
private final static double TOL = 1.0e-10;
// projection parameters
private double lat0deg, lon0deg; // projection origin, degrees
private double lat0rad, lon0rad; // projection origin, radians
private double par1deg, par2deg; // standard parellels, degrees
private double par1rad, par2rad; // standard parellels, radians
private double falseEasting, falseNorthing; // km
// earth shape
private Earth earth;
private double e; // earth.getEccentricitySquared
private double es; // earth.getEccentricitySquared
private double totalScale; // scale to convert cartesian coords in km
private double n;
private double c;
private double rho0;
// spherical vs ellipsoidal
private boolean isSpherical;
@Override
public ProjectionImpl constructCopy() {
ProjectionImpl result = new LambertConformalConicEllipse(getOriginLat(), getOriginLon(), getParallelOne(), getParallelTwo(),
getFalseEasting(), getFalseNorthing(), getEarth());
result.setDefaultMapArea(defaultMapArea);
result.setName(name);
return result;
}
/**
* Constructor with default parameters
*/
public LambertConformalConicEllipse() {
this(23.0, -96.0, 29.5, 45.5, 0, 0, new Earth(6378137.0, 0.0, 298.257222101));
}
/**
* Construct a LambertConformal Projection, two standard parellels.
* For the one standard parellel case, set them both to the same value.
*
* @param lat0 lat origin of the coord. system on the projection plane
* @param lon0 lon origin of the coord. system on the projection plane
* @param par1 standard parallel 1
* @param par2 standard parallel 2
* @param falseEasting natural_x_coordinate + false_easting = x coordinate in km
* @param falseNorthing natural_y_coordinate + false_northing = y coordinate in km
* @param earth shape of the earth
* @throws IllegalArgumentException if lat0, par1, par2 = +/-90 deg
*/
public LambertConformalConicEllipse(double lat0, double lon0, double par1, double par2, double falseEasting, double falseNorthing, Earth earth) {
super("LambertConformalConicEllipse", false);
this.lat0deg = lat0;
this.lon0deg = lon0;
this.lat0rad = Math.toRadians(lat0);
this.lon0rad = Math.toRadians(lat0);
this.par1deg = par1;
this.par2deg = par2;
this.par1rad = Math.toRadians(par1);
this.par2rad = Math.toRadians(par2);
this.falseEasting = falseEasting;
this.falseNorthing = falseNorthing;
this.earth = earth;
this.e = earth.getEccentricity();
this.es = earth.getEccentricitySquared();
this.isSpherical = (e == 0.0);
this.totalScale = earth.getMajor() * .001; // scale factor for cartesion coords in km.
initialize();
addParameter(CF.GRID_MAPPING_NAME, CF.LAMBERT_CONFORMAL_CONIC);
addParameter(CF.LATITUDE_OF_PROJECTION_ORIGIN, lat0);
addParameter(CF.LONGITUDE_OF_CENTRAL_MERIDIAN, lon0);
if (par2 == par1) {
addParameter(CF.STANDARD_PARALLEL, par1);
} else {
double[] data = new double[2];
data[0] = par1;
data[1] = par2;
addParameter(new Parameter(CF.STANDARD_PARALLEL, data));
}
if ((falseEasting != 0.0) || (falseNorthing != 0.0)) {
addParameter(CF.FALSE_EASTING, falseEasting);
addParameter(CF.FALSE_NORTHING, falseNorthing);
addParameter(CDM.UNITS, "km");
}
addParameter(CF.SEMI_MAJOR_AXIS, earth.getMajor());
addParameter(CF.INVERSE_FLATTENING, 1.0/earth.getFlattening());
}
private void initialize() {
if ( par1rad == 0 )
par1rad = par2rad = lat0rad;
if (Math.abs(par1rad + par2rad) < TOL)
throw new IllegalArgumentException("par1rad + par2rad < TOL");
double sinphi = Math.sin(par1rad);
double cosphi = Math.cos(par1rad);
boolean isSecant = Math.abs(par1rad - par2rad) >= TOL;
n = sinphi;
if (!isSpherical) {
double ml1, m1;
m1 = MapMath.msfn(sinphi, cosphi, es);
ml1 = MapMath.tsfn(par1rad, sinphi, e);
if (isSecant) {
n = Math.log(m1 /
MapMath.msfn(sinphi = Math.sin(par2rad), Math.cos(par2rad), es));
n /= Math.log(ml1 / MapMath.tsfn(par2rad, sinphi, e));
}
c = (rho0 = m1 * Math.pow(ml1, -n) / n);
rho0 *= (Math.abs(Math.abs(lat0rad) - MapMath.HALFPI) < TOL) ? 0. :
Math.pow(MapMath.tsfn(lat0rad, Math.sin(lat0rad), e), n);
} else {
if (isSecant)
n = Math.log(cosphi / Math.cos(par2rad)) /
Math.log(Math.tan(MapMath.QUARTERPI + .5 * par2rad) /
Math.tan(MapMath.QUARTERPI + .5 * par1rad));
c = cosphi * Math.pow(Math.tan(MapMath.QUARTERPI + .5 * par1rad), n) / n;
rho0 = (Math.abs(Math.abs(lat0rad) - MapMath.HALFPI) < TOL) ? 0. :
c * Math.pow(Math.tan(MapMath.QUARTERPI + .5 * lat0rad), -n);
}
}
/* private void precalculate() {
if (Math.abs(lat0 - PI_OVER_2) < TOLERANCE) {
throw new IllegalArgumentException("LambertConformal lat0 = 90");
}
if (Math.abs(lat0 + PI_OVER_2) < TOLERANCE) {
throw new IllegalArgumentException("LambertConformal lat0 = -90");
}
if (Math.abs(par1 - 90.0) < TOLERANCE) {
throw new IllegalArgumentException("LambertConformal par1 = 90");
}
if (Math.abs(par1 + 90.0) < TOLERANCE) {
throw new IllegalArgumentException("LambertConformal par1 = -90");
}
if (Math.abs(par2 - 90.0) < TOLERANCE) {
throw new IllegalArgumentException("LambertConformal par2 = 90");
}
if (Math.abs(par2 + 90.0) < TOLERANCE) {
throw new IllegalArgumentException("LambertConformal par2 = -90");
}
double par1r = Math.toRadians(this.par1);
double par2r = Math.toRadians(this.par2);
double t1 = Math.tan(Math.PI / 4 + par1r / 2);
double t2 = Math.tan(Math.PI / 4 + par2r / 2);
if (Math.abs(par2 - par1) < TOLERANCE) { // single parallel
n = Math.sin(par1r);
} else {
n = Math.log(Math.cos(par1r) / Math.cos(par2r))
/ Math.log(t2 / t1);
}
double t1n = Math.pow(t1, n);
F = Math.cos(par1r) * t1n / n;
earthRadiusTimesF = EARTH_RADIUS * F;
double t0n = Math.pow(Math.tan(Math.PI / 4 + lat0 / 2), n);
rho = EARTH_RADIUS * F / t0n;
lon0Degrees = Math.toDegrees(lon0);
// need to know the pole value for crossSeam
//Point2D pt = latLonToProj( 90.0, 0.0);
//maxY = pt.getY();
//System.out.println("LC = " +pt);
} */
// public LambertConformalConicEllipse(double lat0, double lon0, double par1, double par2, double falseEasting, double falseNorthing, Earth earth) {
@Override
public boolean equals(Object o) {
if (this == o) return true;
if (o == null || getClass() != o.getClass()) return false;
LambertConformalConicEllipse that = (LambertConformalConicEllipse) o;
if (Double.compare(that.falseEasting, falseEasting) != 0) return false;
if (Double.compare(that.falseNorthing, falseNorthing) != 0) return false;
if (Double.compare(that.lat0rad, lat0rad) != 0) return false;
if (Double.compare(that.lon0rad, lon0rad) != 0) return false;
if (Double.compare(that.par1rad, par1rad) != 0) return false;
if (Double.compare(that.par2rad, par2rad) != 0) return false;
if (!earth.equals(that.earth)) return false;
return true;
}
@Override
public int hashCode() {
int result;
long temp;
temp = Double.doubleToLongBits(lat0rad);
result = (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(lon0rad);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(par1rad);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(par2rad);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(falseEasting);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(falseNorthing);
result = 31 * result + (int) (temp ^ (temp >>> 32));
result = 31 * result + earth.hashCode();
return result;
}
/*
* Check for equality with the Object in question
*
* @param proj object to check
* @return true if they are equal
*
public boolean equals2(Object proj) {
if (!(proj instanceof LambertConformalConicEllipse))
return false;
LambertConformalConicEllipse oo = (LambertConformalConicEllipse) 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.getParallelOne() == oo.getParallelOne())
&& (this.getParallelTwo() == oo.getParallelTwo())
&& (this.getOriginLat() == oo.getOriginLat())
&& (this.getOriginLon() == oo.getOriginLon())
&& this.earth.equals(oo.earth));
} */
// bean properties
public Earth getEarth() { return earth; }
/**
* Get the second standard parallel
*
* @return the second standard parallel in degrees
*/
public double getParallelTwo() {
return par2deg;
}
/**
* Get the first standard parallel
*
* @return the first standard parallel in degrees
*/
public double getParallelOne() {
return par1deg;
}
/**
* Get the origin longitude.
*
* @return the origin longitude in degrees
*/
public double getOriginLon() {
return lon0deg;
}
/**
* Get the origin latitude.
*
* @return the origin latitude in degrees
*/
public double getOriginLat() {
return lat0deg;
}
/**
* Get the false easting, in km.
*
* @return the false easting in km
*/
public double getFalseEasting() {
return falseEasting;
}
/**
* Get the false northing, in km.
*
* @return the false northing in km
*/
public double getFalseNorthing() {
return falseNorthing;
}
/**
* Get the label to be used in the gui for this type of projection
*
* @return Type label
*/
public String getProjectionTypeLabel() {
return "Lambert Conformal Conic Ellipsoidal Earth";
}
/**
* Create a String of the parameters.
*
* @return a String of the parameters
*/
public String paramsToString() {
Formatter f = new Formatter();
f.format("origin lat,lon=%f,%f parellels=%f,%f earth=%s", lat0deg, lon0deg, par1deg, par2deg, earth );
return f.toString();
}
/*
* Create a WKS string LOOK NOT DONE
*
* @return WKS string
*
public String toWKS() {
Formatter sbuff = new Formatter();
sbuff.format("PROJCS[\"%s\",", getName());
if (isSpherical) {
sbuff.format("GEOGCS[\"Normal Sphere (r=6371007)\",");
sbuff.format("DATUM[\"unknown\",");
sbuff.format("SPHEROID[\"sphere\",6371007,0]],");
} else {
sbuff.format("GEOGCS[\"WGS 84\",");
sbuff.format("DATUM[\"WGS_1984\",");
sbuff.format("SPHEROID[\"WGS 84\",6378137,298.257223563],");
sbuff.format("TOWGS84[0,0,0,0,0,0,0]],");
}
sbuff.format("PRIMEM[\"Greenwich\",0],");
sbuff.format("UNIT[\"degree\",0.0174532925199433]],");
sbuff.format("PROJECTION[\"Lambert_Conformal_Conic_1SP\"],");
sbuff.format("PARAMETER[\"latitude_of_origin\",").append(getOriginLat()).append("],"); // LOOK assumes getOriginLat = getParellel
sbuff.format("PARAMETER[\"central_meridian\",").append(getOriginLon()).append("],");
sbuff.format("PARAMETER[\"scale_factor\",1],");
sbuff.format("PARAMETER[\"false_easting\",").append(falseEasting).append("],");
sbuff.format("PARAMETER[\"false_northing\",").append(falseNorthing).append("],");
return sbuff.toString();
} */
/**
* This returns true when the line between pt1 and pt2 crosses the seam.
* When the cone is flattened, the "seam" is lon0 +- 180.
*
* @param pt1 point 1
* @param pt2 point 2
* @return true when the line between pt1 and pt2 crosses the seam.
*/
public boolean crossSeam(ProjectionPoint pt1, ProjectionPoint pt2) {
if (ProjectionPointImpl.isInfinite(pt1) || ProjectionPointImpl.isInfinite(pt2)) {
return true;
}
return false;
/* opposite signed X values, larger then 5000 km LOOK ????
return (pt1.getX() * pt2.getX() < 0)
&& (Math.abs(pt1.getX() - pt2.getX()) > 5000.0); */
}
/*
public ProjectionPoint latLonToProj(LatLonPoint latLon,
ProjectionPointImpl result) {
double toX, toY;
double fromLat = latLon.getLatitude();
double fromLon = latLon.getLongitude();
fromLat = Math.toRadians(fromLat);
double dlon = LatLonPointImpl.lonNormal(fromLon - lon0Degrees);
double theta = n * Math.toRadians(dlon);
double tn = Math.pow(Math.tan(PI_OVER_4 + fromLat / 2), n);
double r = earthRadiusTimesF / tn;
toX = r * Math.sin(theta);
toY = rho - r * Math.cos(theta);
result.setLocation(toX + falseEasting, toY + falseNorthing);
return result;
}
/*
public Point2D.Double project(double x, double y, Point2D.Double out) {
double rho;
if (Math.abs(Math.abs(y) - MapMath.HALFPI) < 1e-10)
rho = 0.0;
else
rho = c * (spherical ?
Math.pow(Math.tan(MapMath.QUARTERPI + .5 * y), -n) :
Math.pow(MapMath.tsfn(y, Math.sin(y), e), n));
out.x = scaleFactor * (rho * Math.sin(x *= n));
out.y = scaleFactor * (rho0 - rho * Math.cos(x));
return out;
}
*/
private double computeTheta(double lon) {
double dlon = LatLonPointImpl.lonNormal(lon - lon0deg);
return n * Math.toRadians(dlon);
}
// also see Snyder p 101
public ProjectionPoint latLonToProj(LatLonPoint latLon, ProjectionPointImpl result) {
double fromLat = Math.toRadians(latLon.getLatitude());
double theta = computeTheta(latLon.getLongitude());
double rho = 0.0;
if (Math.abs(Math.abs(fromLat) - MapMath.HALFPI) >= TOL) {
double term;
if (isSpherical)
term = Math.pow(Math.tan(MapMath.QUARTERPI + .5 * fromLat), -n);
else
term = Math.pow(MapMath.tsfn(fromLat, Math.sin(fromLat), e), n);
rho = c * term;
}
double toX = (rho * Math.sin(theta));
double toY = (rho0 - rho * Math.cos(theta));
result.setLocation(totalScale * toX + falseEasting, totalScale * toY + falseNorthing);
return result;
}
/*
public LatLonPoint projToLatLon(ProjectionPoint world,
LatLonPointImpl result) {
double toLat, toLon;
double fromX = world.getX() - falseEasting;
double fromY = world.getY() - falseNorthing;
double rhop = rho;
if (n < 0) {
rhop *= -1.0;
fromX *= -1.0;
fromY *= -1.0;
}
double yd = (rhop - fromY);
double theta = Math.atan2(fromX, yd);
double r = Math.sqrt(fromX * fromX + yd * yd);
if (n < 0.0) {
r *= -1.0;
}
toLon = (Math.toDegrees(theta / n + lon0));
if (Math.abs(r) < TOLERANCE) {
toLat = ((n < 0.0)
? -90.0
: 90.0);
} else {
double rn = Math.pow(EARTH_RADIUS * F / r, 1 / n);
toLat = Math.toDegrees(2.0 * Math.atan(rn) - Math.PI / 2);
}
result.setLatitude(toLat);
result.setLongitude(toLon);
return result;
}
public Point2D.Double projectInverse(double x, double y, Point2D.Double out) {
x /= scaleFactor;
y /= scaleFactor;
double rho = MapMath.distance(x, y = rho0 - y);
if (rho != 0) {
if (n < 0.0) {
rho = -rho;
x = -x;
y = -y;
}
if (spherical)
out.y = 2.0 * Math.atan(Math.pow(c / rho, 1.0/n)) - MapMath.HALFPI;
else
out.y = MapMath.phi2(Math.pow(rho / c, 1.0/n), e);
out.x = Math.atan2(x, y) / n;
} else {
out.x = 0.0;
out.y = n > 0.0 ? MapMath.HALFPI : -MapMath.HALFPI;
}
return out;
}
} */
public LatLonPoint projToLatLon(ProjectionPoint world, LatLonPointImpl result) {
double toLat, toLon;
double fromX = (world.getX() - falseEasting) / totalScale; // assumes cartesion coords in km
double fromY = (world.getY() - falseNorthing) / totalScale;
fromY = rho0 - fromY;
double rho = MapMath.distance(fromX, fromY);
if (rho != 0) {
if (n < 0.0) {
rho = -rho;
fromX = -fromX;
fromY = -fromY;
}
if (isSpherical)
toLat = 2.0 * Math.atan(Math.pow(c / rho, 1.0/n)) - MapMath.HALFPI;
else
toLat = MapMath.phi2(Math.pow(rho / c, 1.0/n), e);
toLon = Math.atan2(fromX, fromY) / n;
//coverity[swapped_arguments]
} else {
toLon = 0.0;
toLat = n > 0.0 ? MapMath.HALFPI : -MapMath.HALFPI;
}
result.setLatitude(Math.toDegrees(toLat));
result.setLongitude(Math.toDegrees(toLon) + lon0deg);
return result;
}
////////////////////////////////////////////////////////
// test
/*
proj +inv +proj=aea +lat_0=23.0 +lat_1=29.5 +lat_2=45.5 +a=6378137.0 +rf=298.257222101 +b=6356752.31414 +lon_0=-96.0
Input X,Y: 1730692.593817677 1970917.991173046
results in:
Output Lon,Lat: -75.649278 39.089117
*
*/
private static void toProj(ProjectionImpl p, double lat, double lon) {
System.out.printf("lon,lat = %f %f%n", lon, lat);
ProjectionPoint pt = p.latLonToProj(lat, lon);
System.out.printf("x,y = %f %f%n", pt.getX(), pt.getY());
LatLonPoint ll = p.projToLatLon(pt);
System.out.printf("lon,lat = %f %f%n%n", ll.getLongitude(), ll.getLatitude());
}
private static void fromProj(ProjectionImpl p, double x, double y) {
System.out.printf("x,y = %f %f%n", x,y);
LatLonPoint ll = p.projToLatLon(x, y);
System.out.printf("lon,lat = %f %f%n", ll.getLongitude(), ll.getLatitude());
ProjectionPoint pt = p.latLonToProj(ll);
System.out.printf("x,y = %f %f%n%n", pt.getX(), pt.getY());
}
public static void main(String[] args) {
LambertConformalConicEllipse a = new LambertConformalConicEllipse(23.0, -96.0, 29.5, 45.5, 0, 0, new Earth(6378137.0, 0.0, 298.257222101));
System.out.printf("proj = %s %s%n%n", a.getName(), a.paramsToString());
//fromProj(a, 1730.692593817677, 1970.917991173046);
//toProj(a, 39.089117, -75.649278);
fromProj(a, 5747, 13470);
}
}