root / org.gvsig.jcrs / libJCRS / src / org / geotools / referencing / operation / projection / IdrNewZealandMapGrid.java @ 38
History | View | Annotate | Download (17.6 KB)
| 1 |
/*
|
|---|---|
| 2 |
* GeoTools - OpenSource mapping toolkit
|
| 3 |
* http://geotools.org
|
| 4 |
*
|
| 5 |
* (C) 2005-2006, Geotools Project Managment Committee (PMC)
|
| 6 |
*
|
| 7 |
* This library is free software; you can redistribute it and/or
|
| 8 |
* modify it under the terms of the GNU Lesser General Public
|
| 9 |
* License as published by the Free Software Foundation;
|
| 10 |
* version 2.1 of the License.
|
| 11 |
*
|
| 12 |
* This library is distributed in the hope that it will be useful,
|
| 13 |
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
| 14 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
| 15 |
* Lesser General Public License for more details.
|
| 16 |
*/
|
| 17 |
package org.geotools.referencing.operation.projection; |
| 18 |
|
| 19 |
// J2SE dependencies
|
| 20 |
import java.awt.geom.Point2D; |
| 21 |
import java.util.Collection; |
| 22 |
import javax.units.NonSI; |
| 23 |
import javax.units.SI; |
| 24 |
import javax.units.Unit; |
| 25 |
|
| 26 |
import org.geotools.metadata.iso.citation.CitationImpl; |
| 27 |
import org.geotools.referencing.NamedIdentifier; |
| 28 |
import org.opengis.metadata.Identifier; |
| 29 |
import org.opengis.parameter.ParameterDescriptor; |
| 30 |
import org.opengis.parameter.ParameterDescriptorGroup; |
| 31 |
import org.opengis.parameter.ParameterNotFoundException; |
| 32 |
import org.opengis.parameter.ParameterValueGroup; |
| 33 |
import org.opengis.referencing.operation.MathTransform; |
| 34 |
import org.opengis.spatialschema.geometry.complex.Complex; |
| 35 |
|
| 36 |
|
| 37 |
/**
|
| 38 |
* Implementation of the NZMG (New Zealand Map Grid) projection.
|
| 39 |
* <p>
|
| 40 |
* This is an implementation of algorithm published by
|
| 41 |
* <a href="http://www.govt.nz/record?recordid=28">Land Information New Zealand</a>.
|
| 42 |
* The algorithm is documented <a href="http://www.linz.govt.nz/rcs/linz/6137/">here</a>.
|
| 43 |
*
|
| 44 |
* @since 2.2
|
| 45 |
* @source $URL: http://svn.geotools.org/geotools/tags/2.3.0/module/referencing/src/org/geotools/referencing/operation/projection/NewZealandMapGrid.java $
|
| 46 |
* @version $Id: IdrNewZealandMapGrid.java 12357 2007-06-27 09:05:41Z dguerrero $
|
| 47 |
* @author Justin Deoliveira
|
| 48 |
* @author Martin Desruisseaux
|
| 49 |
*
|
| 50 |
* @todo The algorithm uses complex numbers, which is not very well supported in Java. This
|
| 51 |
* implementation uses {@linkplain Complex} as a support class. Various instances of
|
| 52 |
* {@link Complex} are created once for ever at {@code NewZealandMapGrid} construction
|
| 53 |
* time, in order to avoid creating up to 6 objects for every point to be projected.
|
| 54 |
* The downside is that transformation methods must be synchronized. The cost should
|
| 55 |
* be small for simple applications, but may become important for multi-thread applications.
|
| 56 |
* Furthermore, those fields raise a slight serialization issue.
|
| 57 |
* <p>
|
| 58 |
* The most efficient fix in Java would be to expand inline all {@link Complex} operations
|
| 59 |
* like {@link Complex#add add} (easy), {@link Complex#multiply multiply} (more tedious),
|
| 60 |
* <cite>etc.</cite>, until we get a code using only {@code double} primitives on the stack
|
| 61 |
* and no {@link Complex} objects on the heap (except the {@code A} and {@code B} constants).
|
| 62 |
* But it would make the code significantly more complex and difficult to read.
|
| 63 |
* <p>
|
| 64 |
* An elegant fix would have been "lightweight objects" allocated on the stack (something
|
| 65 |
* similar to {@code struct} in C#), if such thing existed in the Java language.
|
| 66 |
*/
|
| 67 |
public class IdrNewZealandMapGrid extends MapProjection { |
| 68 |
/**
|
| 69 |
* For compatibility with different versions during deserialization.
|
| 70 |
*/
|
| 71 |
private static final long serialVersionUID = 8394817836243729133L; |
| 72 |
|
| 73 |
/**
|
| 74 |
* Coefficients for forward and inverse projection.
|
| 75 |
*/
|
| 76 |
/*private static final Complex[] A = {
|
| 77 |
new Complex( 0.7557853228, 0.0 ),
|
| 78 |
new Complex( 0.249204646, 0.003371507 ),
|
| 79 |
new Complex( -0.001541739, 0.041058560 ),
|
| 80 |
new Complex( -0.10162907, 0.01727609 ),
|
| 81 |
new Complex( -0.26623489, -0.36249218 ),
|
| 82 |
new Complex( -0.6870983, -1.1651967 )
|
| 83 |
};*/
|
| 84 |
|
| 85 |
/**
|
| 86 |
* Coefficients for inverse projection.
|
| 87 |
*/
|
| 88 |
/*private static final Complex[] B = {
|
| 89 |
new Complex( 1.3231270439, 0.0 ),
|
| 90 |
new Complex( -0.577245789, -0.007809598 ),
|
| 91 |
new Complex( 0.508307513, -0.112208952 ),
|
| 92 |
new Complex( -0.15094762, 0.18200602 ),
|
| 93 |
new Complex( 1.01418179, 1.64497696 ),
|
| 94 |
new Complex( 1.9660549, 2.5127645 )
|
| 95 |
};*/
|
| 96 |
|
| 97 |
/**
|
| 98 |
* Coefficients for inverse projection.
|
| 99 |
*/
|
| 100 |
private static final double[] TPHI = new double[] { |
| 101 |
1.5627014243, 0.5185406398, -0.03333098, -0.1052906, -0.0368594, 0.007317, |
| 102 |
0.01220, 0.00394, -0.0013 |
| 103 |
}; |
| 104 |
|
| 105 |
/**
|
| 106 |
* Coefficients for forward projection.
|
| 107 |
*/
|
| 108 |
private static final double[] TPSI = new double[] { |
| 109 |
0.6399175073, -0.1358797613, 0.063294409, -0.02526853, 0.0117879, |
| 110 |
-0.0055161, 0.0026906, -0.001333, 0.00067, -0.00034 |
| 111 |
}; |
| 112 |
|
| 113 |
/**
|
| 114 |
* A temporary complex number used during transform calculation. Created once for
|
| 115 |
* ever in order to avoid new object creation for every point to be transformed.
|
| 116 |
*/
|
| 117 |
/*private transient final Complex theta = new Complex();*/
|
| 118 |
|
| 119 |
/**
|
| 120 |
* An other temporary complex number created once for ever for the same reason than
|
| 121 |
* {@link #theta}. This number is usually equals to some other complex number raised
|
| 122 |
* to some power.
|
| 123 |
*/
|
| 124 |
/*private transient final Complex power = new Complex();*/
|
| 125 |
|
| 126 |
/**
|
| 127 |
* An other temporary complex number created once for ever for the same reason than
|
| 128 |
* {@link #theta}.
|
| 129 |
*
|
| 130 |
* @todo Need to reassign those fields on deserialization.
|
| 131 |
*/
|
| 132 |
|
| 133 |
/*private transient final Complex z=new Complex(), t=new Complex(),
|
| 134 |
num=new Complex(), denom=new Complex();*/
|
| 135 |
|
| 136 |
/**
|
| 137 |
* Constructs a new map projection with default parameter values.
|
| 138 |
*/
|
| 139 |
protected IdrNewZealandMapGrid() {
|
| 140 |
this((ParameterValueGroup) Provider.PARAMETERS.createValue()); |
| 141 |
} |
| 142 |
|
| 143 |
/**
|
| 144 |
* Constructs a new map projection from the supplied parameters.
|
| 145 |
*
|
| 146 |
* @param parameters The parameter values in standard units.
|
| 147 |
* @throws ParameterNotFoundException if a mandatory parameter is missing.
|
| 148 |
*/
|
| 149 |
protected IdrNewZealandMapGrid(final ParameterValueGroup parameters) |
| 150 |
throws ParameterNotFoundException
|
| 151 |
{
|
| 152 |
super(parameters);
|
| 153 |
} |
| 154 |
|
| 155 |
/**
|
| 156 |
* {@inheritDoc}
|
| 157 |
*/
|
| 158 |
public ParameterDescriptorGroup getParameterDescriptors() {
|
| 159 |
return Provider.PARAMETERS; |
| 160 |
} |
| 161 |
|
| 162 |
/**
|
| 163 |
* Must be overrided because {@link Provider} uses instances of
|
| 164 |
* {@link ModifiedParameterDescriptor}. This hack was needed because the New Zeland map
|
| 165 |
* projection uses particular default values for parameters like "False Easting", etc.
|
| 166 |
*/
|
| 167 |
/*final boolean isExpectedParameter(final Collection expected, final ParameterDescriptor param) {
|
| 168 |
return ModifiedParameterDescriptor.contains(expected, param);
|
| 169 |
}*/
|
| 170 |
|
| 171 |
/**
|
| 172 |
* Transforms the specified (<var>x</var>,<var>y</var>) coordinate (units in radians)
|
| 173 |
* and stores the result in {@code ptDst} (linear distance on a unit sphere).
|
| 174 |
*/
|
| 175 |
protected synchronized Point2D transformNormalized(final double x, final double y, |
| 176 |
final Point2D ptDst) |
| 177 |
throws ProjectionException
|
| 178 |
{
|
| 179 |
/*
|
| 180 |
final double dphi = (y - latitudeOfOrigin) * (180/Math.PI * 3600E-5);
|
| 181 |
double dphi_pow_i = dphi;
|
| 182 |
double dpsi = 0;
|
| 183 |
for (int i=0; i<TPSI.length; i++) {
|
| 184 |
dpsi += (TPSI[i] * dphi_pow_i);
|
| 185 |
dphi_pow_i *= dphi;
|
| 186 |
}
|
| 187 |
power.real = theta.real = dpsi;
|
| 188 |
power.imag = theta.imag = x;
|
| 189 |
z.multiply(A[0], power);
|
| 190 |
for (int i=1; i<A.length; i++) {
|
| 191 |
power.multiply(power, theta);
|
| 192 |
z.addMultiply(z, A[i], power);
|
| 193 |
}
|
| 194 |
if (ptDst != null) {
|
| 195 |
ptDst.setLocation(z.imag, z.real);
|
| 196 |
return ptDst;
|
| 197 |
}
|
| 198 |
return new Point2D.Double(z.imag, z.real);
|
| 199 |
*/
|
| 200 |
return new Point2D.Double(0.0, 0.0); |
| 201 |
} |
| 202 |
|
| 203 |
/**
|
| 204 |
* Transforms the specified (<var>x</var>,<var>y</var>) coordinate
|
| 205 |
* and stores the result in {@code ptDst}.
|
| 206 |
*/
|
| 207 |
protected synchronized Point2D inverseTransformNormalized(final double x, final double y, |
| 208 |
final Point2D ptDst) |
| 209 |
throws ProjectionException
|
| 210 |
{
|
| 211 |
/*
|
| 212 |
power.real = z.real = y;
|
| 213 |
power.imag = z.imag = x;
|
| 214 |
theta.multiply(B[0], z);
|
| 215 |
for (int j=1; j<B.length; j++) {
|
| 216 |
power.multiply(power, z);
|
| 217 |
theta.addMultiply(theta, B[j], power);
|
| 218 |
}
|
| 219 |
// increasing the number of iterations through this loop decreases
|
| 220 |
// the error in the calculation, but 3 iterations gives 10-3 accuracy
|
| 221 |
for (int j=0; j<3; j++) {
|
| 222 |
power.power(theta, 2);
|
| 223 |
num.addMultiply(z, A[1], power);
|
| 224 |
for (int k=2; k<A.length; k++) {
|
| 225 |
power.multiply(power, theta);
|
| 226 |
t.multiply(A[k], power);
|
| 227 |
t.multiply(t, k);
|
| 228 |
num.add(num, t);
|
| 229 |
}
|
| 230 |
|
| 231 |
power.real = 1;
|
| 232 |
power.imag = 0;
|
| 233 |
denom.copy(A[0]);
|
| 234 |
for (int k=1; k<A.length; k++) {
|
| 235 |
power.multiply(power, theta);
|
| 236 |
t.multiply(A[k], power);
|
| 237 |
t.multiply(t, k+1);
|
| 238 |
denom.add(denom, t);
|
| 239 |
}
|
| 240 |
theta.divide(num, denom);
|
| 241 |
}
|
| 242 |
|
| 243 |
final double dpsi = theta.real;
|
| 244 |
double dpsi_pow_i = dpsi;
|
| 245 |
double dphi = TPHI[0] * dpsi;
|
| 246 |
for (int i=1; i<TPHI.length; i++) {
|
| 247 |
dpsi_pow_i *= dpsi;
|
| 248 |
dphi += (TPHI[i] * dpsi_pow_i);
|
| 249 |
}
|
| 250 |
|
| 251 |
dphi = dphi / (180/Math.PI * 3600E-5) + latitudeOfOrigin;
|
| 252 |
if (ptDst != null) {
|
| 253 |
ptDst.setLocation(theta.imag, dphi);
|
| 254 |
return ptDst;
|
| 255 |
}
|
| 256 |
return new Point2D.Double(theta.imag, dphi);
|
| 257 |
*/
|
| 258 |
return new Point2D.Double(0.0,0.0); |
| 259 |
} |
| 260 |
|
| 261 |
/**
|
| 262 |
* The {@link org.geotools.referencing.operation.MathTransformProvider} for
|
| 263 |
* {@linkplain IdrNewZealandMapGrid New Zealand Map Grid}.
|
| 264 |
*
|
| 265 |
* @since 2.2
|
| 266 |
* @version $Id: IdrNewZealandMapGrid.java 12357 2007-06-27 09:05:41Z dguerrero $
|
| 267 |
* @author Justin Deoliveira
|
| 268 |
*/
|
| 269 |
public static class Provider extends AbstractProvider { |
| 270 |
/**
|
| 271 |
* For compatibility with different versions during deserialization.
|
| 272 |
*/
|
| 273 |
private static final long serialVersionUID = -7716733400419275656L; |
| 274 |
|
| 275 |
|
| 276 |
/**
|
| 277 |
* The operation parameter descriptor for the {@link #semiMajor semiMajor} parameter value.
|
| 278 |
* Valid values range is from 0 to infinity. This parameter is mandatory.
|
| 279 |
*
|
| 280 |
* @todo Would like to start range from 0 <u>exclusive</u>.
|
| 281 |
*/
|
| 282 |
public static final ParameterDescriptor SEMI_MAJOR_LOCAL = createDescriptor( |
| 283 |
new NamedIdentifier[] { |
| 284 |
new NamedIdentifier(CitationImpl.OGC, "semi_major"), |
| 285 |
new NamedIdentifier(CitationImpl.EPSG, "semi-major axis") //epsg does not specifically define this parameter |
| 286 |
}, |
| 287 |
6378388.0, 6378388.0, 6378388.0, SI.METER); |
| 288 |
|
| 289 |
/**
|
| 290 |
* The operation parameter descriptor for the {@link #semiMinor semiMinor} parameter value.
|
| 291 |
* Valid values range is from 0 to infinity. This parameter is mandatory.
|
| 292 |
*
|
| 293 |
* @todo Would like to start range from 0 <u>exclusive</u>.
|
| 294 |
*/
|
| 295 |
public static final ParameterDescriptor SEMI_MINOR_LOCAL = createDescriptor( |
| 296 |
new NamedIdentifier[] { |
| 297 |
new NamedIdentifier(CitationImpl.OGC, "semi_minor"), |
| 298 |
new NamedIdentifier(CitationImpl.EPSG, "semi-minor axis") //epsg does not specifically define this parameter |
| 299 |
}, |
| 300 |
6378388.0*(1-1/297.0), 6378388.0*(1-1/297.0), 6378388.0*(1-1/297.0), SI.METER); |
| 301 |
|
| 302 |
/**
|
| 303 |
* The operation parameter descriptor for the {@link #centralMeridian centralMeridian}
|
| 304 |
* parameter value. Valid values range is from -180 to 180. Default value is 0.
|
| 305 |
*/
|
| 306 |
public static final ParameterDescriptor CENTRAL_MERIDIAN_LOCAL = createDescriptor( |
| 307 |
new NamedIdentifier[] { |
| 308 |
new NamedIdentifier(CitationImpl.OGC, "central_meridian"), |
| 309 |
new NamedIdentifier(CitationImpl.EPSG, "Longitude of natural origin"), |
| 310 |
new NamedIdentifier(CitationImpl.EPSG, "Longitude of false origin"), |
| 311 |
new NamedIdentifier(CitationImpl.ESRI, "Longitude_Of_Origin"), |
| 312 |
new NamedIdentifier(CitationImpl.ESRI, "Longitude_Of_Center"), //ESRI uses this in orthographic (not to be confused with Longitude_Of_Center in oblique mercator) |
| 313 |
new NamedIdentifier(CitationImpl.GEOTIFF, "NatOriginLong") |
| 314 |
}, |
| 315 |
173.0, 173.0, 173.0, NonSI.DEGREE_ANGLE); |
| 316 |
|
| 317 |
/**
|
| 318 |
* The operation parameter descriptor for the {@link #latitudeOfOrigin latitudeOfOrigin}
|
| 319 |
* parameter value. Valid values range is from -90 to 90. Default value is 0.
|
| 320 |
*/
|
| 321 |
public static final ParameterDescriptor LATITUDE_OF_ORIGIN_LOCAL = createDescriptor( |
| 322 |
new NamedIdentifier[] { |
| 323 |
new NamedIdentifier(CitationImpl.OGC, "latitude_of_origin"), |
| 324 |
new NamedIdentifier(CitationImpl.EPSG, "Latitude of false origin"), |
| 325 |
new NamedIdentifier(CitationImpl.EPSG, "Latitude of natural origin"), |
| 326 |
new NamedIdentifier(CitationImpl.ESRI, "Latitude_Of_Center"), //ESRI uses this in orthographic |
| 327 |
new NamedIdentifier(CitationImpl.GEOTIFF, "NatOriginLat") |
| 328 |
}, |
| 329 |
-41.0, -41.0, -41.0, NonSI.DEGREE_ANGLE); |
| 330 |
|
| 331 |
/**
|
| 332 |
* The operation parameter descriptor for the {@link #falseEasting falseEasting}
|
| 333 |
* parameter value. Valid values range is unrestricted. Default value is 0.
|
| 334 |
*/
|
| 335 |
public static final ParameterDescriptor FALSE_EASTING_LOCAL = createDescriptor( |
| 336 |
new NamedIdentifier[] { |
| 337 |
new NamedIdentifier(CitationImpl.OGC, "false_easting"), |
| 338 |
new NamedIdentifier(CitationImpl.EPSG, "False easting"), |
| 339 |
new NamedIdentifier(CitationImpl.EPSG, "Easting at false origin"), |
| 340 |
new NamedIdentifier(CitationImpl.GEOTIFF, "FalseEasting") |
| 341 |
}, |
| 342 |
2510000.0, 2510000.0, 2510000.0, SI.METER); |
| 343 |
|
| 344 |
/**
|
| 345 |
* The operation parameter descriptor for the {@link #falseNorthing falseNorthing}
|
| 346 |
* parameter value. Valid values range is unrestricted. Default value is 0.
|
| 347 |
*/
|
| 348 |
public static final ParameterDescriptor FALSE_NORTHING_LOCAL = createDescriptor( |
| 349 |
new NamedIdentifier[] { |
| 350 |
new NamedIdentifier(CitationImpl.OGC, "false_northing"), |
| 351 |
new NamedIdentifier(CitationImpl.EPSG, "False northing"), |
| 352 |
new NamedIdentifier(CitationImpl.EPSG, "Northing at false origin"), |
| 353 |
new NamedIdentifier(CitationImpl.GEOTIFF, "FalseNorthing") |
| 354 |
}, |
| 355 |
6023150.0, 6023150.0, 6023150.0, SI.METER); |
| 356 |
|
| 357 |
/**
|
| 358 |
* The parameters group.
|
| 359 |
*/
|
| 360 |
static final ParameterDescriptorGroup PARAMETERS = createDescriptorGroup(new Identifier[] { |
| 361 |
new NamedIdentifier(CitationImpl.OGC, "New_Zealand_Map_Grid"), |
| 362 |
new NamedIdentifier(CitationImpl.EPSG, "New Zealand Map Grid"), |
| 363 |
new NamedIdentifier(CitationImpl.EPSG, "27200"), |
| 364 |
new NamedIdentifier(new CitationImpl("IDR"), "IDR") |
| 365 |
}, |
| 366 |
new ParameterDescriptor[] { |
| 367 |
SEMI_MAJOR_LOCAL, SEMI_MINOR_LOCAL, |
| 368 |
LATITUDE_OF_ORIGIN_LOCAL, CENTRAL_MERIDIAN_LOCAL, |
| 369 |
FALSE_EASTING_LOCAL, FALSE_NORTHING_LOCAL |
| 370 |
}); |
| 371 |
/*new ParameterDescriptor[] {
|
| 372 |
new ModifiedParameterDescriptor(SEMI_MAJOR_L, 6378388.0),
|
| 373 |
new ModifiedParameterDescriptor(SEMI_MINOR, 6378388.0*(1-1/297.0)),
|
| 374 |
new ModifiedParameterDescriptor(LATITUDE_OF_ORIGIN, -41.0),
|
| 375 |
new ModifiedParameterDescriptor(CENTRAL_MERIDIAN, 173.0),
|
| 376 |
new ModifiedParameterDescriptor(FALSE_EASTING, 2510000.0),
|
| 377 |
new ModifiedParameterDescriptor(FALSE_NORTHING, 6023150.0)*/
|
| 378 |
|
| 379 |
|
| 380 |
/**
|
| 381 |
* Constructs a new provider.
|
| 382 |
*/
|
| 383 |
public Provider() { |
| 384 |
super(PARAMETERS);
|
| 385 |
} |
| 386 |
|
| 387 |
/**
|
| 388 |
* Creates a transform from the specified group of parameter values. This method doesn't
|
| 389 |
* check for the spherical case, since the New Zealand Map Grid projection is used with
|
| 390 |
* the International 1924 ellipsoid.
|
| 391 |
*
|
| 392 |
* @param parameters The group of parameter values.
|
| 393 |
* @return The created math transform.
|
| 394 |
* @throws ParameterNotFoundException if a required parameter was not found.
|
| 395 |
*/
|
| 396 |
public MathTransform createMathTransform(final ParameterValueGroup parameters) |
| 397 |
throws ParameterNotFoundException
|
| 398 |
{
|
| 399 |
return new IdrNewZealandMapGrid(parameters); |
| 400 |
} |
| 401 |
} |
| 402 |
} |