Search in sources :

Example 6 with Length

use of javax.measure.quantity.Length in project sis by apache.

the class PositionalAccuracyConstant method getLinearAccuracy.

/**
 * Convenience method returning the accuracy in meters for the specified operation.
 * This method tries each of the following procedures and returns the first successful one:
 *
 * <ul>
 *   <li>If at least one {@link QuantitativeResult} is found with a linear unit, then the largest
 *       accuracy estimate is converted to {@linkplain Units#METRE metres} and returned.</li>
 *   <li>Otherwise, if the operation is a {@link Conversion}, then returns 0 since a conversion
 *       is by definition accurate up to rounding errors.</li>
 *   <li>Otherwise, if the operation is a {@link Transformation}, then checks if the datum shift
 *       were applied with the help of Bursa-Wolf parameters. This procedure looks for SIS-specific
 *       {@link #DATUM_SHIFT_APPLIED} and {@link #DATUM_SHIFT_OMITTED DATUM_SHIFT_OMITTED} constants.</li>
 *   <li>Otherwise, if the operation is a {@link ConcatenatedOperation}, returns the sum of the accuracy
 *       of all components. This is a conservative scenario where we assume that errors cumulate linearly.
 *       Note that this is not necessarily the "worst case" scenario since the accuracy could be worst
 *       if the math transforms are highly non-linear.</li>
 * </ul>
 *
 * If the above is modified, please update {@code AbstractCoordinateOperation.getLinearAccuracy()} javadoc.
 *
 * @param  operation  the operation to inspect for accuracy.
 * @return the accuracy estimate (always in meters), or NaN if unknown.
 *
 * @see org.apache.sis.referencing.operation.AbstractCoordinateOperation#getLinearAccuracy()
 */
public static double getLinearAccuracy(final CoordinateOperation operation) {
    double accuracy = Double.NaN;
    final Collection<PositionalAccuracy> accuracies = operation.getCoordinateOperationAccuracy();
    for (final PositionalAccuracy metadata : accuracies) {
        for (final Result result : metadata.getResults()) {
            if (result instanceof QuantitativeResult) {
                final QuantitativeResult quantity = (QuantitativeResult) result;
                final Collection<? extends Record> records = quantity.getValues();
                if (records != null) {
                    final Unit<?> unit = quantity.getValueUnit();
                    if (Units.isLinear(unit)) {
                        final Unit<Length> unitOfLength = unit.asType(Length.class);
                        for (final Record record : records) {
                            for (final Object value : record.getAttributes().values()) {
                                if (value instanceof Number) {
                                    double v = ((Number) value).doubleValue();
                                    v = unitOfLength.getConverterTo(Units.METRE).convert(v);
                                    if (v >= 0 && !(v <= accuracy)) {
                                        // '!' is for replacing the NaN value.
                                        accuracy = v;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    if (Double.isNaN(accuracy)) {
        /*
             * No quantitative (linear) accuracy were found. If the coordinate operation is actually
             * a conversion, the accuracy is up to rounding error (i.e. conceptually 0) by definition.
             */
        if (operation instanceof Conversion) {
            return 0;
        }
        /*
             * If the coordinate operation is actually a transformation, checks if Bursa-Wolf parameters
             * were available for the datum shift. This is SIS-specific. See field javadoc for a rational
             * about the return values chosen.
             */
        if (operation instanceof Transformation) {
            if (accuracies.contains(DATUM_SHIFT_APPLIED)) {
                return DATUM_SHIFT_ACCURACY;
            }
            if (accuracies.contains(DATUM_SHIFT_OMITTED)) {
                return UNKNOWN_ACCURACY;
            }
        }
        /*
             * If the coordinate operation is a compound of other coordinate operations, returns the sum of their accuracy,
             * skipping unknown ones. Making the sum is a conservative approach (not exactly the "worst case" scenario,
             * since it could be worst if the transforms are highly non-linear).
             */
        if (operation instanceof ConcatenatedOperation) {
            for (final CoordinateOperation op : ((ConcatenatedOperation) operation).getOperations()) {
                final double candidate = Math.abs(getLinearAccuracy(op));
                if (!Double.isNaN(candidate)) {
                    if (Double.isNaN(accuracy)) {
                        accuracy = candidate;
                    } else {
                        accuracy += candidate;
                    }
                }
            }
        }
    }
    return accuracy;
}
Also used : Transformation(org.opengis.referencing.operation.Transformation) PositionalAccuracy(org.opengis.metadata.quality.PositionalAccuracy) DefaultAbsoluteExternalPositionalAccuracy(org.apache.sis.metadata.iso.quality.DefaultAbsoluteExternalPositionalAccuracy) CoordinateOperation(org.opengis.referencing.operation.CoordinateOperation) Conversion(org.opengis.referencing.operation.Conversion) Result(org.opengis.metadata.quality.Result) QuantitativeResult(org.opengis.metadata.quality.QuantitativeResult) DefaultConformanceResult(org.apache.sis.metadata.iso.quality.DefaultConformanceResult) QuantitativeResult(org.opengis.metadata.quality.QuantitativeResult) Length(javax.measure.quantity.Length) Record(org.opengis.util.Record) ConcatenatedOperation(org.opengis.referencing.operation.ConcatenatedOperation)

Example 7 with Length

use of javax.measure.quantity.Length 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;
}
Also used : Path(java.nio.file.Path) NoninvertibleTransformException(org.opengis.referencing.operation.NoninvertibleTransformException) Parameters(org.apache.sis.parameter.Parameters) Angle(javax.measure.quantity.Angle) MathTransform(org.opengis.referencing.operation.MathTransform) InvalidParameterValueException(org.opengis.parameter.InvalidParameterValueException) Length(javax.measure.quantity.Length) FactoryException(org.opengis.util.FactoryException)

Example 8 with Length

use of javax.measure.quantity.Length 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);
}
Also used : Path(java.nio.file.Path) IOException(java.io.IOException) NoninvertibleTransformException(org.opengis.referencing.operation.NoninvertibleTransformException) Angle(javax.measure.quantity.Angle) Length(javax.measure.quantity.Length) BufferedReader(java.io.BufferedReader) Cache(org.apache.sis.util.collection.Cache)

Example 9 with Length

use of javax.measure.quantity.Length in project sis by apache.

the class EPSGDataAccess method createEllipsoid.

/**
 * Creates a geometric figure that can be used to describe the approximate shape of the earth.
 * In mathematical terms, it is a surface formed by the rotation of an ellipse about its minor axis.
 *
 * <div class="note"><b>Example:</b>
 * some EPSG codes for ellipsoids are:
 *
 * <table class="sis" summary="EPSG codes examples">
 *   <tr><th>Code</th> <th>Description</th></tr>
 *   <tr><td>7030</td> <td>WGS 84</td></tr>
 *   <tr><td>7034</td> <td>Clarke 1880</td></tr>
 *   <tr><td>7048</td> <td>GRS 1980 Authalic Sphere</td></tr>
 * </table></div>
 *
 * @param  code  value allocated by EPSG.
 * @return the ellipsoid for the given code.
 * @throws NoSuchAuthorityCodeException if the specified {@code code} was not found.
 * @throws FactoryException if the object creation failed for some other reason.
 *
 * @see #createGeodeticDatum(String)
 * @see #createEllipsoidalCS(String)
 * @see org.apache.sis.referencing.datum.DefaultEllipsoid
 */
@Override
public synchronized Ellipsoid createEllipsoid(final String code) throws NoSuchAuthorityCodeException, FactoryException {
    ArgumentChecks.ensureNonNull("code", code);
    Ellipsoid returnValue = null;
    try (ResultSet result = executeQuery("Ellipsoid", "ELLIPSOID_CODE", "ELLIPSOID_NAME", "SELECT ELLIPSOID_CODE," + " ELLIPSOID_NAME," + " SEMI_MAJOR_AXIS," + " INV_FLATTENING," + " SEMI_MINOR_AXIS," + " UOM_CODE," + " REMARKS," + " DEPRECATED" + " FROM [Ellipsoid]" + " WHERE ELLIPSOID_CODE = ?", code)) {
        while (result.next()) {
            /*
                 * One of 'semiMinorAxis' and 'inverseFlattening' values can be NULL in the database.
                 * Consequently, we don't use 'getString(ResultSet, int)' for those parameters because
                 * we do not want to thrown an exception if a NULL value is found.
                 */
            final Integer epsg = getInteger(code, result, 1);
            final String name = getString(code, result, 2);
            final double semiMajorAxis = getDouble(code, result, 3);
            final double inverseFlattening = getOptionalDouble(result, 4);
            final double semiMinorAxis = getOptionalDouble(result, 5);
            final String unitCode = getString(code, result, 6);
            final String remarks = getOptionalString(result, 7);
            final boolean deprecated = getOptionalBoolean(result, 8);
            final Unit<Length> unit = owner.createUnit(unitCode).asType(Length.class);
            final Map<String, Object> properties = createProperties("Ellipsoid", name, epsg, remarks, deprecated);
            final Ellipsoid ellipsoid;
            if (Double.isNaN(inverseFlattening)) {
                if (Double.isNaN(semiMinorAxis)) {
                    // Both are null, which is not allowed.
                    final String column = result.getMetaData().getColumnName(3);
                    throw new FactoryDataException(error().getString(Errors.Keys.NullValueInTable_3, code, column));
                } else {
                    // We only have semiMinorAxis defined. It is OK
                    ellipsoid = owner.datumFactory.createEllipsoid(properties, semiMajorAxis, semiMinorAxis, unit);
                }
            } else {
                if (!Double.isNaN(semiMinorAxis)) {
                    // Both 'inverseFlattening' and 'semiMinorAxis' are defined.
                    // Log a warning and create the ellipsoid using the inverse flattening.
                    final LogRecord record = resources().getLogRecord(Level.WARNING, Resources.Keys.AmbiguousEllipsoid_1, Constants.EPSG + DefaultNameSpace.DEFAULT_SEPARATOR + code);
                    record.setLoggerName(Loggers.CRS_FACTORY);
                    Logging.log(EPSGDataAccess.class, "createEllipsoid", record);
                }
                ellipsoid = owner.datumFactory.createFlattenedSphere(properties, semiMajorAxis, inverseFlattening, unit);
            }
            returnValue = ensureSingleton(ellipsoid, returnValue, code);
        }
    } catch (SQLException exception) {
        throw databaseFailure(Ellipsoid.class, code, exception);
    }
    if (returnValue == null) {
        throw noSuchAuthorityCode(Ellipsoid.class, code);
    }
    return returnValue;
}
Also used : SQLException(java.sql.SQLException) InternationalString(org.opengis.util.InternationalString) SimpleInternationalString(org.apache.sis.util.iso.SimpleInternationalString) FactoryDataException(org.apache.sis.referencing.factory.FactoryDataException) Length(javax.measure.quantity.Length) LogRecord(java.util.logging.LogRecord) ResultSet(java.sql.ResultSet) AbstractIdentifiedObject(org.apache.sis.referencing.AbstractIdentifiedObject) IdentifiedObject(org.opengis.referencing.IdentifiedObject)

Example 10 with Length

use of javax.measure.quantity.Length in project sis by apache.

the class CommonAuthorityFactory method createAuto.

/**
 * Creates a projected CRS from parameters in the {@code AUTO(2)} namespace.
 *
 * @param  code        the user-specified code, used only for error reporting.
 * @param  projection  the projection code (e.g. 42001).
 * @param  isLegacy    {@code true} if the code was found in {@code "AUTO"} or {@code "AUTO1"} namespace.
 * @param  factor      the multiplication factor for the unit of measurement.
 * @param  longitude   a longitude in the desired projection zone.
 * @param  latitude    a latitude in the desired projection zone.
 * @return the projected CRS for the given projection and parameters.
 */
@SuppressWarnings("null")
private ProjectedCRS createAuto(final String code, final int projection, final boolean isLegacy, final double factor, final double longitude, final double latitude) throws FactoryException {
    Boolean isUTM = null;
    String method = null;
    String param = null;
    switch(projection) {
        /*
             * 42001: Universal Transverse Mercator   —   central meridian must be in the center of a UTM zone.
             * 42002: Transverse Mercator             —   like 42001 except that central meridian can be anywhere.
             * 42003: WGS 84 / Auto Orthographic      —   defined by "Central_Meridian" and "Latitude_of_Origin".
             * 42004: WGS 84 / Auto Equirectangular   —   defined by "Central_Meridian" and "Standard_Parallel_1".
             * 42005: WGS 84 / Auto Mollweide         —   defined by "Central_Meridian" only.
             */
        case 42001:
            isUTM = true;
            break;
        case 42002:
            isUTM = (latitude == 0) && (Zoner.UTM.centralMeridian(Zoner.UTM.zone(0, longitude)) == longitude);
            break;
        case 42003:
            method = "Orthographic";
            param = Constants.LATITUDE_OF_ORIGIN;
            break;
        case 42004:
            method = "Equirectangular";
            param = Constants.STANDARD_PARALLEL_1;
            break;
        case 42005:
            method = "Mollweide";
            break;
        default:
            throw noSuchAuthorityCode(String.valueOf(projection), code, null);
    }
    /*
         * For the (Universal) Transverse Mercator case (AUTO:42001 and 42002), we delegate to the CommonCRS
         * enumeration if possible because CommonCRS will itself delegate to the EPSG factory if possible.
         * The Math.signum(latitude) instruction is for preventing "AUTO:42001" to handle the UTM special cases
         * (Norway and Svalbard) or to switch on the Universal Polar Stereographic projection for high latitudes,
         * because the WMS specification does not said that we should.
         */
    final CommonCRS datum = CommonCRS.WGS84;
    // To be set, directly or indirectly, to WGS84.geographic().
    final GeographicCRS baseCRS;
    // Temporary UTM projection, for extracting other properties.
    final ProjectedCRS crs;
    // Coordinate system with (E,N) axes in metres.
    CartesianCS cs;
    try {
        if (isUTM != null && isUTM) {
            crs = datum.universal(Math.signum(latitude), longitude);
            if (factor == (isLegacy ? Constants.EPSG_METRE : 1)) {
                return crs;
            }
            baseCRS = crs.getBaseCRS();
            cs = crs.getCoordinateSystem();
        } else {
            cs = projectedCS;
            if (cs == null) {
                crs = datum.universal(Math.signum(latitude), longitude);
                projectedCS = cs = crs.getCoordinateSystem();
                baseCRS = crs.getBaseCRS();
            } else {
                crs = null;
                baseCRS = datum.geographic();
            }
        }
        /*
             * At this point we got a coordinate system with axes in metres.
             * If the user asked for another unit of measurement, change the axes now.
             */
        Unit<Length> unit;
        if (isLegacy) {
            unit = createUnitFromEPSG(factor).asType(Length.class);
        } else {
            unit = Units.METRE;
            if (factor != 1)
                unit = unit.multiply(factor);
        }
        if (!Units.METRE.equals(unit)) {
            cs = (CartesianCS) CoordinateSystems.replaceLinearUnit(cs, unit);
        }
        /*
             * Set the projection name, operation method and parameters. The parameters for the Transverse Mercator
             * projection are a little bit more tedious to set, so we use a convenience method for that.
             */
        final GeodeticObjectBuilder builder = new GeodeticObjectBuilder();
        if (isUTM != null) {
            if (isUTM && crs != null) {
                builder.addName(crs.getName());
            }
            // else default to the conversion name, which is "Transverse Mercator".
            builder.setTransverseMercator(isUTM ? Zoner.UTM : Zoner.ANY, latitude, longitude);
        } else {
            builder.setConversionMethod(method).addName(PROJECTION_NAMES[projection - FIRST_PROJECTION_CODE]).setParameter(Constants.CENTRAL_MERIDIAN, longitude, Units.DEGREE);
            if (param != null) {
                builder.setParameter(param, latitude, Units.DEGREE);
            }
        }
        return builder.createProjectedCRS(baseCRS, cs);
    } catch (IllegalArgumentException e) {
        throw noSuchAuthorityCode(String.valueOf(projection), code, e);
    }
}
Also used : CartesianCS(org.opengis.referencing.cs.CartesianCS) ProjectedCRS(org.opengis.referencing.crs.ProjectedCRS) Length(javax.measure.quantity.Length) GeodeticObjectBuilder(org.apache.sis.internal.referencing.GeodeticObjectBuilder) InternationalString(org.opengis.util.InternationalString) SimpleInternationalString(org.apache.sis.util.iso.SimpleInternationalString) CommonCRS(org.apache.sis.referencing.CommonCRS) GeographicCRS(org.opengis.referencing.crs.GeographicCRS)

Aggregations

Length (javax.measure.quantity.Length)44 Angle (javax.measure.quantity.Angle)7 Test (org.junit.Test)7 Color (java.awt.Color)5 Unit (javax.measure.Unit)5 CoordinateSystem (com.revolsys.geometry.cs.CoordinateSystem)4 AffineTransform (java.awt.geom.AffineTransform)4 ArrayList (java.util.ArrayList)4 GeographicCoordinateSystem (com.revolsys.geometry.cs.GeographicCoordinateSystem)3 BoundingBox (com.revolsys.geometry.model.BoundingBox)3 LineString (com.revolsys.geometry.model.LineString)3 Point (com.revolsys.geometry.model.Point)3 Shape (java.awt.Shape)3 Rectangle2D (java.awt.geom.Rectangle2D)3 Path (java.nio.file.Path)3 AtomicLong (java.util.concurrent.atomic.AtomicLong)3 Quantity (javax.measure.Quantity)3 CloseableAffineTransform (com.revolsys.awt.CloseableAffineTransform)2 ProjectedCoordinateSystem (com.revolsys.geometry.cs.ProjectedCoordinateSystem)2 ChainedCoordinatesOperation (com.revolsys.geometry.cs.projection.ChainedCoordinatesOperation)2