use of org.opengis.referencing.operation.NoninvertibleTransformException in project sis by apache.
the class Shapes2D method transform.
/**
* Transforms a rectangular envelope using the given coordinate operation.
* The transformation is only approximative: the returned envelope may be bigger
* than the smallest possible bounding box, but should not be smaller in most cases.
*
* <p>This method can handle the case where the rectangle contains the North or South pole,
* or when it cross the ±180° longitude.</p>
*
* @param operation the operation to use. Source and target dimension must be 2.
* @param envelope the rectangle to transform (may be {@code null}).
* @param destination the destination rectangle (may be {@code envelope}).
* If {@code null}, a new rectangle will be created and returned.
* @return {@code destination}, or a new rectangle if {@code destination} was non-null and {@code envelope} was null.
* @throws TransformException if a transform failed.
*
* @see #transform(MathTransform2D, Rectangle2D, Rectangle2D)
* @see Envelopes#transform(CoordinateOperation, Envelope)
*/
@SuppressWarnings("null")
public static Rectangle2D transform(final CoordinateOperation operation, final Rectangle2D envelope, Rectangle2D destination) throws TransformException {
ArgumentChecks.ensureNonNull("operation", operation);
if (envelope == null) {
return null;
}
final MathTransform transform = operation.getMathTransform();
if (!(transform instanceof MathTransform2D)) {
throw new MismatchedDimensionException(Errors.format(Errors.Keys.IllegalPropertyValueClass_3, "transform", MathTransform2D.class, MathTransform.class));
}
MathTransform2D mt = (MathTransform2D) transform;
final double[] center = new double[2];
destination = transform(mt, envelope, destination, center);
/*
* If the source envelope crosses the expected range of valid coordinates, also projects
* the range bounds as a safety. See the comments in transform(Envelope, ...).
*/
final CoordinateReferenceSystem sourceCRS = operation.getSourceCRS();
if (sourceCRS != null) {
final CoordinateSystem cs = sourceCRS.getCoordinateSystem();
if (cs != null && cs.getDimension() == 2) {
// Paranoiac check.
CoordinateSystemAxis axis = cs.getAxis(0);
double min = envelope.getMinX();
double max = envelope.getMaxX();
Point2D.Double pt = null;
for (int i = 0; i < 4; i++) {
if (i == 2) {
axis = cs.getAxis(1);
min = envelope.getMinY();
max = envelope.getMaxY();
}
final double v = (i & 1) == 0 ? axis.getMinimumValue() : axis.getMaximumValue();
if (!(v > min && v < max)) {
continue;
}
if (pt == null) {
pt = new Point2D.Double();
}
if ((i & 2) == 0) {
pt.x = v;
pt.y = envelope.getCenterY();
} else {
pt.x = envelope.getCenterX();
pt.y = v;
}
destination.add(mt.transform(pt, pt));
}
}
}
/*
* Now take the target CRS in account.
*/
final CoordinateReferenceSystem targetCRS = operation.getTargetCRS();
if (targetCRS == null) {
return destination;
}
final CoordinateSystem targetCS = targetCRS.getCoordinateSystem();
if (targetCS == null || targetCS.getDimension() != 2) {
// It should be an error, but we keep this method tolerant.
return destination;
}
/*
* Checks for singularity points. See the Envelopes.transform(CoordinateOperation, Envelope)
* method for comments about the algorithm. The code below is the same algorithm adapted for
* the 2D case and the related objects (Point2D, Rectangle2D, etc.).
*
* The 'border' variable in the loop below is used in order to compress 2 dimensions
* and 2 extremums in a single loop, in this order: (xmin, xmax, ymin, ymax).
*/
TransformException warning = null;
Point2D sourcePt = null;
Point2D targetPt = null;
// A bitmask for each (dimension, extremum) pairs.
int includedBoundsValue = 0;
for (int border = 0; border < 4; border++) {
// 2 dimensions and 2 extremums compacted in a flag.
// The dimension index being examined.
final int dimension = border >>> 1;
final CoordinateSystemAxis axis = targetCS.getAxis(dimension);
if (axis == null) {
// Should never be null, but check as a paranoiac safety.
continue;
}
final double extremum = (border & 1) == 0 ? axis.getMinimumValue() : axis.getMaximumValue();
if (Double.isInfinite(extremum) || Double.isNaN(extremum)) {
continue;
}
if (targetPt == null) {
try {
mt = mt.inverse();
} catch (NoninvertibleTransformException exception) {
Envelopes.recoverableException(Shapes2D.class, exception);
return destination;
}
targetPt = new Point2D.Double();
}
switch(dimension) {
case 0:
targetPt.setLocation(extremum, center[1]);
break;
case 1:
targetPt.setLocation(center[0], extremum);
break;
default:
throw new AssertionError(border);
}
try {
sourcePt = mt.transform(targetPt, sourcePt);
} catch (TransformException exception) {
if (warning == null) {
warning = exception;
} else {
warning.addSuppressed(exception);
}
continue;
}
if (envelope.contains(sourcePt)) {
destination.add(targetPt);
includedBoundsValue |= (1 << border);
}
}
/*
* Iterate over all dimensions of type "WRAPAROUND" for which minimal or maximal axis
* values have not yet been included in the envelope. We could inline this check inside
* the above loop, but we don't in order to have a chance to exclude the dimensions for
* which the point have already been added.
*
* See transform(CoordinateOperation, Envelope) for more comments about the algorithm.
*/
if (includedBoundsValue != 0) {
/*
* Bits mask transformation:
* 1) Swaps the two dimensions (YyXx → XxYy)
* 2) Insert a space between each bits (XxYy → X.x.Y.y.)
* 3) Fill the space with duplicated values (X.x.Y.y. → XXxxYYyy)
*
* In terms of bit positions 1,2,4,8 (not bit values), we have:
*
* 8421 → 22881144
* i.e. (ymax, ymin, xmax, xmin) → (xmax², ymax², xmin², ymin²)
*
* Now look at the last part: (xmin², ymin²). The next step is to perform a bitwise
* AND operation in order to have only both of the following conditions:
*
* Borders not yet added to the envelope: ~(ymax, ymin, xmax, xmin)
* Borders in which a singularity exists: (xmin, xmin, ymin, ymin)
*
* The same operation is repeated on the next 4 bits for (xmax, xmax, ymax, ymax).
*/
int toTest = ((includedBoundsValue & 1) << 3) | ((includedBoundsValue & 4) >>> 1) | ((includedBoundsValue & 2) << 6) | ((includedBoundsValue & 8) << 2);
// Duplicate the bit values.
toTest |= (toTest >>> 1);
toTest &= ~(includedBoundsValue | (includedBoundsValue << 4));
/*
* Forget any axes that are not of kind "WRAPAROUND". Then get the final
* bit pattern indicating which points to test. Iterate over that bits.
*/
if ((toTest & 0x33333333) != 0 && !CoordinateOperations.isWrapAround(targetCS.getAxis(0)))
toTest &= 0xCCCCCCCC;
if ((toTest & 0xCCCCCCCC) != 0 && !CoordinateOperations.isWrapAround(targetCS.getAxis(1)))
toTest &= 0x33333333;
while (toTest != 0) {
final int border = Integer.numberOfTrailingZeros(toTest);
final int bitMask = 1 << border;
// Clear now the bit, for the next iteration.
toTest &= ~bitMask;
final int dimensionToAdd = (border >>> 1) & 1;
final CoordinateSystemAxis toAdd = targetCS.getAxis(dimensionToAdd);
final CoordinateSystemAxis added = targetCS.getAxis(dimensionToAdd ^ 1);
double x = (border & 1) == 0 ? toAdd.getMinimumValue() : toAdd.getMaximumValue();
double y = (border & 4) == 0 ? added.getMinimumValue() : added.getMaximumValue();
if (dimensionToAdd != 0) {
final double t = x;
x = y;
y = t;
}
targetPt.setLocation(x, y);
try {
sourcePt = mt.transform(targetPt, sourcePt);
} catch (TransformException exception) {
if (warning == null) {
warning = exception;
} else {
warning.addSuppressed(exception);
}
continue;
}
if (envelope.contains(sourcePt)) {
destination.add(targetPt);
}
}
}
/*
* At this point we finished envelope transformation. Verify if some ordinates 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. In such case we
* set the rectangle to the full range (we do not use the mechanism documented in Envelope2D) because most
* Rectangle2D implementations do not support spanning the anti-meridian. This results in larger rectangle
* than what would be possible with GeneralEnvelope or Envelope2D, but we try to limit the situation where
* this expansion is applied.
*/
final Set<Integer> wrapAroundChanges;
if (operation instanceof AbstractCoordinateOperation) {
wrapAroundChanges = ((AbstractCoordinateOperation) operation).getWrapAroundChanges();
} else {
wrapAroundChanges = CoordinateOperations.wrapAroundChanges(sourceCRS, targetCS);
}
for (int dim : wrapAroundChanges) {
// Empty in the vast majority of cases.
final CoordinateSystemAxis axis = targetCS.getAxis(dim);
final double minimum = axis.getMinimumValue();
final double maximum = axis.getMaximumValue();
final double o1, o2;
if (dim == 0) {
o1 = destination.getMinX();
o2 = destination.getMaxX();
} else {
o1 = destination.getMinY();
o2 = destination.getMaxY();
}
if (o1 < minimum || o2 > maximum) {
final double span = maximum - minimum;
if (dim == 0) {
destination.setRect(minimum, destination.getY(), span, destination.getHeight());
} else {
destination.setRect(destination.getX(), minimum, destination.getWidth(), span);
}
}
}
if (warning != null) {
Envelopes.recoverableException(Shapes2D.class, warning);
}
return destination;
}
use of org.opengis.referencing.operation.NoninvertibleTransformException in project sis by apache.
the class FranceGeocentricInterpolation method createMathTransform.
/**
* Creates a transform from the specified group of parameter values.
* This method creates the transform from <em>target</em> to <em>source</em>
* (which is the direction that use the interpolation grid directly without iteration),
* then inverts the transform.
*
* @param factory the factory to use if this constructor needs to create other math transforms.
* @param values the group of parameter values.
* @return the created math transform.
* @throws ParameterNotFoundException if a required parameter was not found.
* @throws FactoryException if an error occurred while loading the grid.
*/
@Override
public MathTransform createMathTransform(final MathTransformFactory factory, final ParameterValueGroup values) throws ParameterNotFoundException, FactoryException {
boolean withHeights = false;
final Parameters pg = Parameters.castOrWrap(values);
final Integer dim = pg.getValue(Molodensky.DIMENSION);
if (dim != null)
switch(dim) {
case 2:
break;
case 3:
withHeights = true;
break;
default:
throw new InvalidParameterValueException(Errors.format(Errors.Keys.IllegalArgumentValue_2, "dim", dim), "dim", dim);
}
final Path file = pg.getMandatoryValue(FILE);
final DatumShiftGridFile<Angle, Length> grid = getOrLoad(file, isRecognized(file) ? new double[] { TX, TY, TZ } : null, PRECISION);
MathTransform tr = createGeodeticTransformation(factory, createEllipsoid(pg, Molodensky.TGT_SEMI_MAJOR, Molodensky.TGT_SEMI_MINOR, // GRS 1980 ellipsoid
CommonCRS.ETRS89.ellipsoid()), createEllipsoid(pg, Molodensky.SRC_SEMI_MAJOR, Molodensky.SRC_SEMI_MINOR, // Clarke 1880 (IGN) ellipsoid
null), withHeights, grid);
try {
tr = tr.inverse();
} catch (NoninvertibleTransformException e) {
// Should never happen.
throw new FactoryException(e);
}
return tr;
}
use of org.opengis.referencing.operation.NoninvertibleTransformException in project sis by apache.
the class FranceGeocentricInterpolation method getOrLoad.
/**
* Returns the grid of the given name. This method returns the cached instance if it still exists,
* or load the grid otherwise.
*
* @param file name of the datum shift grid file to load.
* @param averages an "average" value for the offset in each dimension, or {@code null} if unknown.
* @param scale the factor by which to multiply each compressed value before to add to the average value.
*/
@SuppressWarnings("null")
static DatumShiftGridFile<Angle, Length> getOrLoad(final Path file, final double[] averages, final double scale) throws FactoryException {
final Path resolved = DataDirectory.DATUM_CHANGES.resolve(file).toAbsolutePath();
DatumShiftGridFile<?, ?> grid = DatumShiftGridFile.CACHE.peek(resolved);
if (grid == null) {
final Cache.Handler<DatumShiftGridFile<?, ?>> handler = DatumShiftGridFile.CACHE.lock(resolved);
try {
grid = handler.peek();
if (grid == null) {
try (BufferedReader in = Files.newBufferedReader(resolved)) {
DatumShiftGridLoader.log(FranceGeocentricInterpolation.class, file);
final DatumShiftGridFile.Float<Angle, Length> g = load(in, file);
grid = DatumShiftGridCompressed.compress(g, averages, scale);
} catch (IOException | NoninvertibleTransformException | RuntimeException e) {
// NumberFormatException, ArithmeticException, NoSuchElementException, possibly other.
throw DatumShiftGridLoader.canNotLoad(HEADER, file, e);
}
grid = grid.useSharedData();
}
} finally {
handler.putAndUnlock(grid);
}
}
return grid.castTo(Angle.class, Length.class);
}
use of org.opengis.referencing.operation.NoninvertibleTransformException in project sis by apache.
the class GeocentricAffineBetweenGeographic method createMathTransform.
/**
* Creates a math transform from the specified group of parameter values.
* This method wraps the affine operation into Geographic/Geocentric conversions.
*
* @param factory the factory to use for creating concatenated transforms.
* @param values the group of parameter values.
* @return the created math transform.
* @throws FactoryException if a transform can not be created.
*/
@Override
public MathTransform createMathTransform(final MathTransformFactory factory, final ParameterValueGroup values) throws FactoryException {
final Parameters pv = Parameters.castOrWrap(values);
final MathTransform affine = super.createMathTransform(factory, pv);
/*
* Create a "Geographic to Geocentric" conversion with ellipsoid axis length units converted to metres
* (the unit implied by SRC_SEMI_MAJOR) because it is the unit of Bursa-Wolf parameters that we created above.
*/
MathTransform toGeocentric = EllipsoidToCentricTransform.createGeodeticConversion(factory, pv.doubleValue(SRC_SEMI_MAJOR), pv.doubleValue(SRC_SEMI_MINOR), Units.METRE, getSourceDimensions() >= 3, EllipsoidToCentricTransform.TargetType.CARTESIAN);
/*
* Create a "Geocentric to Geographic" conversion with ellipsoid axis length units converted to metres
* because this is the unit of the Geocentric CRS used above.
*/
MathTransform toGeographic = EllipsoidToCentricTransform.createGeodeticConversion(factory, pv.doubleValue(TGT_SEMI_MAJOR), pv.doubleValue(TGT_SEMI_MINOR), Units.METRE, getTargetDimensions() >= 3, EllipsoidToCentricTransform.TargetType.CARTESIAN);
try {
toGeographic = toGeographic.inverse();
} catch (NoninvertibleTransformException e) {
// Should never happen with SIS implementation.
throw new FactoryException(e);
}
/*
* The Geocentric → Affine → Geographic chain.
*/
return factory.createConcatenatedTransform(toGeocentric, factory.createConcatenatedTransform(affine, toGeographic));
}
use of org.opengis.referencing.operation.NoninvertibleTransformException in project sis by apache.
the class NTv2 method getOrLoad.
/**
* Returns the grid of the given name. This method returns the cached instance if it still exists,
* or load the grid otherwise.
*
* @param file name of the datum shift grid file to load.
*/
@SuppressWarnings("null")
static DatumShiftGridFile<Angle, Angle> getOrLoad(final Path file) throws FactoryException {
final Path resolved = DataDirectory.DATUM_CHANGES.resolve(file).toAbsolutePath();
DatumShiftGridFile<?, ?> grid = DatumShiftGridFile.CACHE.peek(resolved);
if (grid == null) {
final Cache.Handler<DatumShiftGridFile<?, ?>> handler = DatumShiftGridFile.CACHE.lock(resolved);
try {
grid = handler.peek();
if (grid == null) {
try (ReadableByteChannel in = Files.newByteChannel(resolved)) {
DatumShiftGridLoader.log(NTv2.class, file);
final Loader loader = new Loader(in, file);
grid = loader.readGrid();
loader.reportWarnings();
} catch (IOException | NoninvertibleTransformException | RuntimeException e) {
throw DatumShiftGridLoader.canNotLoad("NTv2", file, e);
}
grid = grid.useSharedData();
}
} finally {
handler.putAndUnlock(grid);
}
}
return grid.castTo(Angle.class, Angle.class);
}
Aggregations