I added the Bonne projection, which is a generalized form of Sinusoidal and Werner (I don't have Werner). As it turns out, this was one of the go-to map projections back in the days before Albers and Van der Grinten had been invented.
This commit is contained in:
Justin Kunimune 2019-09-24 19:59:29 -04:00
parent 3e819d9fdb
commit c4df39aa6c
4 changed files with 128 additions and 6 deletions

View File

@ -197,9 +197,11 @@ public class Conic {
lat = -lat;
lon = -lon;
}
final double s = reversed ? -1 : 1;
final double r = Math.sqrt(C - 2*n*Math.sin(lat));
return new double[] { s*r*Math.sin(n*lon), s*(y0 - r*Math.cos(n*lon)) };
final double x = r*Math.sin(n*lon);
final double y = -r*Math.cos(n*lon);
if (reversed) return new double[] {-x,-y};
else return new double[] { x, y};
}
public double[] inverse(double x, double y) {
@ -219,7 +221,6 @@ public class Conic {
};
/**
* A base for all conic projections
* @author jkunimune

View File

@ -233,6 +233,89 @@ public class Misc {
};
public static final Projection BONNE =
new Projection(
"Bonne", "A traditional pseudoconic projection, also known as the Sylvanus projection.",
10, 10, 0b1111, Type.PSEUDOCONIC, Property.EQUAL_AREA, 1,
new String[] {"Std. Parallel"}, new double[][] {{-90, 90, 45}}) {
private double r0;
private double yC; // the y coordinate of the centers of the parallels
private boolean reversed; // if the standard parallel is southern
public void setParameters(double... params) {
double lat0 = Math.toRadians(params[0]);
this.reversed = (lat0 < 0);
if (reversed)
lat0 = -lat0;
this.r0 = 1/Math.tan(lat0) + lat0;
if (Double.isInfinite(r0)) {
this.width = Pseudocylindrical.SINUSOIDAL.getWidth();
this.height = Pseudocylindrical.SINUSOIDAL.getHeight();
}
else { // for such a simple map projection...
double argmaxX = NumericalAnalysis.bisectionFind(
(p) -> (Math.PI*(1/(r0-p)*Math.cos(p) - Math.sin(p))*Math.cos(Math.PI/(r0-p)*Math.cos(p)) - Math.sin(Math.PI/(r0-p)*Math.cos(p))),
-Math.PI/2, 0, 1e-3); // it sure is complicated to find its dimensions!
double maxX = (r0 - argmaxX)*Math.sin(Math.PI/(r0 - argmaxX)*Math.cos(argmaxX));
this.width = 2*maxX;
double argmaxY;
try {
argmaxY = NumericalAnalysis.bisectionFind(
(p) -> (Math.PI*(1/(r0-p)*Math.cos(p) - Math.sin(p))*Math.sin(Math.PI/(r0-p)*Math.cos(p)) + Math.cos(Math.PI/(r0-p)*Math.cos(p))),
0, Math.PI/4, 1e-3);
} catch (IllegalArgumentException e) {
argmaxY = Math.PI/2;
}
double maxY = Math.max(-r0 + Math.PI/2,
-(r0 - argmaxY)*Math.cos(Math.PI/(r0 - argmaxY)*Math.cos(argmaxY)));
this.height = r0 + Math.PI/2 + maxY;
this.yC = (r0 + Math.PI/2 - maxY)/2;
}
}
public double[] project(double lat, double lon) {
if (Double.isInfinite(r0))
return Pseudocylindrical.SINUSOIDAL.project(lat, lon);
if (reversed) {
lat = -lat;
lon = -lon;
}
double r = r0 - lat;
double th = lon*Math.cos(lat)/r;
double x = r*Math.sin(th);
double y = yC - r*Math.cos(th);
if (reversed)
return new double[] {-x,-y};
else
return new double[] { x, y};
}
public double[] inverse(double x, double y) {
if (Double.isInfinite(r0))
return Pseudocylindrical.SINUSOIDAL.inverse(x, y);
if (reversed) {
x = -x;
y = -y;
}
double r = Math.hypot(x, y-yC);
if (r < r0 - Math.PI/2 || r > r0 + Math.PI/2)
return null;
double th = Math.atan2(x, -(y-yC));
double lat = r0 - r;
double lon = th*r/Math.cos(lat);
if (reversed)
return new double[] {-lat,-lon};
else
return new double[] { lat, lon};
}
};
public static final Projection T_SHIRT =
new Projection(
"T-Shirt", "A conformal projection onto a torso.",

View File

@ -560,6 +560,14 @@ public abstract class Projection {
return this.property;
}
/**
* @return the rating:
* 0 if I hate it with a burning passion;
* 1 if it is bad and you shouldn't use it;;
* 2 if it has its use cases but isn't very good outside of them;
* 3 if it is a solid, defensible choice; and
* 4 if I love it with a burning passion.
*/
public final int getRating() {
return this.rating;
}
@ -612,10 +620,10 @@ public abstract class Projection {
public static enum Type {
CYLINDRICAL("Cylindrical"), CONIC("Conic"), AZIMUTHAL("Azimuthal"),
PSEUDOCYLINDRICAL("Pseudocylindrical"), PSEUDOCONIC("Pseudoconic"),
PSEUDOAZIMUTHAL("Pseudoazimuthal"), QUASIAZIMUTHAL("Quasiazimuthal"),
PSEUDOAZIMUTHAL("Pseudoazimuthal"),
TETRAHEDRAL("Tetrahedral"), OCTOHEDRAL("Octohedral"),
TETRADECAHEDRAL("Truncated Octohedral"), ICOSOHEDRAL("Icosohedral"),
POLYNOMIAL("Polynomial"), STREBE("Strebe Blend"), PLANAR("Planar"), OTHER("Other");
POLYNOMIAL("Polynomial"), STREBE("Strebe blend"), PLANAR("Planar"), OTHER("Other");
private String name;

View File

@ -138,6 +138,36 @@ public class NumericalAnalysis {
}
/**
* Applies a bisection zero-finding scheme to find x such that f(x)=y
* @param f The function whose zero must be found
* @param xMin The lower bound for the zero
* @param xMax The upper bound for the zero
* @param tolerence The maximum error in x that this can return
* @return The value of x that sets f to zero
*/
public static final double bisectionFind(DoubleUnaryOperator f,
double xMin, double xMax, double tolerance) {
double yMin = f.applyAsDouble(xMin);
double yMax = f.applyAsDouble(xMax);
if ((yMin < 0) == (yMax < 0))
throw new IllegalArgumentException("Bisection failed; bounds "+xMin+" and "+xMax+" do not necessarily straddle a zero.");
while (Math.abs(xMax - xMin) > tolerance) {
double x = (xMax + xMin)/2;
double y = f.applyAsDouble(x);
if ((y < 0) == (yMin < 0)) {
xMin = x;
yMin = y;
}
else {
xMax = x;
yMax = y;
}
}
return (xMax + xMin)/2;
}
/**
* Applies Newton's method in one dimension to solve for x such that f(x)=y
* @param y Desired value for f
@ -145,7 +175,7 @@ public class NumericalAnalysis {
* @param f The error in terms of x
* @param dfdx The derivative of f with respect to x
* @param tolerance The maximum error that this can return
* @return The value of x that puts f near 0, or NaN if it does not converge in 5 iterations
* @return The value of x that puts f near y, or NaN if it does not converge in 5 iterations
*/
public static final double newtonRaphsonApproximation(
double y, double x0, DoubleUnaryOperator f, DoubleUnaryOperator dfdx, double tolerance) {