blob: 119af5219e74a6409e5ed61e755f828cfe4917ec [file] [log] [blame]
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You 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.
*/
package org.apache.sis.metadata.iso.extent;
import java.util.Date;
import java.util.List;
import java.util.ArrayList;
import javax.measure.Unit;
import org.opengis.geometry.Envelope;
import org.opengis.geometry.DirectPosition;
import org.opengis.temporal.TemporalPrimitive;
import org.opengis.metadata.extent.Extent;
import org.opengis.metadata.extent.VerticalExtent;
import org.opengis.metadata.extent.TemporalExtent;
import org.opengis.metadata.extent.BoundingPolygon;
import org.opengis.metadata.extent.GeographicExtent;
import org.opengis.metadata.extent.GeographicBoundingBox;
import org.opengis.metadata.extent.GeographicDescription;
import org.opengis.referencing.cs.CoordinateSystemAxis;
import org.opengis.referencing.cs.AxisDirection;
import org.opengis.referencing.crs.VerticalCRS;
import org.opengis.referencing.crs.GeographicCRS;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.datum.VerticalDatum;
import org.opengis.referencing.datum.VerticalDatumType;
import org.opengis.referencing.operation.TransformException;
import org.apache.sis.internal.metadata.ReferencingServices;
import org.apache.sis.metadata.InvalidMetadataException;
import org.apache.sis.measure.Longitude;
import org.apache.sis.measure.MeasurementRange;
import org.apache.sis.measure.Range;
import org.apache.sis.util.resources.Vocabulary;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.ComparisonMode;
import org.apache.sis.util.Static;
import static java.lang.Math.*;
import static org.apache.sis.internal.metadata.ReferencingServices.AUTHALIC_RADIUS;
/**
* Convenience static methods for extracting information from {@link Extent} objects.
* This class provides methods for:
*
* <ul>
* <li>{@linkplain #getGeographicBoundingBox Fetching geographic},
* {@linkplain #getVerticalRange vertical} or
* {@linkplain #getDate temporal components} in a convenient form.</li>
* <li>Computing {@linkplain #intersection intersection} of bounding boxes</li>
* <li>Computing {@linkplain #area area} estimations.</li>
* </ul>
*
* @author Martin Desruisseaux (Geomatys)
* @version 0.8
*
* @see org.apache.sis.geometry.Envelopes
*
* @since 0.3
* @module
*/
public final class Extents extends Static {
/**
* Do no allow instantiation of this class.
*/
private Extents() {
}
/**
* A geographic extent ranging from 180°W to 180°E and 90°S to 90°N.
* This extent has no vertical and no temporal components.
*/
public static final Extent WORLD;
static {
final DefaultGeographicBoundingBox box = new DefaultGeographicBoundingBox(-180, 180, -90, 90);
box.transitionTo(DefaultGeographicBoundingBox.State.FINAL);
final DefaultExtent world = new DefaultExtent(
Vocabulary.formatInternational(Vocabulary.Keys.World), box, null, null);
world.transitionTo(DefaultExtent.State.FINAL);
WORLD = world;
}
/**
* Returns a single geographic bounding box from the specified extent.
* This method tries to find the bounding box in the cheapest way
* before to fallback on more expansive computations:
*
* <ol>
* <li>First, this method searches geographic elements that are instance of {@link GeographicBoundingBox}.<ul>
* <li>If exactly one such instance is found, then this method returns that instance directly (no copy).</li>
* <li>If more than one instance is found, then this method computes and returns the
* {@linkplain DefaultGeographicBoundingBox#add union} of all bounding boxes.</li>
* </ul></li>
* <li>If above step found no {@code GeographicBoundingBox}, then this method inspects geographic elements
* that are instance of {@link BoundingPolygon}, taking in account only the envelopes associated to a
* coordinate reference system of kind {@link GeographicCRS}. If such envelopes are found, then this
* method computes and returns their union.</li>
* <li>If above step found no polygon's envelope associated to a geographic CRS, then in last resort this
* method uses all polygon's envelopes regardless their coordinate reference system (provided that the
* CRS is not null), applying coordinate transformations if needed.</li>
* <li>If above step found no polygon's envelope, then this method returns {@code null}.</li>
* </ol>
*
* @param extent the extent to convert to a geographic bounding box, or {@code null}.
* @return a geographic bounding box extracted from the given extent, or {@code null} if none.
*
* @see org.apache.sis.referencing.CRS#getDomainOfValidity(CoordinateReferenceSystem)
*/
public static GeographicBoundingBox getGeographicBoundingBox(final Extent extent) {
GeographicBoundingBox bounds = null;
if (extent != null) {
DefaultGeographicBoundingBox modifiable = null;
final List<Envelope> fallbacks = new ArrayList<>();
for (final GeographicExtent element : extent.getGeographicElements()) {
/*
* If a geographic bounding box can be obtained, add it to the previous boxes (if any).
* All exclusion boxes before the first inclusion box are ignored.
*/
if (element instanceof GeographicBoundingBox) {
final GeographicBoundingBox item = (GeographicBoundingBox) element;
if (bounds == null) {
/*
* We use DefaultGeographicBoundingBox.getInclusion(Boolean) below because
* add(…) method that we use cares about the case where inclusion is false.
*/
if (DefaultGeographicBoundingBox.getInclusion(item.getInclusion())) {
bounds = item;
}
} else {
if (modifiable == null) {
bounds = modifiable = new DefaultGeographicBoundingBox(bounds);
}
modifiable.add(item);
}
}
}
/*
* If we found not explicit GeographicBoundingBox element, use the information that we
* collected in BoundingPolygon elements. This may involve coordinate transformations.
*/
if (bounds == null) try {
for (final Envelope envelope : fallbacks) {
final DefaultGeographicBoundingBox item = new DefaultGeographicBoundingBox();
item.setBounds(envelope);
if (bounds == null) {
bounds = item;
} else {
if (modifiable == null) {
bounds = modifiable = new DefaultGeographicBoundingBox(bounds);
}
modifiable.add(item);
}
}
} catch (TransformException e) {
throw new InvalidMetadataException(Errors.format(Errors.Keys.CanNotTransformEnvelope), e);
}
}
return bounds;
}
/**
* Returns the union of chosen vertical ranges found in the given extent, or {@code null} if none.
* This method gives preference to heights above the Mean Sea Level when possible.
* Depths have negative height values: if the
* {@linkplain org.apache.sis.referencing.cs.DefaultCoordinateSystemAxis#getDirection() axis direction}
* is toward down, then this method reverses the sign of minimum and maximum values.
*
* <h4>Multi-occurrences</h4>
* If the given {@code Extent} object contains more than one vertical extent, then this method
* performs a choice based on the vertical datum and the unit of measurement:
*
* <ul class="verbose">
* <li><p><b>Choice based on vertical datum</b><br>
* Only the extents associated (indirectly, through their CRS) to the same non-null {@link VerticalDatumType}
* will be taken in account. If all datum types are null, then this method conservatively uses only the first
* vertical extent. Otherwise the datum type used for filtering the vertical extents is:</p>
*
* <ul>
* <li>{@link VerticalDatumType#GEOIDAL} or {@link VerticalDatumType#DEPTH DEPTH} if at least one extent
* uses those datum types. For this method, {@code DEPTH} is considered as equivalent to {@code GEOIDAL}
* except for the axis direction.</li>
* <li>Otherwise, the first non-null datum type found in iteration order.</li>
* </ul>
*
* <div class="note"><b>Rational:</b> like {@linkplain #getGeographicBoundingBox(Extent) geographic bounding box},
* the vertical range is an approximated information; the range returned by this method does not carry any
* information about the vertical CRS and this method does not attempt to perform coordinate transformation.
* But this method is more useful if the returned ranges are close to a frequently used surface, like the
* Mean Sea Level. The same simplification is applied in the
* <a href="http://docs.opengeospatial.org/is/12-063r5/12-063r5.html#31">{@code VerticalExtent} element of
* Well Known Text (WKT) format</a>, which specifies that <cite>“Vertical extent is an approximate description
* of location; heights are relative to an unspecified mean sea level.”</cite></div></li>
*
* <li><p><b>Choice based on units of measurement</b><br>
* If, after the choice based on the vertical datum described above, there is still more than one vertical
* extent to consider, then the next criterion checks for the units of measurement.</p>
* <ul>
* <li>If no range specify a unit of measurement, return the first range and ignore all others.</li>
* <li>Otherwise take the first range having a unit of measurement. Then:<ul>
* <li>All other ranges having an incompatible unit of measurement will be ignored.</li>
* <li>All other ranges having a compatible unit of measurement will be converted to
* the unit of the first retained range, and their union will be computed.</li>
* </ul></li>
* </ul>
*
* <div class="note"><b>Example:</b>
* Heights or depths are often measured using some pressure units, for example hectopascals (hPa).
* An {@code Extent} could contain two vertical elements: one with the height measurements in hPa,
* and the other element with heights transformed to metres using an empirical formula.
* In such case this method will select the first vertical element on the assumption that it is
* the "main" one that the metadata producer intended to show. Next, this method will search for
* other vertical elements using pressure unit. In our example there is none, but if such elements
* were found, this method would compute their union.</div></li>
* </ul>
*
* @param extent the extent to convert to a vertical measurement range, or {@code null}.
* @return a vertical measurement range created from the given extent, or {@code null} if none.
*
* @since 0.4
*/
public static MeasurementRange<Double> getVerticalRange(final Extent extent) {
MeasurementRange<Double> range = null;
VerticalDatumType selectedType = null;
if (extent != null) {
for (final VerticalExtent element : extent.getVerticalElements()) {
double min = element.getMinimumValue();
double max = element.getMaximumValue();
final VerticalCRS crs = element.getVerticalCRS();
VerticalDatumType type = null;
Unit<?> unit = null;
if (crs != null) {
final VerticalDatum datum = crs.getDatum();
if (datum != null) {
type = datum.getVerticalDatumType();
if (VerticalDatumType.DEPTH.equals(type)) {
type = VerticalDatumType.GEOIDAL;
}
}
final CoordinateSystemAxis axis = crs.getCoordinateSystem().getAxis(0);
unit = axis.getUnit();
if (AxisDirection.DOWN.equals(axis.getDirection())) {
final double tmp = min;
min = -max;
max = -tmp;
}
}
if (range != null) {
/*
* If the new range does not specify any datum type or unit, then we do not know how to
* convert the values before to perform the union operation. Conservatively do nothing.
*/
if (type == null || unit == null) {
continue;
}
/*
* If the new range is not measured relative to the same kind of surface than the previous range,
* then we do not know how to combine those ranges. Do nothing, unless the new range is a Mean Sea
* Level Height in which case we forget all previous ranges and use the new one instead.
*/
if (!type.equals(selectedType)) {
if (!type.equals(VerticalDatumType.GEOIDAL)) {
continue;
}
} else if (selectedType != null) {
/*
* If previous range did not specify any unit, then unconditionally replace it by
* the new range since it provides more information. If both ranges specify units,
* then we will compute the union if we can, or ignore the new range otherwise.
*/
final Unit<?> previous = range.unit();
if (previous != null) {
if (previous.isCompatible(unit)) {
range = (MeasurementRange<Double>) range.union(
MeasurementRange.create(min, true, max, true, unit));
}
continue;
}
}
}
range = MeasurementRange.create(min, true, max, true, unit);
selectedType = type;
}
}
return range;
}
/**
* Returns the union of all time ranges found in the given extent, or {@code null} if none.
*
* @param extent the extent to convert to a time range, or {@code null}.
* @return a time range created from the given extent, or {@code null} if none.
*
* @since 0.4
*/
public static Range<Date> getTimeRange(final Extent extent) {
Date min = null;
Date max = null;
if (extent != null) {
for (final TemporalExtent t : extent.getTemporalElements()) {
final Date startTime, endTime;
if (t instanceof DefaultTemporalExtent) {
final DefaultTemporalExtent dt = (DefaultTemporalExtent) t;
startTime = dt.getStartTime(); // Maybe user has overridden those methods.
endTime = dt.getEndTime();
} else {
final TemporalPrimitive p = t.getExtent();
startTime = DefaultTemporalExtent.getTime(p, true);
endTime = DefaultTemporalExtent.getTime(p, false);
}
if (startTime != null && (min == null || startTime.before(min))) min = startTime;
if ( endTime != null && (max == null || endTime.after (max))) max = endTime;
}
}
if (min == null && max == null) {
return null;
}
return new Range<>(Date.class, min, true, max, true);
}
/**
* Returns an instant in the {@linkplain Extent#getTemporalElements() temporal elements} of the given extent,
* or {@code null} if none. First, this method computes the union of all temporal elements. Then this method
* computes the linear interpolation between the start and end time as in the following pseudo-code:
*
* {@preformat java
* return new Date(startTime + (endTime - startTime) * location);
* }
*
* Special cases:
* <ul>
* <li>If {@code location} is 0, then this method returns the {@linkplain DefaultTemporalExtent#getStartTime() start time}.</li>
* <li>If {@code location} is 1, then this method returns the {@linkplain DefaultTemporalExtent#getEndTime() end time}.</li>
* <li>If {@code location} is 0.5, then this method returns the average of start time and end time.</li>
* <li>If {@code location} is outside the [0 … 1] range, then the result will be outside the temporal extent.</li>
* </ul>
*
* @param extent the extent from which to get an instant, or {@code null}.
* @param location 0 for the start time, 1 for the end time, 0.5 for the average time, or the
* coefficient (usually in the [0 … 1] range) for interpolating an instant.
* @return an instant interpolated at the given location, or {@code null} if none.
*
* @since 0.4
*/
public static Date getDate(final Extent extent, final double location) {
ArgumentChecks.ensureFinite("location", location);
Date min = null;
Date max = null;
if (extent != null) {
for (final TemporalExtent t : extent.getTemporalElements()) {
Date startTime = null;
Date endTime = null;
if (t instanceof DefaultTemporalExtent) {
final DefaultTemporalExtent dt = (DefaultTemporalExtent) t;
if (location != 1) startTime = dt.getStartTime(); // Maybe user has overridden those methods.
if (location != 0) endTime = dt.getEndTime();
} else {
final TemporalPrimitive p = t.getExtent();
if (location != 1) startTime = DefaultTemporalExtent.getTime(p, true);
if (location != 0) endTime = DefaultTemporalExtent.getTime(p, false);
}
if (startTime != null && (min == null || startTime.before(min))) min = startTime;
if ( endTime != null && (max == null || endTime.after (max))) max = endTime;
}
}
if (min == null) return max;
if (max == null) return min;
final long startTime = min.getTime();
return new Date(Math.addExact(startTime, Math.round((max.getTime() - startTime) * location)));
}
/**
* Returns the position at the median longitude and latitude values of the given bounding box.
* This method does not check the {@linkplain DefaultGeographicBoundingBox#getInclusion() inclusion} status.
* This method takes in account bounding boxes that cross the anti-meridian.
*
* @param bbox the bounding box for which to get the median longitude and latitude values, or {@code null}.
* @return a median position of the given bounding box, or {@code null} if none.
*/
public static DirectPosition centroid(final GeographicBoundingBox bbox) {
if (bbox != null) {
double y = (bbox.getNorthBoundLatitude() + bbox.getSouthBoundLatitude()) / 2;
double x = bbox.getWestBoundLongitude();
double xmax = bbox.getEastBoundLongitude();
if (xmax < x) {
xmax += (Longitude.MAX_VALUE - Longitude.MIN_VALUE);
}
x = Longitude.normalize((x + xmax) / 2);
if (Double.isFinite(x) || Double.isFinite(y)) {
return ReferencingServices.getInstance().geographic(x, y);
}
}
return null;
}
/**
* Returns an <em>estimation</em> of the area (in square metres) of the given bounding box.
* Since {@code GeographicBoundingBox} provides only approximated information (for example
* it does not specify the datum), the value returned by this method is also approximated.
*
* <p>The current implementation performs its computation on the
* {@linkplain org.apache.sis.referencing.CommonCRS#SPHERE GRS 1980 Authalic Sphere}.
* However this may change in any future SIS version.</p>
*
* @param box the geographic bounding box for which to compute the area, or {@code null}.
* @return an estimation of the area in the given bounding box (m²),
* or {@linkplain Double#NaN NaN} if the given box was null.
*
* @since 0.4
*/
public static double area(final GeographicBoundingBox box) {
if (box == null) {
return Double.NaN;
}
double Δλ = box.getEastBoundLongitude() - box.getWestBoundLongitude();
final double span = Longitude.MAX_VALUE - Longitude.MIN_VALUE;
if (Δλ > span) {
Δλ = span;
} else if (Δλ < 0) {
if (Δλ < -span) {
Δλ = -span;
} else {
Δλ += span;
}
}
return (AUTHALIC_RADIUS * AUTHALIC_RADIUS) * toRadians(Δλ) *
max(0, sin(toRadians(box.getNorthBoundLatitude())) -
sin(toRadians(box.getSouthBoundLatitude())));
}
/**
* Returns the intersection of the given geographic bounding boxes. If any of the arguments is {@code null},
* then this method returns the other argument (which may be null). Otherwise this method returns a box which
* is the intersection of the two given boxes.
*
* <p>This method never modify the given boxes, but may return directly one of the given arguments if it
* already represents the intersection result.</p>
*
* @param b1 the first bounding box, or {@code null}.
* @param b2 the second bounding box, or {@code null}.
* @return the intersection (may be any of the {@code b1} or {@code b2} argument if unchanged),
* or {@code null} if the two given boxes are null.
* @throws IllegalArgumentException if the {@linkplain DefaultGeographicBoundingBox#getInclusion() inclusion status}
* is not the same for both boxes.
*
* @see DefaultGeographicBoundingBox#intersect(GeographicBoundingBox)
*
* @since 0.4
*/
public static GeographicBoundingBox intersection(final GeographicBoundingBox b1, final GeographicBoundingBox b2) {
if (b1 == null) return b2;
if (b2 == null || b2 == b1) return b1;
final DefaultGeographicBoundingBox box = new DefaultGeographicBoundingBox(b1);
box.intersect(b2);
if (box.equals(b1, ComparisonMode.BY_CONTRACT)) return b1;
if (box.equals(b2, ComparisonMode.BY_CONTRACT)) return b2;
return box;
}
/**
* May compute an intersection between the given geographic extents.
* Current implementation supports only {@link GeographicBoundingBox};
* all other kinds are handled as if they were {@code null}.
*
* <p>We may improve this method in future Apache SIS version, but it is not yet clear how.
* For example how to handle {@link GeographicDescription} or {@link BoundingPolygon}?
* This method should not be public before we find a better contract.</p>
*/
static GeographicExtent intersection(final GeographicExtent e1, final GeographicExtent e2) {
return intersection(e1 instanceof GeographicBoundingBox ? (GeographicBoundingBox) e1 : null,
e2 instanceof GeographicBoundingBox ? (GeographicBoundingBox) e2 : null);
}
/**
* Returns the intersection of the given vertical extents. If any of the arguments is {@code null},
* then this method returns the other argument (which may be null). Otherwise this method returns a
* vertical extent which is the intersection of the two given extents.
*
* <p>This method never modify the given extents, but may return directly one of the given arguments
* if it already represents the intersection result.</p>
*
* <h4>Advantage and inconvenient of this method</h4>
* This method can not intersect extents defined with different datums because height transformations
* generally require the geodetic positions (latitudes and longitudes) of the heights to transform.
* For more general transformations, it is better to convert all extent components into a single envelope,
* then {@linkplain org.apache.sis.geometry.Envelopes#transform(CoordinateOperation, Envelope) transform
* the envelope at once}. On the other hand, this {@code intersect(…)} method preserves better
* the {@link org.apache.sis.xml.NilReason} (if any).
*
* @param e1 the first extent, or {@code null}.
* @param e2 the second extent, or {@code null}.
* @return the intersection (may be any of the {@code e1} or {@code e2} argument if unchanged),
* or {@code null} if the two given extents are null.
* @throws IllegalArgumentException if the two extents do not use the same datum, ignoring metadata.
*
* @see DefaultVerticalExtent#intersect(VerticalExtent)
*
* @since 0.8
*/
public static VerticalExtent intersection(final VerticalExtent e1, final VerticalExtent e2) {
if (e1 == null) return e2;
if (e2 == null || e2 == e1) return e1;
final DefaultVerticalExtent extent = new DefaultVerticalExtent(e1);
extent.intersect(e2);
if (extent.equals(e1, ComparisonMode.BY_CONTRACT)) return e1;
if (extent.equals(e2, ComparisonMode.BY_CONTRACT)) return e2;
return extent;
}
/**
* Returns the intersection of the given temporal extents. If any of the arguments is {@code null},
* then this method returns the other argument (which may be null). Otherwise this method returns a
* temporal extent which is the intersection of the two given extents.
*
* <p>This method never modify the given extents, but may return directly one of the given arguments
* if it already represents the intersection result.</p>
*
* @param e1 the first extent, or {@code null}.
* @param e2 the second extent, or {@code null}.
* @return the intersection (may be any of the {@code e1} or {@code e2} argument if unchanged),
* or {@code null} if the two given extents are null.
* @throws UnsupportedOperationException if no implementation of {@code TemporalFactory} has been found
* on the classpath.
*
* @see DefaultTemporalExtent#intersect(TemporalExtent)
*
* @since 0.8
*/
public static TemporalExtent intersection(final TemporalExtent e1, final TemporalExtent e2) {
if (e1 == null) return e2;
if (e2 == null || e2 == e1) return e1;
final DefaultTemporalExtent extent = new DefaultTemporalExtent(e1);
extent.intersect(e2);
if (extent.equals(e1, ComparisonMode.BY_CONTRACT)) return e1;
if (extent.equals(e2, ComparisonMode.BY_CONTRACT)) return e2;
return extent;
}
/**
* Returns the intersection of the given extents. If any of the arguments is {@code null},
* then this method returns the other argument (which may be null). Otherwise this method
* returns an extent which is the intersection of all geographic, vertical and temporal
* elements in the two given extents.
*
* <p>This method never modify the given extents, but may return directly one of the given
* arguments if it already represents the intersection result.</p>
*
* @param e1 the first extent, or {@code null}.
* @param e2 the second extent, or {@code null}.
* @return the intersection (may be any of the {@code e1} or {@code e2} argument if unchanged),
* or {@code null} if the two given extents are null.
* @throws IllegalArgumentException if two elements to intersect are not compatible (e.g. mismatched
* {@linkplain DefaultGeographicBoundingBox#getInclusion() bounding box inclusion status} or
* mismatched {@linkplain DefaultVerticalExtent#getVerticalCRS() vertical datum}).
* @throws UnsupportedOperationException if a {@code TemporalFactory} is required but no implementation
* has been found on the classpath.
*
* @see DefaultExtent#intersect(Extent)
*
* @since 0.8
*/
public static Extent intersection(final Extent e1, final Extent e2) {
if (e1 == null) return e2;
if (e2 == null || e2 == e1) return e1;
final DefaultExtent extent = new DefaultExtent(e1);
extent.intersect(e2);
if (extent.equals(e1, ComparisonMode.BY_CONTRACT)) return e1;
if (extent.equals(e2, ComparisonMode.BY_CONTRACT)) return e2;
return extent;
}
}