in endorsed/src/org.apache.sis.referencing/main/org/apache/sis/geometry/Envelopes.java [614:947]
public static GeneralEnvelope transform(final CoordinateOperation operation, Envelope envelope)
throws TransformException
{
ArgumentChecks.ensureNonNull("operation", operation);
if (envelope == null) {
return null;
}
boolean isOperationComplete = true;
final CoordinateReferenceSystem sourceCRS = operation.getSourceCRS();
if (sourceCRS != null) {
final CoordinateReferenceSystem crs = envelope.getCoordinateReferenceSystem();
if (crs != null && !Utilities.equalsIgnoreMetadata(crs, sourceCRS)) {
/*
* Argument-check: the envelope CRS seems inconsistent with the given operation.
* However, we need to push the check a little bit further, since 3D-GeographicCRS
* are considered not equal to CompoundCRS[2D-GeographicCRS + ellipsoidal height].
* Checking for identity MathTransform is a more powerfull (but more costly) check.
* Since we have the MathTransform, perform an opportunist envelope transform if it
* happen to be required.
*/
final MathTransform mt;
try {
mt = CRS.findOperation(crs, sourceCRS, null).getMathTransform();
} catch (FactoryException e) {
throw new TransformException(Errors.format(Errors.Keys.CanNotTransformEnvelope), e);
}
if (!mt.isIdentity()) {
isOperationComplete = false;
envelope = transform(mt, envelope);
}
}
}
final MathTransform mt = operation.getMathTransform();
final double[] centerPt = new double[mt.getTargetDimensions()];
final GeneralEnvelope transformed = transform(mt, envelope, centerPt, null);
/*
* If the source envelope crosses the expected range of valid coordinates, also projects
* the range bounds as a safety. Example: if the source envelope goes from 150 to 200°E,
* some map projections will interpret 200° as if it was -160°, and consequently produce
* an envelope which do not include the 180°W extremum. We will add those extremum points
* explicitly as a safety. It may leads to bigger than necessary target envelope, but the
* contract is to include at least the source envelope, not to return the smallest one.
*/
if (sourceCRS != null) {
final CoordinateSystem cs = sourceCRS.getCoordinateSystem();
if (cs != null) { // Should never be null, but check as a paranoiac safety.
DirectPosition sourcePt = null;
DirectPosition targetPt = null;
final int dimension = cs.getDimension();
for (int i=0; i<dimension; i++) {
final CoordinateSystemAxis axis = cs.getAxis(i);
if (axis == null) { // Should never be null, but check as a paranoiac safety.
continue;
}
final double min = envelope.getMinimum(i);
final double max = envelope.getMaximum(i);
final double v1 = axis.getMinimumValue();
final double v2 = axis.getMaximumValue();
final boolean b1 = (v1 > min && v1 < max);
final boolean b2 = (v2 > min && v2 < max);
if (!b1 && !b2) {
continue;
}
if (sourcePt == null) {
sourcePt = new GeneralDirectPosition(dimension);
for (int j=0; j<dimension; j++) {
sourcePt.setOrdinate(j, envelope.getMedian(j));
}
}
if (b1) {
sourcePt.setOrdinate(i, v1);
transformed.add(targetPt = mt.transform(sourcePt, targetPt));
}
if (b2) {
sourcePt.setOrdinate(i, v2);
transformed.add(targetPt = mt.transform(sourcePt, targetPt));
}
sourcePt.setOrdinate(i, envelope.getMedian(i));
}
}
}
/*
* Now takes the target CRS in account...
*/
final CoordinateReferenceSystem targetCRS = operation.getTargetCRS();
if (targetCRS == null) {
return transformed;
}
transformed.setCoordinateReferenceSystem(targetCRS);
final CoordinateSystem targetCS = targetCRS.getCoordinateSystem();
if (targetCS == null) {
// It should be an error, but we keep this method tolerant.
return transformed;
}
/*
* Checks for singularity points. For example, the south pole is a singularity point in
* geographic CRS because it is located at the maximal value allowed by one particular
* axis, namely latitude. This point is not a singularity in the stereographic projection,
* because axes extends toward infinity in all directions (mathematically) and because the
* South pole has nothing special apart being the origin (0,0).
*
* Algorithm:
*
* 1) Inspect the target axis, looking if there is any bounds. If bounds are found, get
* the coordinates of singularity points and project them from target to source CRS.
*
* Example: If the transformed envelope above is (80 … 85°S, 10 … 50°W), and if the
* latitude in the target CRS is bounded at 90°S, then project (90°S, 30°W)
* to the source CRS. Note that the longitude is set to the center of
* the envelope longitude range (more on this below).
*
* 2) If the singularity point computed above is inside the source envelope,
* add that point to the target (transformed) envelope.
*
* 3) For each singularity point found at step #2, check if there is a wraparound axis
* in a dimension other than the dimension that was set to an axis bounds value.
* If such wraparound axis is found, test for inclusion of the same point as the
* point tested at step #1, except for the coordinate of the wraparound axis which
* is set to the minimal and maximal values (reminder: step #1 used the center value).
*
* Example: If the above steps found that the point (90°S, 30°W) need to be included,
* then this step #3 will also test the points (90°S, 180°W) and (90°S, 180°E).
*
* NOTE: we test (-180°, centerY), (180°, centerY), (centerX, -90°) and (centerX, 90°)
* at step #1 before to test (-180°, -90°), (180°, -90°), (-180°, 90°) and (180°, 90°)
* at step #3 because the latter may not be supported by every projections. For example
* if the target envelope is located between 20°N and 40°N, then a Mercator projection
* may fail to transform the (-180°, 90°) coordinate while the (-180°, 30°) coordinate
* is a valid point.
*/
MathTransform inverse = null;
TransformException warning = null;
AbstractEnvelope sourceBox = null;
DirectPosition sourcePt = null;
DirectPosition targetPt = null;
DirectPosition revertPt = null;
long includedMinValue = 0; // A bitmask for each dimension.
long includedMaxValue = 0;
long isWrapAroundAxis = 0;
final int dimension = targetCS.getDimension();
poles: for (int i=0; i<dimension; i++) {
final long dimensionBitMask = Numerics.bitmask(i);
final CoordinateSystemAxis axis = targetCS.getAxis(i);
if (axis == null) { // Should never be null, but check as a paranoiac safety.
continue;
}
boolean testMax = false; // Tells if we are testing the minimal or maximal value.
do {
final double extremum = testMax ? axis.getMaximumValue() : axis.getMinimumValue();
if (!Double.isFinite(extremum)) {
/*
* The axis is unbounded. It should always be the case when the target CRS is
* a map projection, in which case this loop will finish soon and this method
* will do nothing more (no object instantiated, no MathTransform inversed...)
*/
continue;
}
if (inverse == null) {
try {
inverse = mt.inverse();
} catch (NoninvertibleTransformException exception) {
/*
* If the transform is non-invertible, this "poles" loop cannot do anything.
* This is not a fatal error because the envelope has already be transformed
* by the code before this loop. We lost the check below for singularity points,
* but it makes no difference in the common case where all those points would
* have been outside the source envelope anyway.
*
* Note that this exception is normal if the number of target dimension is smaller
* than the number of source dimension, because the math transform cannot guess
* coordinates in the lost dimensions. So we do not log any warning in this case.
*/
if (dimension >= mt.getSourceDimensions()) {
warning = exception;
}
break poles;
}
targetPt = new GeneralDirectPosition(centerPt.clone());
sourceBox = AbstractEnvelope.castOrCopy(envelope);
}
targetPt.setOrdinate(i, extremum);
try {
sourcePt = inverse.transform(targetPt, sourcePt);
if (sourceBox.contains(sourcePt)) {
/*
* The point is inside the source envelope and consequently should be added in the
* transformed envelope. However, there is a possible confusion: if the axis that we
* tested is a wraparound axis, then (for example) +180° and -180° of longitude may
* be the same point in source CRS, despite being 2 very different points in target
* CRS. We do yet another projection in opposite direction for checking if we really
* have the point that we wanted to test.
*/
if (CoordinateOperations.isWrapAround(axis)) {
revertPt = mt.transform(sourcePt, revertPt);
final double delta = Math.abs(revertPt.getOrdinate(i) - extremum);
if (!(delta < SPAN_FRACTION_AS_BOUND * (axis.getMaximumValue() - axis.getMinimumValue()))) {
continue;
}
}
transformed.add(targetPt);
if (testMax) includedMaxValue |= dimensionBitMask;
else includedMinValue |= dimensionBitMask;
}
} catch (TransformException exception) {
/*
* This exception may be normal. For example if may occur when projecting
* the latitude extremums with a cylindrical Mercator projection. Do not
* log any message (unless logging level is fine) and try the other points.
*/
if (warning == null) {
warning = exception;
} else {
warning.addSuppressed(exception);
}
}
} while ((testMax = !testMax) == true);
/*
* Keep trace of axes of kind WRAPAROUND, except if the two extremum values of that
* axis have been included in the envelope (in which case the next step after this
* loop does not need to be executed for this axis).
*/
if ((includedMinValue & includedMaxValue & dimensionBitMask) == 0 && CoordinateOperations.isWrapAround(axis)) {
isWrapAroundAxis |= dimensionBitMask;
}
// Restore `targetPt` to its initial state, which is equal to `centerPt`.
if (targetPt != null) {
targetPt.setOrdinate(i, centerPt[i]);
}
}
/*
* Step #3 described in the above "Algorithm" section: iterate over all dimensions
* of type "WRAPAROUND" for which minimal or maximal axis values have not yet been
* included in the envelope. The set of axes is specified by a bitmask computed in
* the above loop. We examine only the points that have not already been included
* in the envelope.
*/
final long includedBoundsValue = (includedMinValue | includedMaxValue);
if (includedBoundsValue != 0) {
while (isWrapAroundAxis != 0) {
final int wrapAroundDimension = Long.numberOfTrailingZeros(isWrapAroundAxis);
final long dimensionBitMask = Numerics.bitmask(wrapAroundDimension);
isWrapAroundAxis &= ~dimensionBitMask; // Clear now the bit, for the next iteration.
final CoordinateSystemAxis wrapAroundAxis = targetCS.getAxis(wrapAroundDimension);
final double min = wrapAroundAxis.getMinimumValue();
final double max = wrapAroundAxis.getMaximumValue();
/*
* Iterate over all axes for which a singularity point has been previously found,
* excluding the "wrap around axis" currently under consideration.
*/
for (long am=(includedBoundsValue & ~dimensionBitMask), bm; am != 0; am &= ~bm) {
bm = Long.lowestOneBit(am);
final int axisIndex = Long.numberOfTrailingZeros(bm);
final CoordinateSystemAxis axis = targetCS.getAxis(axisIndex);
/*
* switch (c) {
* case 0: targetPt = (..., singularityMin, ..., wrapAroundMin, ...)
* case 1: targetPt = (..., singularityMin, ..., wrapAroundMax, ...)
* case 2: targetPt = (..., singularityMax, ..., wrapAroundMin, ...)
* case 3: targetPt = (..., singularityMax, ..., wrapAroundMax, ...)
* }
*/
for (int c=0; c<4; c++) {
/*
* Set the coordinate value along the axis having the singularity point
* (cases c=0 and c=2). If the envelope did not included that point,
* then skip completely this case and the next one, i.e. skip c={0,1}
* or skip c={2,3}.
*/
final double value;
if ((c & 1) == 0) { // `true` if we are testing "wrapAroundMin".
if (((c == 0 ? includedMinValue : includedMaxValue) & bm) == 0) {
c++; // Skip also the case for "wrapAroundMax".
continue;
}
targetPt.setOrdinate(axisIndex, (c == 0) ? axis.getMinimumValue() : axis.getMaximumValue());
value = min;
} else {
value = max;
}
targetPt.setOrdinate(wrapAroundDimension, value);
try {
sourcePt = inverse.transform(targetPt, sourcePt);
if (sourceBox.contains(sourcePt)) {
/*
* The `value` limit is typically the 180°E or 180°W longitude value. We could test
* its validity as below (see similar code in other loop above for explanation):
*
* revertPt = mt.transform(sourcePt, revertPt);
* final double delta = Math.abs(revertPt.getCoordinate(wrapAroundDimension) - value);
* if (delta < SPAN_FRACTION_AS_BOUND * (max - min)) {
* transformed.add(targetPt);
* }
*
* But we don't because this block is executed when another coordinate is at its bounds,
* and that other coordinate is usually the latitude at 90°S or 90°N. In such case, all
* longitude values are the same point on Earth (a pole) and can be anything. Comparing
* `revertPt` coordinate with `value` is meaningless. In order to keep above check, we
* would need to determine if `targetPt` is at a pole. But we have fully reliable way.
*/
transformed.add(targetPt);
}
} catch (TransformException exception) {
if (warning == null) {
warning = exception;
} else {
warning.addSuppressed(exception);
}
}
}
targetPt.setOrdinate(axisIndex, centerPt[axisIndex]);
}
targetPt.setOrdinate(wrapAroundDimension, centerPt[wrapAroundDimension]);
}
}
/*
* At this point we finished envelope transformation. Verify if some coordinates need to be "wrapped around"
* as a result of the coordinate operation. This is usually the longitude axis where the source CRS uses
* the [-180 … +180]° range and the target CRS uses the [0 … 360]° range, or the converse. We do not wrap
* around if the source and target axes use the same range (e.g. the longitude stay [-180 … +180]°) in order
* to reduce the risk of discontinuities. If the user really wants unconditional wrap around, (s)he can call
* `GeneralEnvelope.normalize()`.
*/
final Set<Integer> wrapAroundChanges;
if (isOperationComplete && operation instanceof AbstractCoordinateOperation) {
wrapAroundChanges = ((AbstractCoordinateOperation) operation).getWrapAroundChanges();
} else {
wrapAroundChanges = CoordinateOperations.wrapAroundChanges(sourceCRS, targetCS);
}
transformed.normalize(targetCS, 0, wrapAroundChanges.size(), wrapAroundChanges.iterator());
if (warning != null) {
recoverableException(Envelopes.class, warning);
}
return transformed;
}