ucar.unidata.geoloc.projection.proj4.PolyconicProjection Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of cdm Show documentation
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 (c) 1998 - 2012. 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.
*/
/**
* Polyconic Projection.
* This file was semi-automatically converted from the public-domain USGS PROJ source.
*
* Bernhard Jenny, 19 September 2010: fixed spherical inverse.
* @author [email protected] 1/8/2012
*/
package ucar.unidata.geoloc.projection.proj4;
import java.util.Formatter;
import ucar.nc2.constants.CDM;
import ucar.nc2.constants.CF;
import ucar.unidata.geoloc.Earth;
import ucar.unidata.geoloc.LatLonPoint;
import ucar.unidata.geoloc.LatLonPointImpl;
import ucar.unidata.geoloc.ProjectionImpl;
import ucar.unidata.geoloc.ProjectionPoint;
import ucar.unidata.geoloc.ProjectionPointImpl;
public class PolyconicProjection extends ProjectionImpl {
private double ml0;
private double[] en;
private final static double TOL = 1e-10;
private final static double CONV = 1e-10;
private final static int N_ITER = 10;
private final static int I_ITER = 20;
private final static double ITOL = 1.e-12;
//New variables added
private boolean spherical = true;
private double projectionLatitude;
private double projectionLongitude;
private Earth ellipsoid;
private double es;
private double falseEasting;
private double falseNorthing;
private double totalScale;
public PolyconicProjection() {
this(23.56, 76.54);
}
public PolyconicProjection(double lat0, double lon0) {
this(lat0, lon0, new Earth());
}
public PolyconicProjection(double lat0, double lon0, Earth ellipsoid) {
this(lat0, lon0, 0.0, 0.0, ellipsoid);
}
public PolyconicProjection(double lat0, double lon0, double falseEasting, double falseNorthing, Earth ellipsoid) {
super("Polyconic", false);
//Initialization
this.projectionLatitude = Math.toRadians(lat0);
this.projectionLongitude = Math.toRadians(lon0);
this.ellipsoid = ellipsoid;
this.falseEasting = falseEasting;
this.falseNorthing = falseNorthing;
this.es = this.ellipsoid.getEccentricitySquared();
this.totalScale = this.ellipsoid.getMajor() * .001;
initialize();
//Adding parameters
addParameter(CF.GRID_MAPPING_NAME, name);
addParameter(CF.LATITUDE_OF_PROJECTION_ORIGIN, lat0);
addParameter(CF.LONGITUDE_OF_CENTRAL_MERIDIAN, lon0);
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, ellipsoid.getMajor());
addParameter(CF.INVERSE_FLATTENING, 1.0 / ellipsoid.getFlattening());
}
private ProjectionPoint project(double lplam, double lpphi, ProjectionPointImpl out) {
if (spherical) {
double cot, E;
if (Math.abs(lpphi) <= TOL) {
out.setLocation(lplam, ml0);
} else {
cot = 1. / Math.tan(lpphi);
double x = Math.sin(E = lplam * Math.sin(lpphi)) * cot;
double y = lpphi - projectionLatitude + cot * (1. - Math.cos(E));
out.setLocation(x, y);
}
} else {
/*
FORWARD(e_forward); ellipsoid
double ms, sp, cp;
if (fabs(lp.phi) <= TOL) {
xy.x = lp.lam; xy.y = -P->ml0; }
else {
sp = sin(lp.phi);
ms = fabs(cp = cos(lp.phi)) > TOL ? pj_msfn(sp, cp, P->es) / sp : 0.;
xy.x = ms * sin(lp.lam *= sp); // LOOK
xy.y = (pj_mlfn(lp.phi, sp, cp, P->en) - P->ml0) + ms * (1. - cos(lp.lam));
}
return (xy);
}
*/
double ms, sp, cp;
if (Math.abs(lpphi) <= TOL) {
out.setLocation(lplam, -ml0);
} else {
sp = Math.sin(lpphi);
ms = Math.abs(cp = Math.cos(lpphi)) > TOL ? MapMath.msfn(sp, cp, es) / sp : 0.;
lplam *= sp; // LOOK
double x = ms * Math.sin(lplam);
double y = (MapMath.mlfn(lpphi, sp, cp, en) - ml0) + ms * (1. - Math.cos(lplam));
out.setLocation(x, y);
}
}
return out;
}
private ProjectionPoint projectInverse(double xyx, double xyy, ProjectionPointImpl out) {
double lpphi;
if (spherical) {
double B, dphi, tp;
int i;
if (Math.abs(xyy = projectionLatitude + xyy) <= TOL) {
out.setLocation(xyx, 0);
} else {
lpphi = xyy;
B = xyx * xyx + xyy * xyy;
i = N_ITER;
do {
tp = Math.tan(lpphi);
lpphi -= (dphi = (xyy * (lpphi * tp + 1.) - lpphi
- .5 * (lpphi * lpphi + B) * tp)
/ ((lpphi - xyy) / tp - 1.));
} while (Math.abs(dphi) > CONV && --i > 0);
if (i == 0) {
out.setLocation(Double.NaN, Double.NaN);
}
double x = Math.asin(xyx * Math.tan(lpphi)) / Math.sin(lpphi);
double y = lpphi;
out.setLocation(x, y);
}
} else {
xyy += ml0;
if (Math.abs(xyy) <= TOL) {
out.setLocation(xyx, 0);
} else {
double r, c, sp, cp, s2ph, ml, mlb, mlp, dPhi;
int i;
r = xyy * xyy + xyx * xyx;
for (lpphi = xyy, i = I_ITER; i > 0; --i) {
sp = Math.sin(lpphi);
s2ph = sp * (cp = Math.cos(lpphi));
if (Math.abs(cp) < ITOL) {
throw new RuntimeException("I");
}
c = sp * (mlp = Math.sqrt(1. - es * sp * sp)) / cp;
ml = MapMath.mlfn(lpphi, sp, cp, en);
mlb = ml * ml + r;
mlp = (1.0 / es) / (mlp * mlp * mlp);
lpphi += (dPhi =
(ml + ml + c * mlb - 2. * xyy * (c * ml + 1.)) / (es * s2ph * (mlb - 2. * xyy * ml) / c
+ 2. * (xyy - ml) * (c * mlp - 1. / s2ph) - mlp - mlp));
if (Math.abs(dPhi) <= ITOL) {
break;
}
}
if (i == 0) {
out.setLocation(Double.NaN, Double.NaN);
}
c = Math.sin(lpphi);
double x = Math.asin(xyx * Math.tan(lpphi) * Math.sqrt(1. - es * c * c)) / Math.sin(lpphi);
double y = lpphi;
out.setLocation(x, y);
}
}
return out;
}
private void initialize() {
if (!spherical) {
en = MapMath.enfn(es);
if (en == null) {
throw new RuntimeException("E");
}
ml0 = MapMath.mlfn(projectionLatitude, Math.sin(projectionLatitude), Math.cos(projectionLatitude), en);
} else {
ml0 = -projectionLatitude;
}
}
// public PolyconicProjection(double lat0, double lon0, double falseEasting, double falseNorthing, Earth ellipsoid) {
@Override
public boolean equals(Object o) {
if (this == o) return true;
if (o == null || getClass() != o.getClass()) return false;
PolyconicProjection that = (PolyconicProjection) o;
if (Double.compare(that.falseEasting, falseEasting) != 0) return false;
if (Double.compare(that.falseNorthing, falseNorthing) != 0) return false;
if (Double.compare(that.projectionLatitude, projectionLatitude) != 0) return false;
if (Double.compare(that.projectionLongitude, projectionLongitude) != 0) return false;
if (!ellipsoid.equals(that.ellipsoid)) return false;
return true;
}
@Override
public int hashCode() {
int result;
long temp;
temp = Double.doubleToLongBits(projectionLatitude);
result = (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(projectionLongitude);
result = 31 * result + (int) (temp ^ (temp >>> 32));
result = 31 * result + ellipsoid.hashCode();
temp = Double.doubleToLongBits(falseEasting);
result = 31 * result + (int) (temp ^ (temp >>> 32));
temp = Double.doubleToLongBits(falseNorthing);
result = 31 * result + (int) (temp ^ (temp >>> 32));
return result;
}
/*
* Check for equality with the Object in question
*
* @param proj object to check
* @return true if they are equal
*
public boolean equals(Object proj) {
if (!(proj instanceof PolyconicProjection))
return false;
PolyconicProjection oo = (PolyconicProjection) 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.getOriginLatitude() == oo.getOriginLatitude())
&& (this.getOriginLongitude() == oo.getOriginLongitude())
&& this.ellipsoid.equals(oo.getEarth()));
} */
//Beans properties
public Earth getEarth() {
return ellipsoid;
}
/**
* Set the origin latitude.
*
* @param lat the origin latitude.
*/
public void setOriginLatitude(double lat) { //suppply in degrees
projectionLatitude = Math.toRadians(lat);
}
/**
* Get the origin longitude.
*
* @return the origin longitude in degrees
*/
public double getOriginLatitude() {
return Math.toDegrees(projectionLatitude);
}
/**
* Set the origin longitude.
*
* @param lon the origin longitude.
*/
public void setOriginLongitude(double lon) {//suppply in degrees
projectionLongitude = Math.toRadians(lon);
}
/**
* Get the origin longitude.
*
* @return the origin longitude in degrees
*/
public double getOriginLongitude() {
return Math.toDegrees(projectionLongitude);
}
/**
* Get the false easting, in km.
*
* @return the false easting in km
*/
public double getFalseEasting() {
return falseEasting;
}
/**
* Set the false_easting, in km.
* natural_x_coordinate + false_easting = x coordinate
*
* @param falseEasting x offset
*/
public void setFalseEasting(double falseEasting) {
this.falseEasting = falseEasting;
}
/**
* Get the false northing, in km.
*
* @return the false northing in km
*/
public double getFalseNorthing() {
return falseNorthing;
}
/**
* Set the false northing, in km.
* natural_y_coordinate + false_northing = y coordinate
*
* @param falseNorthing y offset
*/
public void setFalseNorthing(double falseNorthing) {
this.falseNorthing = falseNorthing;
}
/**
* Get the label to be used in the gui for this type of projection
*
* @return Type label
*/
public String getProjectionTypeLabel() {
return "Polyconic Projection";
}
/**
* Create a String of the parameters.
*
* @return a String of the parameters
*/
@Override
public String paramsToString() {
Formatter f = new Formatter();
f.format("origin lat=%f, origin lon=%f earth=%s", Math.toDegrees(projectionLatitude), Math.toDegrees(projectionLongitude), ellipsoid);
return f.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.
*/
@Override
public boolean crossSeam(ProjectionPoint pt1, ProjectionPoint pt2) {
if (ProjectionPointImpl.isInfinite(pt1) || ProjectionPointImpl.isInfinite(pt2)) {
return true;
}
// opposite signed X values, larger then 20,000 kml... similar to LambertConformal Conic
return (pt1.getX() * pt2.getX() < 0)
&& (Math.abs(pt1.getX() - pt2.getX()) > 20000.0);
}
/**
* Convert a LatLonPoint to projection coordinates
*
* @param latlon convert from these lat, lon coordinates
* @param result the object to write to
* @return the given result
*/
@Override
public ProjectionPoint latLonToProj(LatLonPoint latlon, ProjectionPointImpl result) {
double fromLat = Math.toRadians(latlon.getLatitude());
double theta = Math.toRadians(latlon.getLongitude());
if (projectionLongitude != 0 && !Double.isNaN(theta)) {
theta = MapMath.normalizeLongitude(theta - projectionLongitude);
}
ProjectionPointImpl out = new ProjectionPointImpl();
project(theta, fromLat, out);
result.setLocation(totalScale * out.getX() + falseEasting, totalScale * out.getY() + falseNorthing);
return result;
}
/**
* Convert projection coordinates to a LatLonPoint
* Note: a new object is not created on each call for the return value.
*
* @param world convert from these projection coordinates
* @param result the object to write to
* @return LatLonPoint convert to these lat/lon coordinates
*/
@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;
ProjectionPointImpl pp = new ProjectionPointImpl();
projectInverse(fromX, fromY, pp);
if (pp.getX() < -Math.PI) {
pp.setX(-Math.PI);
} else if (pp.getX() > Math.PI) {
pp.setX(Math.PI);
}
if (projectionLongitude != 0 && !Double.isNaN(pp.getX())) {
pp.setX(MapMath.normalizeLongitude(pp.getX() + projectionLongitude));
}
result.setLatitude( Math.toDegrees(pp.getY()));
result.setLongitude( Math.toDegrees(pp.getX()));
return result;
}
@Override
public ProjectionImpl constructCopy() {
ProjectionImpl result = new PolyconicProjection(getOriginLatitude(), getOriginLongitude(), getFalseEasting(), getFalseNorthing(), getEarth());
result.setDefaultMapArea(defaultMapArea);
result.setName(name);
return result;
}
}
© 2015 - 2024 Weber Informatics LLC | Privacy Policy