Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/CRSFactory.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/CRSFactory.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/CRSFactory.java (revision 27902) @@ -8,4 +8,6 @@ * from a variety of ways * of specifying them. + * This is the primary way of creating coordinate systems + * for carrying out projections transformations. *

* CoordinateReferenceSystems can be used to @@ -42,20 +44,26 @@ /** - * Creates a {@link CoordinateReferenceSystem} from a well-known name. - * Names are of the form: "authority:code". + * Creates a {@link CoordinateReferenceSystem} (CRS) from a well-known name. + * CRS names are of the form: "authority:code", + * with the components being: *

- * An example of a valid name is EPSG:3005. - * + * An example of a valid CRS name is EPSG:3005. + *

* @param name the name of a coordinate system, with optional authority prefix - * @return a {@link CoordinateReferenceSystem} + * @return the {@link CoordinateReferenceSystem} corresponding to the given name * @throws UnsupportedParameterException if a PROJ.4 parameter is not supported * @throws InvalidValueException if a parameter value is invalid @@ -73,15 +81,15 @@ /** * Creates a {@link CoordinateReferenceSystem} - * defined by a PROJ.4 parameter string. + * from a PROJ.4 projection parameter string. *

- * An example of a valid PROJ.4 parameter string is: + * An example of a valid PROJ.4 projection parameter string is: * 

    * +proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +ellps=GRS80 +units=m
    *
- * @param name a name for this coordinate system (may be null) - * @param paramStr a PROJ.4 parameter string - * @return a {@link CoordinateReferenceSystem} - * @throws UnsupportedParameterException if a PROJ.4 parameter is not supported - * @throws InvalidValueException if a parameter value is invalid + * @param name a name for this coordinate system (may be null for an anonymous coordinate system) + * @param paramStr a PROJ.4 projection parameter string + * @return the specified {@link CoordinateReferenceSystem} + * @throws UnsupportedParameterException if a given PROJ.4 parameter is not supported + * @throws InvalidValueException if a supplied parameter value is invalid */ public CoordinateReferenceSystem createFromParameters(String name, String paramStr) @@ -93,20 +101,22 @@ /** * Creates a {@link CoordinateReferenceSystem} - * defined by an array of PROJ.4 parameters. + * defined by an array of PROJ.4 projection parameters. + * PROJ.4 parameters are generally of the form + * "+name=value". * * @param name a name for this coordinate system (may be null) - * @param args an array of PROJ.4 parameters + * @param params an array of PROJ.4 projection parameters * @return a {@link CoordinateReferenceSystem} * @throws UnsupportedParameterException if a PROJ.4 parameter is not supported * @throws InvalidValueException if a parameter value is invalid */ - public CoordinateReferenceSystem createFromParameters(String name, String[] args) + public CoordinateReferenceSystem createFromParameters(String name, String[] params) throws UnsupportedParameterException, InvalidValueException { - if (args == null) + if (params == null) return null; Proj4Parser parser = new Proj4Parser(registry); - return parser.parse(name, args); + return parser.parse(name, params); } Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/ProjCoordinate.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/ProjCoordinate.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/ProjCoordinate.java (revision 27902) @@ -1,3 +1,5 @@ package org.osgeo.proj4j; + +import java.text.DecimalFormat; /** @@ -14,4 +16,7 @@ public class ProjCoordinate { + public static String DECIMAL_FORMAT_PATTERN = "0.0###############"; + public static DecimalFormat DECIMAL_FORMAT = new DecimalFormat(DECIMAL_FORMAT_PATTERN); + /** * The X ordinate for this point. @@ -317,7 +322,7 @@ StringBuilder builder = new StringBuilder(); builder.append("["); - builder.append(this.x); + builder.append(DECIMAL_FORMAT.format(x)); builder.append(", "); - builder.append(this.y); + builder.append(DECIMAL_FORMAT.format(y)); if (! Double.isNaN(z)) { builder.append(", "); Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/ProjectionException.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/ProjectionException.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/ProjectionException.java (revision 27902) @@ -21,5 +21,5 @@ /** * Signals that an erroneous situation has - * occurred during the computation of + * occured during the computation of * a projected coordinate system value. * Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/Registry.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/Registry.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/Registry.java (revision 27902) @@ -9,5 +9,5 @@ /** - * Records predefined values for various library classes + * Supplies predefined values for various library classes * such as {@link Ellipsoid}, {@link Datum}, and {@link Projection}. * @@ -176,5 +176,5 @@ // register( "kav7", Projection.class, "Kavraisky VII" ); // register( "labrd", Projection.class, "Laborde" ); -// register( "laea", Projection.class, "Lambert Azimuthal Equal Area" ); + register( "laea", LambertAzimuthalEqualAreaProjection.class, "Lambert Azimuthal Equal Area" ); register( "lagrng", LagrangeProjection.class, "Lagrange" ); register( "larr", LarriveeProjection.class, "Larrivee" ); @@ -226,5 +226,5 @@ register( "rpoly", RectangularPolyconicProjection.class, "Rectangular Polyconic" ); register( "sinu", SinusoidalProjection.class, "Sinusoidal (Sanson-Flamsteed)" ); -// register( "somerc", Projection.class, "Swiss. Obl. Mercator" ); + register( "somerc", SwissObliqueMercatorProjection.class, "Swiss Oblique Mercator" ); register( "stere", StereographicAzimuthalProjection.class, "Stereographic" ); register( "sterea", ObliqueStereographicAlternativeProjection.class, "Oblique Stereographic Alternative" ); Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/datum/Datum.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/datum/Datum.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/datum/Datum.java (revision 27902) @@ -36,5 +36,8 @@ *
  • A grid-shift conversion * - * + * In order to be able to transform between any two datums, + * the parameter-based transforms are provided as a transform to + * the common WGS84 datum. The WGS transforms of two arbitrary datum transforms can + * be concatenated to provide a transform between the two datums. *

    * Notable datums in common use include {@link #NAD83} and {@link #WGS84}. @@ -121,15 +124,30 @@ if (transform == null) return TYPE_WGS84; - if (isZero(transform)) return TYPE_WGS84; + if (isIdentity(transform)) return TYPE_WGS84; if (transform.length == 3) return TYPE_3PARAM; if (transform.length == 7) return TYPE_7PARAM; + return TYPE_WGS84; } - private static boolean isZero(double[] transform) + /** + * Tests whether the datum parameter-based transform + * is the identity transform + * (in which case datum transformation can be short-circuited, + * thus avoiding some loss of numerical precision). + * + * @param transform + * @return + */ + private static boolean isIdentity(double[] transform) { for (int i = 0; i < transform.length; i++) { - if (transform[i] != 0.0) return false; + // scale factor will normally be 1 for an identity transform + if (i == 6) { + if (transform[i] != 1.0 && transform[i] != 0.0) + return false; + } + else if (transform[i] != 0.0) return false; } return true; @@ -138,5 +156,6 @@ public boolean hasTransformToWGS84() { - return getTransformType() == TYPE_3PARAM || getTransformType() == TYPE_7PARAM; + int transformType = getTransformType(); + return transformType == TYPE_3PARAM || transformType == TYPE_7PARAM; } @@ -187,5 +206,4 @@ } - public void transformFromGeocentricToWgs84(ProjCoordinate p) { Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/io/Proj4FileReader.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/io/Proj4FileReader.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/io/Proj4FileReader.java (revision 27902) @@ -16,5 +16,5 @@ } - public String[] readFile( String file, String name ) + public String[] readParametersFromFile( String authorityCode, String name ) throws IOException { @@ -23,8 +23,8 @@ // TODO: parse CSes line-at-a-time (this allows preserving CS param string for later access) - String filename = "/nad/" + file.toLowerCase(); + String filename = "/nad/" + authorityCode.toLowerCase(); InputStream inStr = Proj4FileReader.class.getResourceAsStream( filename ); if (inStr == null) { - throw new IllegalStateException("Unable to access resource " + filename); + throw new IllegalStateException("Unable to access CRS file: " + filename); } BufferedReader reader = new BufferedReader( @@ -41,6 +41,5 @@ } - private String[] readFile( BufferedReader reader, String name) - throws IOException + private StreamTokenizer createTokenizer(BufferedReader reader) { StreamTokenizer t = new StreamTokenizer( reader ); @@ -59,4 +58,12 @@ t.wordChars( ',', ',' ); t.wordChars( '@', '@' ); + return t; + } + + private String[] readFile( BufferedReader reader, String name) + throws IOException + { + StreamTokenizer t = createTokenizer(reader); + t.nextToken(); while ( t.ttype == '<' ) { @@ -64,5 +71,5 @@ if ( t.ttype != StreamTokenizer.TT_WORD ) throw new IOException( t.lineno()+": Word expected after '<'" ); - String cname = t.sval; + String crsName = t.sval; t.nextToken(); if ( t.ttype != '>' ) @@ -99,5 +106,7 @@ throw new IOException( t.lineno()+": '<>' expected" ); t.nextToken(); - if ( cname.equals( name ) ) { + + // found requested CRS? + if ( crsName.equals( name ) ) { String[] args = (String[]) v.toArray( new String[0] ); return args; @@ -117,23 +126,23 @@ else v.add(plusKey); - - } - public String[] getParameters(String name) { + + /** + * Gets the list of PROJ.4 parameters which define + * the coordinate system specified by name. + * + * @param crsName the name of the coordinate system + * @return the PROJ.4 projection parameters which define the coordinate system + */ + public String[] getParameters(String crsName) { try { - int p = name.indexOf(':'); - if (p >= 0) - return readFile(name.substring(0, p), name.substring(p + 1)); - - /* - // not sure this is needed or a good idea - String[] files = { "world", "nad83", "nad27", "esri", "epsg", }; - for (int i = 0; i < files.length; i++) { - Projection projection = readProjectionFile(files[i], name); - if (projection != null) - return projection; + int p = crsName.indexOf(':'); + if (p >= 0) { + String auth = crsName.substring(0, p); + String id = crsName.substring(p + 1); + return readParametersFromFile(auth, id); } - */ - } catch (IOException e) { + } + catch (IOException e) { e.printStackTrace(); } Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/parser/DatumParameters.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/parser/DatumParameters.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/parser/DatumParameters.java (revision 27902) @@ -11,5 +11,5 @@ * This class also implements the policies for * which parameters take precedence - * when multiple inconsisent ones are present. + * when multiple inconsistent ones are present. * * @author Martin Davis Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/parser/Proj4Parser.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/parser/Proj4Parser.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/parser/Proj4Parser.java (revision 27902) @@ -16,4 +16,8 @@ public class Proj4Parser { + /* SECONDS_TO_RAD = Pi/180/3600 */ + private static final double SECONDS_TO_RAD = 4.84813681109535993589914102357e-6; + private static final double MILLION = 1000000.0; + private Registry registry; @@ -169,5 +173,5 @@ String towgs84 = (String) params.get(Proj4Keyword.towgs84); if (towgs84 != null) { - double[] datumConvParams = parseDatumTransform(towgs84); + double[] datumConvParams = parseToWGS84(towgs84); datumParam.setDatumTransform(datumConvParams); } @@ -183,5 +187,5 @@ } - private double[] parseDatumTransform(String paramList) + private double[] parseToWGS84(String paramList) { String[] numStr = paramList.split(","); @@ -195,4 +199,27 @@ param[i] = Double.parseDouble(numStr[i]); } + + // optimization to detect 3-parameter transform + if (param[3] == 0.0 + && param[4] == 0.0 + && param[5] == 0.0 + && param[6] == 0.0 + ) { + param = new double[] { param[0], param[1], param[2] }; + } + + /** + * PROJ4 towgs84 7-parameter transform uses + * units of arc-seconds for the rotation factors, + * and parts-per-million for the scale factor. + * These need to be converted to radians and a scale factor. + */ + if (param.length > 3) { + param[3] *= SECONDS_TO_RAD; + param[4] *= SECONDS_TO_RAD; + param[5] *= SECONDS_TO_RAD; + param[6] = (param[6]/MILLION) + 1; + } + return param; } Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/AlbersProjection.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/AlbersProjection.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/AlbersProjection.java (revision 27902) @@ -182,9 +182,4 @@ } -/* - public String toString() { - return "Lambert Equal Area Conic"; - } -*/ } Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/LambertAzimuthalEqualAreaProjection.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/LambertAzimuthalEqualAreaProjection.java (revision 27902) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/LambertAzimuthalEqualAreaProjection.java (revision 27902) @@ -0,0 +1,314 @@ +/* +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 org.osgeo.proj4j.proj; + +import org.osgeo.proj4j.ProjCoordinate; +import org.osgeo.proj4j.ProjectionException; +import org.osgeo.proj4j.util.ProjectionMath; + +public class LambertAzimuthalEqualAreaProjection extends Projection { + + private static final int N_POLE = 0; + private static final int S_POLE = 1; + private static final int EQUIT = 2; + private static final int OBLIQ = 3; + + private int mode = 0; + private double phi0; + private double sinb1; + private double cosb1; + private double xmf; + private double ymf; + private double mmf; + private double qp; + private double dd; + private double rq; + private double[] apa; + private double sinph0; + private double cosph0; + + public LambertAzimuthalEqualAreaProjection() { + this( false ); + } + + public LambertAzimuthalEqualAreaProjection( boolean south ) { + //minLatitude = Math.toRadians(0); + //maxLatitude = Math.toRadians(90); + //projectionLatitude1 = south ? -ProjectionMath.QUARTERPI : ProjectionMath.QUARTERPI; + //projectionLatitude2 = south ? -ProjectionMath.HALFPI : ProjectionMath.HALFPI; + //initialize(); + } + + public void initialize() { + super.initialize(); + double t; + + phi0 = projectionLatitude; + + if (Math.abs((t = Math.abs(phi0)) - ProjectionMath.HALFPI) < ProjectionMath.EPS10) { + mode = phi0 < 0. ? S_POLE : N_POLE; + } + else if (Math.abs(t) < ProjectionMath.EPS10) { + mode = EQUIT; + } + else { + mode = OBLIQ; + } + if (!spherical) { + double sinphi; + + e = Math.sqrt(es); + qp = ProjectionMath.qsfn(1., e, one_es); + mmf = .5 / (1. - es); + apa = ProjectionMath.authset(es); + switch (mode) { + case N_POLE: + case S_POLE: + dd = 1.; + break; + case EQUIT: + dd = 1. / (rq = Math.sqrt(.5 * qp)); + xmf = 1.; + ymf = .5 * qp; + break; + case OBLIQ: + rq = Math.sqrt(.5 * qp); + sinphi = Math.sin(phi0); + sinb1 = ProjectionMath.qsfn(sinphi, e, one_es) / qp; + cosb1 = Math.sqrt(1. - sinb1 * sinb1); + dd = Math.cos(phi0) / (Math.sqrt(1. - es * sinphi * sinphi) * + rq * cosb1); + ymf = (xmf = rq) / dd; + xmf *= dd; + break; + } + } else { + if (mode == OBLIQ) { + sinph0 = Math.sin(phi0); + cosph0 = Math.cos(phi0); + } + } + } + + public ProjCoordinate project(double lplam, double lpphi, ProjCoordinate out) { + if (spherical) { + project_s(lplam, lpphi, out); + } + else { + project_e(lplam, lpphi, out); + } + return out; + } + + public void project_s(double lplam, double lpphi, ProjCoordinate out) { + double coslam, cosphi, sinphi; + + sinphi = Math.sin(lpphi); + cosphi = Math.cos(lpphi); + coslam = Math.cos(lplam); + switch (mode) { + case EQUIT: + case OBLIQ: + if (mode == EQUIT) + out.y = 1. + cosphi * coslam; + else + out.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam; + + if (out.y <= EPS10) throw new ProjectionException("F"); + out.x = (out.y = Math.sqrt(2. / out.y)) * cosphi * Math.sin(lplam); + out.y *= mode == EQUIT ? sinphi : + cosph0 * sinphi - sinph0 * cosphi * coslam; + break; + case N_POLE: + coslam = -coslam; + case S_POLE: + if (Math.abs(lpphi + phi0) < EPS10) throw new ProjectionException("F");; + out.y = ProjectionMath.QUARTERPI - lpphi * .5; + out.y = 2. * (mode == S_POLE ? Math.cos(out.y) : Math.sin(out.y)); + out.x = out.y * Math.sin(lplam); + out.y *= coslam; + break; + } + } + + public void project_e(double lplam, double lpphi, ProjCoordinate out) { + double coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0; + + coslam = Math.cos(lplam); + sinlam = Math.sin(lplam); + sinphi = Math.sin(lpphi); + q = ProjectionMath.qsfn(sinphi, e, one_es); + if (mode == OBLIQ || mode == EQUIT) { + sinb = q / qp; + cosb = Math.sqrt(1. - sinb * sinb); + } + switch (mode) { + case OBLIQ: + b = 1. + sinb1 * sinb + cosb1 * cosb * coslam; + break; + case EQUIT: + b = 1. + cosb * coslam; + break; + case N_POLE: + b = ProjectionMath.HALFPI + lpphi; + q = qp - q; + break; + case S_POLE: + b = lpphi - ProjectionMath.HALFPI; + q = qp + q; + break; + } + if (Math.abs(b) < EPS10) throw new ProjectionException("F"); + + switch (mode) { + case OBLIQ: + case EQUIT: + if (mode == OBLIQ) { + out.y = ymf * ( b = Math.sqrt(2. / b) ) + * (cosb1 * sinb - sinb1 * cosb * coslam); + } + else { + out.y = (b = Math.sqrt(2. / (1. + cosb * coslam))) * sinb * ymf; + } + out.x = xmf * b * cosb * sinlam; + break; + case N_POLE: + case S_POLE: + if (q >= 0.) { + out.x = (b = Math.sqrt(q)) * sinlam; + out.y = coslam * (mode == S_POLE ? b : -b); + } else + out.x = out.y = 0.; + break; + } + } + + public ProjCoordinate projectInverse(double xyx, double xyy, ProjCoordinate out) { + if (spherical) { + projectInverse_s(xyx, xyy, out); + } + else { + projectInverse_e(xyx, xyy, out); + } + return out; + } + + public void projectInverse_s(double xyx, double xyy, ProjCoordinate out) { + double cosz=0.0, rh, sinz=0.0; + double lpphi, lplam; + + rh = Math.hypot(xyx, xyy); + if ((lpphi = rh * .5 ) > 1.) throw new ProjectionException("I_ERROR"); + lpphi = 2. * Math.asin(lpphi); + if (mode == OBLIQ || mode == EQUIT) { + sinz = Math.sin(lpphi); + cosz = Math.cos(lpphi); + } + switch (mode) { + case EQUIT: + lpphi = Math.abs(rh) <= EPS10 ? 0. : Math.asin(xyy * sinz / rh); + xyx *= sinz; + xyy = cosz * rh; + break; + case OBLIQ: + lpphi = Math.abs(rh) <= EPS10 ? phi0 : + Math.asin(cosz * sinph0 + xyy * sinz * cosph0 / rh); + xyx *= sinz * cosph0; + xyy = (cosz - Math.sin(lpphi) * sinph0) * rh; + break; + case N_POLE: + xyy = -xyy; + lpphi = ProjectionMath.HALFPI - lpphi; + break; + case S_POLE: + lpphi -= ProjectionMath.HALFPI; + break; + } + lplam = (xyy == 0. && (mode == EQUIT || mode == OBLIQ)) ? + 0. : Math.atan2(xyx, xyy); + out.x = lplam; + out.y = lpphi; + } + + public void projectInverse_e(double xyx, double xyy, ProjCoordinate out) { + double lpphi, lplam; + double cCe, sCe, q, rho, ab=0.0; + + switch (mode) { + case EQUIT: + case OBLIQ: + if ((rho = Math.hypot(xyx /= dd, xyy *= dd)) < EPS10) { + lplam = 0.; + lpphi = phi0; + out.x = lplam; + out.y = lpphi; + return; + } + cCe = Math.cos(sCe = 2. * Math.asin(.5 * rho / rq)); + xyx *= (sCe = Math.sin(sCe)); + if (mode == OBLIQ) { + q = qp * (ab = cCe * sinb1 + xyy * sCe * cosb1 / rho); + xyy = rho * cosb1 * cCe - xyy * sinb1 * sCe; + } else { + q = qp * (ab = xyy * sCe / rho); + xyy = rho * cCe; + } + break; + case N_POLE: + xyy = -xyy; + case S_POLE: + if (0 == (q = (xyx * xyx + xyy * xyy)) ) { + lplam = 0.; + lpphi = phi0; + out.x = lplam; + out.y = lpphi; + return; + } + /* + q = P->qp - q; + */ + ab = 1. - q / qp; + if (mode == S_POLE) + ab = - ab; + break; + } + lplam = Math.atan2(xyx, xyy); + lpphi = ProjectionMath.authlat(Math.asin(ab), apa); + out.x = lplam; + out.y = lpphi; + } + + /** + * Returns true if this projection is equal area + */ + public boolean isEqualArea() { + return true; + } + + public boolean hasInverse() { + return true; + } + + public String toString() { + return "Lambert Azimuthal Equal Area"; + } + + +} Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/SwissObliqueMercatorProjection.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/SwissObliqueMercatorProjection.java (revision 27902) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/SwissObliqueMercatorProjection.java (revision 27902) @@ -0,0 +1,116 @@ +/* +Copyright 2006 Martin Davis + +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 converted from the PROJ.4 source. + */ +package org.osgeo.proj4j.proj; + +import org.osgeo.proj4j.ProjCoordinate; +import org.osgeo.proj4j.ProjectionException; +import org.osgeo.proj4j.datum.Ellipsoid; +import org.osgeo.proj4j.util.ProjectionMath; + +/** +* Swiss Oblique Mercator Projection algorithm is taken from the USGS PROJ.4 package. +*/ +public class SwissObliqueMercatorProjection extends Projection { + + private static final int NITER = 6; + + private double K, c, hlf_e, kR, cosp0, sinp0; + private double phi0; + + public SwissObliqueMercatorProjection() { + //initialize(); + } + + public void initialize() { + super.initialize(); + double cp, phip0, sp; + + phi0 = projectionLatitude; + + hlf_e = 0.5 * e; + cp = Math.cos(phi0); + cp *= cp; + c = Math.sqrt(1 + es * cp * cp * rone_es); + sp = Math.sin(phi0); + cosp0 = Math.cos( phip0 = Math.asin(sinp0 = sp / c) ); + sp *= e; + K = Math.log(Math.tan(ProjectionMath.FORTPI + 0.5 * phip0)) - c * ( + Math.log(Math.tan(ProjectionMath.FORTPI + 0.5 * phi0)) - hlf_e * + Math.log((1. + sp) / (1. - sp))); + kR = scaleFactor * Math.sqrt(one_es) / (1. - sp * sp); + } + + public ProjCoordinate project(double lplam, double lpphi, ProjCoordinate xy) { + double phip, lamp, phipp, lampp, sp, cp; + + sp = e * Math.sin(lpphi); + phip = 2.* Math.atan( Math.exp( c * ( + Math.log(Math.tan(ProjectionMath.FORTPI + 0.5 * lpphi)) - hlf_e * Math.log((1. + sp)/(1. - sp))) + + K)) - ProjectionMath.HALFPI; + lamp = c * lplam; + cp = Math.cos(phip); + phipp = Math.asin(cosp0 * Math.sin(phip) - sinp0 * cp * Math.cos(lamp)); + lampp = Math.asin(cp * Math.sin(lamp) / Math.cos(phipp)); + xy.x = kR * lampp; + xy.y = kR * Math.log(Math.tan(ProjectionMath.FORTPI + 0.5 * phipp)); + return xy; + } + + public ProjCoordinate projectInverse(double xyx, double xyy, ProjCoordinate lp) { + double phip, lamp, phipp, lampp, cp, esp, con, delp; + int i; + double lplam, lpphi; + + phipp = 2. * (Math.atan(Math.exp(xyy / kR)) - ProjectionMath.FORTPI); + lampp = xyx / kR; + cp = Math.cos(phipp); + phip = Math.asin(cosp0 * Math.sin(phipp) + sinp0 * cp * Math.cos(lampp)); + lamp = Math.asin(cp * Math.sin(lampp) / Math.cos(phip)); + con = (K - Math.log(Math.tan(ProjectionMath.FORTPI + 0.5 * phip)))/c; + for (i = NITER; i != 0 ; --i) { + esp = e * Math.sin(phip); + delp = (con + Math.log(Math.tan(ProjectionMath.FORTPI + 0.5 * phip)) - hlf_e * + Math.log((1. + esp)/(1. - esp)) ) * + (1. - esp * esp) * Math.cos(phip) * rone_es; + phip -= delp; + if (Math.abs(delp) < ProjectionMath.EPS10) + break; + } + if (i != 0) { + lpphi = phip; + lplam = lamp / c; + } else { + throw new ProjectionException("I_ERROR"); + } + lp.x = lplam; + lp.y = lpphi; + + return lp; + } + + public boolean hasInverse() { + return true; + } + + public String toString() { + return "Swiss Oblique Mercator"; + } + +} Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/TransverseMercatorProjection.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/TransverseMercatorProjection.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/proj/TransverseMercatorProjection.java (revision 27902) @@ -17,9 +17,4 @@ /* * This file was semi-automatically converted from the public-domain USGS PROJ source. - * - * Bernhard Jenny, February 2 2010: Corrected code for spherical case in - * projectInverse, added isConformal. - * 27 September 2010: added missing tests to forward spherical, removed - * initialization code in constructor. */ package org.osgeo.proj4j.proj; @@ -160,22 +155,11 @@ public ProjCoordinate projectInverse(double x, double y, ProjCoordinate out) { if (spherical) { - /* - Original code - x = Math.exp(x / scaleFactor); - y = .5 * (x - 1. / x); - x = Math.cos(projectionLatitude + y / scaleFactor); - out.y = MapMath.asin(Math.sqrt((1. - x * x) / (1. + y * y))); - if (y < 0) { - out.y = -out.y; - } - out.x = Math.atan2(y, x); - */ - - // new code by Bernhard Jenny, February 2 2010 - double D = y / scaleFactor + projectionLatitude; - double xp = x / scaleFactor; - - out.y = Math.asin(Math.sin(D) / Math.cosh(xp)); - out.x = Math.atan2(Math.sinh(xp), Math.cos(D)); + double h = Math.exp(x / scaleFactor); + double g = .5 * (h - 1. / h); + h = Math.cos(projectionLatitude + y / scaleFactor); + out.y = ProjectionMath.asin(Math.sqrt((1. - h*h) / (1. + g*g))); + if (y < 0) + out.y = -out.y; + out.x = Math.atan2(g, h); } else { double n, con, cosphi, d, ds, sinphi, t; @@ -213,8 +197,4 @@ return true; } - - public boolean isConformal() { - return true; - } public String toString() { Index: /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/util/ProjectionMath.java =================================================================== --- /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/util/ProjectionMath.java (revision 27901) +++ /applications/editors/josm/plugins/proj4j/src/org/osgeo/proj4j/util/ProjectionMath.java (revision 27902) @@ -27,5 +27,6 @@ public final static double PI = Math.PI; public final static double HALFPI = Math.PI/2.0; - public final static double QUARTERPI = Math.PI/4.0; + public final static double QUARTERPI = Math.PI/4.0; + public final static double FORTPI = QUARTERPI; public final static double TWOPI = Math.PI*2.0; public final static double RTD = 180.0/Math.PI; @@ -34,4 +35,5 @@ public final static Rectangle2D WORLD_BOUNDS = new Rectangle2D.Double(-180, -90, 360, 180); + public final static double EPS10 = 1.0e-10; /** * Degree versions of trigonometric functions