/** * Copyright (c) 2013-2020 Contributors to the Eclipse Foundation * * <p> See the NOTICE file distributed with this work for additional information regarding copyright * ownership. All rights reserved. This program and the accompanying materials are made available * under the terms of the Apache License, Version 2.0 which accompanies this distribution and is * available at http://www.apache.org/licenses/LICENSE-2.0.txt */ package org.locationtech.geowave.core.geotime.util; import java.awt.geom.Point2D; import java.lang.reflect.Field; import java.net.MalformedURLException; import java.util.ArrayList; import java.util.Arrays; import java.util.BitSet; import java.util.Collections; import java.util.HashMap; import java.util.HashSet; import java.util.LinkedList; import java.util.List; import java.util.Map; import java.util.Set; import javax.annotation.Nullable; import javax.measure.Unit; import javax.measure.quantity.Length; import org.apache.commons.lang3.tuple.Pair; import org.geotools.factory.CommonFactoryFinder; import org.geotools.factory.GeoTools; import org.geotools.feature.simple.SimpleFeatureBuilder; import org.geotools.geometry.jts.JTS; import org.geotools.geometry.jts.ReferencedEnvelope; import org.geotools.referencing.CRS; import org.geotools.referencing.GeodeticCalculator; import org.geotools.referencing.crs.DefaultGeographicCRS; import org.locationtech.geowave.core.geotime.index.dimension.LatitudeDefinition; import org.locationtech.geowave.core.geotime.index.dimension.LongitudeDefinition; import org.locationtech.geowave.core.geotime.store.dimension.CustomCRSBoundedSpatialDimensionX; import org.locationtech.geowave.core.geotime.store.dimension.CustomCRSBoundedSpatialDimensionY; import org.locationtech.geowave.core.geotime.store.dimension.CustomCRSUnboundedSpatialDimensionX; import org.locationtech.geowave.core.geotime.store.dimension.CustomCRSUnboundedSpatialDimensionY; import org.locationtech.geowave.core.geotime.store.dimension.CustomCrsIndexModel; import org.locationtech.geowave.core.index.GeoWaveSerializationException; import org.locationtech.geowave.core.index.dimension.NumericDimensionDefinition; import org.locationtech.geowave.core.index.sfc.data.BasicNumericDataset; import org.locationtech.geowave.core.index.sfc.data.MultiDimensionalNumericData; import org.locationtech.geowave.core.index.sfc.data.NumericData; import org.locationtech.geowave.core.index.sfc.data.NumericRange; import org.locationtech.geowave.core.index.sfc.data.NumericValue; import org.locationtech.geowave.core.store.api.Index; import org.locationtech.geowave.core.store.data.field.FieldUtils; import org.locationtech.geowave.core.store.query.constraints.BasicQueryByClass.ConstraintData; import org.locationtech.geowave.core.store.query.constraints.BasicQueryByClass.ConstraintSet; import org.locationtech.geowave.core.store.query.constraints.BasicQueryByClass.ConstraintsByClass; import org.locationtech.geowave.core.store.query.constraints.Constraints; import org.locationtech.geowave.core.store.util.ClasspathUtils; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Envelope; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.GeometryCollection; import org.locationtech.jts.geom.GeometryFactory; import org.locationtech.jts.geom.MultiPolygon; import org.locationtech.jts.geom.Point; import org.locationtech.jts.geom.Polygon; import org.locationtech.jts.io.ParseException; import org.locationtech.jts.io.WKBReader; import org.locationtech.jts.io.WKBWriter; import org.opengis.feature.simple.SimpleFeature; import org.opengis.feature.simple.SimpleFeatureType; import org.opengis.filter.FilterFactory2; import org.opengis.filter.spatial.SpatialOperator; import org.opengis.geometry.MismatchedDimensionException; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.opengis.referencing.cs.CoordinateSystem; import org.opengis.referencing.cs.CoordinateSystemAxis; import org.opengis.referencing.operation.MathTransform; import org.opengis.referencing.operation.TransformException; import org.slf4j.Logger; import org.slf4j.LoggerFactory; import com.google.uzaygezen.core.BitSetMath; import si.uom.NonSI; import si.uom.SI; import systems.uom.common.USCustomary; import tec.uom.se.AbstractUnit; import tec.uom.se.unit.AlternateUnit; import tec.uom.se.unit.BaseUnit; import tec.uom.se.unit.Units; /** * This class contains a set of Geometry utility methods that are generally useful throughout the * GeoWave core codebase */ public class GeometryUtils { public static final GeometryFactory GEOMETRY_FACTORY = new GeometryFactory(); private static final Logger LOGGER = LoggerFactory.getLogger(GeometryUtils.class); private static final Object MUTEX = new Object(); private static final Object MUTEX_DEFAULT_CRS = new Object(); private static final int DEFAULT_DIMENSIONALITY = 2; public static final String DEFAULT_CRS_STR = "EPSG:4326"; private static CoordinateReferenceSystem defaultCrsSingleton; private static Set<ClassLoader> initializedClassLoaders = new HashSet<>(); public static final Integer MAX_GEOMETRY_PRECISION = Integer.valueOf(TWKBUtils.MAX_COORD_PRECISION); public static SpatialOperator geometryToSpatialOperator( final Geometry jtsGeom, final String geometryAttributeName) { final FilterFactory2 factory = CommonFactoryFinder.getFilterFactory2(); if (jtsGeom.equalsTopo(jtsGeom.getEnvelope())) { return factory.bbox( factory.property(geometryAttributeName), new ReferencedEnvelope(jtsGeom.getEnvelopeInternal(), GeometryUtils.getDefaultCRS())); } // there apparently is no way to associate a CRS with a poly // intersection operation so it will have to assume the same CRS as the // feature type return factory.intersects(factory.property(geometryAttributeName), factory.literal(jtsGeom)); } @edu.umd.cs.findbugs.annotations.SuppressFBWarnings() public static CoordinateReferenceSystem getDefaultCRS() { if (defaultCrsSingleton == null) { // avoid sync penalty if we can synchronized (MUTEX_DEFAULT_CRS) { // have to do this inside the sync to avoid double init if (defaultCrsSingleton == null) { try { initClassLoader(); defaultCrsSingleton = CRS.decode(DEFAULT_CRS_STR, true); } catch (final Exception e) { LOGGER.error("Unable to decode " + DEFAULT_CRS_STR + " CRS", e); defaultCrsSingleton = DefaultGeographicCRS.WGS84; } } } } return defaultCrsSingleton; } public static void initClassLoader() throws MalformedURLException { synchronized (MUTEX) { final ClassLoader myCl = GeometryUtils.class.getClassLoader(); if (initializedClassLoaders.contains(myCl)) { return; } final ClassLoader classLoader = ClasspathUtils.transformClassLoader(myCl); if (classLoader != null) { GeoTools.addClassLoader(classLoader); } initializedClassLoaders.add(myCl); } } public static ConstraintsByClass basicConstraintsFromGeometry(final Geometry geometry) { final List<ConstraintSet> set = new LinkedList<>(); constructListOfConstraintSetsFromGeometry(geometry, set, false); return new ConstraintsByClass(set); } /** * This utility method will convert a JTS geometry to contraints that can be used in a GeoWave * query. * * @return Constraints as a mapping of NumericData objects representing ranges for a latitude * dimension and a longitude dimension */ public static GeoConstraintsWrapper basicGeoConstraintsWrapperFromGeometry( final Geometry geometry) { final List<ConstraintSet> set = new LinkedList<>(); final boolean geometryConstraintsExactMatch = constructListOfConstraintSetsFromGeometry(geometry, set, true); return new GeoConstraintsWrapper( new ConstraintsByClass(set), geometryConstraintsExactMatch, geometry); } /** * Recursively decompose geometry into a set of envelopes to create a single set. * * @param geometry * @param destinationListOfSets * @param checkTopoEquality */ private static boolean constructListOfConstraintSetsFromGeometry( final Geometry geometry, final List<ConstraintSet> destinationListOfSets, final boolean checkTopoEquality) { // Get the envelope of the geometry being held final int n = geometry.getNumGeometries(); boolean retVal = true; if (n > 1) { retVal = false; for (int gi = 0; gi < n; gi++) { constructListOfConstraintSetsFromGeometry( geometry.getGeometryN(gi), destinationListOfSets, checkTopoEquality); } } else { final Envelope env = geometry.getEnvelopeInternal(); destinationListOfSets.add(basicConstraintSetFromEnvelope(env)); if (checkTopoEquality) { retVal = new GeometryFactory().toGeometry(env).equalsTopo(geometry); } } return retVal; } /** * This utility method will convert a JTS envelope to contraints that can be used in a GeoWave * query. * * @return Constraints as a mapping of NumericData objects representing ranges for a latitude * dimension and a longitude dimension */ public static ConstraintSet basicConstraintSetFromEnvelope(final Envelope env) { // Create a NumericRange object using the x axis final NumericRange rangeLongitude = new NumericRange(env.getMinX(), env.getMaxX()); // Create a NumericRange object using the y axis final NumericRange rangeLatitude = new NumericRange(env.getMinY(), env.getMaxY()); final Map<Class<? extends NumericDimensionDefinition>, ConstraintData> constraintsPerDimension = new HashMap<>(); // Create and return a new IndexRange array with an x and y axis // range final ConstraintData xRange = new ConstraintData(rangeLongitude, false); final ConstraintData yRange = new ConstraintData(rangeLatitude, false); constraintsPerDimension.put(CustomCRSUnboundedSpatialDimensionX.class, xRange); constraintsPerDimension.put(CustomCRSUnboundedSpatialDimensionY.class, yRange); constraintsPerDimension.put(CustomCRSBoundedSpatialDimensionX.class, xRange); constraintsPerDimension.put(CustomCRSBoundedSpatialDimensionY.class, yRange); constraintsPerDimension.put(LongitudeDefinition.class, xRange); constraintsPerDimension.put(LatitudeDefinition.class, yRange); return new ConstraintSet(constraintsPerDimension); } /** * This utility method will convert a JTS envelope to contraints that can be used in a GeoWave * query. * * @return Constraints as a mapping of NumericData objects representing ranges for a latitude * dimension and a longitude dimension */ public static Constraints basicConstraintsFromEnvelope(final Envelope env) { return new ConstraintsByClass(basicConstraintSetFromEnvelope(env)); } /** * This utility method will convert a JTS envelope to that can be used in a GeoWave query. * * @return Constraints as a mapping of NumericData objects representing ranges for a latitude * dimension and a longitude dimension */ public static ConstraintSet basicConstraintsFromPoint( final double latitudeDegrees, final double longitudeDegrees) { // Create a NumericData object using the x axis final NumericData latitude = new NumericValue(latitudeDegrees); // Create a NumericData object using the y axis final NumericData longitude = new NumericValue(longitudeDegrees); final Map<Class<? extends NumericDimensionDefinition>, ConstraintData> constraintsPerDimension = new HashMap<>(); // Create and return a new IndexRange array with an x and y axis // range constraintsPerDimension.put(LongitudeDefinition.class, new ConstraintData(longitude, false)); constraintsPerDimension.put(LatitudeDefinition.class, new ConstraintData(latitude, false)); return new ConstraintSet(constraintsPerDimension); } public static MultiDimensionalNumericData getBoundsFromEnvelope(final Envelope envelope) { final NumericRange[] boundsPerDimension = new NumericRange[2]; boundsPerDimension[0] = new NumericRange(envelope.getMinX(), envelope.getMaxX()); boundsPerDimension[1] = new NumericRange(envelope.getMinY(), envelope.getMaxY()); return new BasicNumericDataset(boundsPerDimension); } /** * Generate a longitude range from a JTS geometry * * @param geometry The JTS geometry * @return The x range */ public static NumericData xRangeFromGeometry(final Geometry geometry) { if ((geometry == null) || geometry.isEmpty()) { return new NumericRange(0, 0); } // Get the envelope of the geometry being held final Envelope env = geometry.getEnvelopeInternal(); // Create a NumericRange object using the x axis return new NumericRange(env.getMinX(), env.getMaxX()); } /** * Generate a latitude range from a JTS geometry * * @param geometry The JTS geometry * @return The y range */ public static NumericData yRangeFromGeometry(final Geometry geometry) { if ((geometry == null) || geometry.isEmpty()) { return new NumericRange(0, 0); } // Get the envelope of the geometry being held final Envelope env = geometry.getEnvelopeInternal(); // Create a NumericRange object using the y axis return new NumericRange(env.getMinY(), env.getMaxY()); } /** * Converts a JTS geometry to binary using JTS a Well Known Binary writer * * @param geometry The JTS geometry * @return The binary representation of the geometry */ public static byte[] geometryToBinary( final Geometry geometry, final @Nullable Integer precision) { if (precision == null) { return new WKBWriter().write(geometry); } return new TWKBWriter(precision).write(geometry); } /** * Converts a byte array as well-known binary to a JTS geometry * * @param binary The well known binary * @return The JTS geometry */ public static Geometry geometryFromBinary( final byte[] binary, final @Nullable Integer precision) { try { if (precision == null) { return new WKBReader().read(binary); } return new TWKBReader().read(binary); } catch (final ParseException e) { throw new GeoWaveSerializationException("Unable to deserialize geometry data", e); } } /** * Converts a byte array as well-known binary to a JTS geometry * * @param binary The well known binary * @return The JTS geometry */ public static Geometry geometryFromBinary( final byte[] binary, final @Nullable Integer precision, final byte serializationVersion) { if (serializationVersion < FieldUtils.SERIALIZATION_VERSION) { try { return new WKBReader().read(binary); } catch (final ParseException e) { LOGGER.warn("Unable to deserialize geometry data", e); throw new GeoWaveSerializationException(e); } } return geometryFromBinary(binary, precision); } /** * This mehtod returns an envelope between negative infinite and positive inifinity in both x and * y * * @return the infinite bounding box */ public static Geometry infinity() { // unless we make this synchronized, we will want to instantiate a new // geometry factory because geometry factories are not thread safe return new GeometryFactory().toGeometry( new Envelope( Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY)); } public static class GeoConstraintsWrapper { private final ConstraintsByClass constraints; private final boolean constraintsMatchGeometry; private final Geometry jtsBounds; public GeoConstraintsWrapper( final ConstraintsByClass constraints, final boolean constraintsMatchGeometry, final Geometry jtsBounds) { this.constraints = constraints; this.constraintsMatchGeometry = constraintsMatchGeometry; this.jtsBounds = jtsBounds; } public ConstraintsByClass getConstraints() { return constraints; } public boolean isConstraintsMatchGeometry() { return constraintsMatchGeometry; } public Geometry getGeometry() { return jtsBounds; } } public static CoordinateReferenceSystem getIndexCrs(final Index[] indices) { CoordinateReferenceSystem indexCrs = null; for (final Index primaryindx : indices) { // for first iteration if (indexCrs == null) { indexCrs = getIndexCrs(primaryindx); } else { if (primaryindx.getIndexModel() instanceof CustomCrsIndexModel) { // check if indexes have different CRS if (!indexCrs.equals(((CustomCrsIndexModel) primaryindx.getIndexModel()).getCrs())) { LOGGER.error("Multiple indices with different CRS is not supported"); throw new RuntimeException("Multiple indices with different CRS is not supported"); } else { if (!indexCrs.equals(getDefaultCRS())) { LOGGER.error("Multiple indices with different CRS is not supported"); throw new RuntimeException("Multiple indices with different CRS is not supported"); } } } } } return indexCrs; } public static CoordinateReferenceSystem getIndexCrs(final Index index) { CoordinateReferenceSystem indexCrs = null; if (index.getIndexModel() instanceof CustomCrsIndexModel) { indexCrs = ((CustomCrsIndexModel) index.getIndexModel()).getCrs(); } else { indexCrs = getDefaultCRS(); } return indexCrs; } public static String getCrsCode(final CoordinateReferenceSystem crs) { return (CRS.toSRS(crs)); } /** * Build a buffer around a geometry * * @param crs * @param geometry * @param distanceUnits * @param distance * @return the buffered geometry and the degrees that it was buffered * @throws TransformException */ public static final Pair<Geometry, Double> buffer( final CoordinateReferenceSystem crs, final Geometry geometry, final String distanceUnits, final double distance) throws TransformException { Unit<Length> unit; try { unit = lookup(distanceUnits); } catch (final Exception e) { unit = Units.METRE; LOGGER.warn("Cannot lookup unit of measure " + distanceUnits, e); } final double meterDistance = unit.getConverterTo(Units.METRE).convert(distance); final double degrees = distanceToDegrees(crs, geometry, meterDistance); // buffer does not respect the CRS; it uses simple cartesian math. // nor does buffer handle dateline boundaries return Pair.of(adjustGeo(crs, geometry.buffer(degrees)), degrees); } public static Unit<Length> lookup(final String name) { final String lowerCaseName = name.toLowerCase(); Unit<Length> unit = lookup(SI.class, lowerCaseName); if (unit != null) { return unit; } unit = lookup(NonSI.class, lowerCaseName); if (unit != null) { return unit; } if (lowerCaseName.endsWith("s")) { return lookup(lowerCaseName.substring(0, lowerCaseName.length() - 1)); } if (lowerCaseName.startsWith("kilo") && (lowerCaseName.length() > 4)) { final Unit<Length> u = lookup(lowerCaseName.substring(4)); if (u != null) { return u.multiply(1000); } } // if we get here, try some aliases if (lowerCaseName.equals("feet")) { return USCustomary.FOOT; } // if we get here, try some aliases if (lowerCaseName.equals("meter")) { return Units.METRE; } if (lowerCaseName.equals("unity")) { return (Unit) AbstractUnit.ONE; } return null; } private static Unit<Length> lookup(final Class class1, final String name) { Unit<Length> unit = null; final Field[] fields = class1.getDeclaredFields(); for (int i = 0; i < fields.length; i++) { final Field field = fields[i]; final String name2 = field.getName(); if ((field.getType().isAssignableFrom(BaseUnit.class) || field.getType().isAssignableFrom(AlternateUnit.class)) && name2.equalsIgnoreCase(name)) { try { unit = (Unit<Length>) field.get(unit); return unit; } catch (final Exception e) { } } } return unit; } /** * Consume a geometry that may be over the ranges of the CRS (e.g date-line crossing). Adjust for * crossings with a multi-polygon instance where each contained polygon represents a portion of * the provided geometry longitude value. Clip hemisphere crossings (fix TBD). * * @param crs * @param geometry * @return the adjusted geometry */ public static Geometry adjustGeo(final CoordinateReferenceSystem crs, final Geometry geometry) { final List<Polygon> polygons = fixRangeOfCoordinates(crs, geometry); if (polygons.size() == 1) { return polygons.get(0); } return geometry.getFactory().createMultiPolygon(polygons.toArray(new Polygon[polygons.size()])); } /** * Adjust geometry so that coordinates fit into long/lat bounds. * * <p> Split date-line crossing polygons. * * <p> For now, clip hemisphere crossing portions of the polygon. * * @param geometry * @return list valid polygons */ public static List<Polygon> fixRangeOfCoordinates( final CoordinateReferenceSystem crs, final Geometry geometry) { final List<Polygon> replacements = new ArrayList<>(); if (geometry instanceof MultiPolygon) { final MultiPolygon multi = (MultiPolygon) geometry; for (int i = 0; i < multi.getNumGeometries(); i++) { final Geometry geo = multi.getGeometryN(i); replacements.addAll(fixRangeOfCoordinates(crs, geo)); } return replacements; } // collection is more general than multi-polygon else if (geometry instanceof GeometryCollection) { final GeometryCollection multi = (GeometryCollection) geometry; for (int i = 0; i < multi.getNumGeometries(); i++) { final Geometry geo = multi.getGeometryN(i); replacements.addAll(fixRangeOfCoordinates(crs, geo)); } return replacements; } final Coordinate[] geoCoords = geometry.getCoordinates(); final Coordinate modifier = findModifier(crs, geoCoords); replacements.addAll(constructGeometriesOverMapRegions(modifier, geometry)); return replacements; } /** * update modifier for each axis of the coordinate where the modifier's axis is less extreme than * the provides coordinate * * @param modifier * @param cood */ private static void updateModifier(final Coordinate coord, final Coordinate modifier) { for (int i = 0; i < 3; i++) { double coordOrdinateValue, modifierOrdinateValue; switch (i) { case 1: coordOrdinateValue = coord.getY(); modifierOrdinateValue = modifier.getY(); break; case 2: coordOrdinateValue = coord.getZ(); modifierOrdinateValue = modifier.getZ(); break; default: case 0: coordOrdinateValue = coord.getX(); modifierOrdinateValue = modifier.getX(); break; } if (!Double.isNaN(coordOrdinateValue) && !Double.isNaN(modifierOrdinateValue)) { if (Math.abs(modifierOrdinateValue) < Math.abs(coordOrdinateValue)) { modifier.setOrdinate(i, coord.getOrdinate(i)); } } } } /** * Build a modifier that, when added to the coordinates of a polygon, moves invalid sections of * the polygon to a valid portion of the map. * * @param crs * @param coords * @return */ private static Coordinate findModifier( final CoordinateReferenceSystem crs, final Coordinate[] coords) { final Coordinate maxModifier = new Coordinate(0, 0, 0); for (final Coordinate coord : coords) { final Coordinate modifier = diff(adjustCoordinateToFitInRange(crs, coord), coord); updateModifier(modifier, maxModifier); } return maxModifier; } /** * Produce a set of polygons for each region of the map corrected for date line and hemisphere * crossings. Due to the complexity of going around the hemisphere, clip the range. * * <p> Consider a polygon that cross both the hemisphere in the north and the date line in the * west (-182 92, -182 88, -178 88, -178 92, -182 92). The result is two polygons: (-180 90, -180 * 88, -178 88, -178 90, -180 90) (180 90, 180 88, 178 88, 178 90, 180 90) * * @param modifier * @param geometry - a geometry that may cross date line and/or hemispheres. * @return the set of polygons */ public static List<Polygon> constructGeometriesOverMapRegions( final Coordinate modifier, final Geometry geometry) { final Coordinate[] geoCoords = geometry.getCoordinates(); final List<Polygon> polygons = new LinkedList<>(); final Geometry world = world(geometry.getFactory(), GeometryUtils.getDefaultCRS()); // First do the polygon unchanged world final Geometry worldIntersections = world.intersection(geometry); for (int i = 0; i < worldIntersections.getNumGeometries(); i++) { final Polygon polyToAdd = (Polygon) worldIntersections.getGeometryN(i); if (!polygons.contains(polyToAdd)) { polygons.add(polyToAdd); } } // now use the modifier...but just the x axis for longitude // optimization...do not modify if 0 if (Math.abs(modifier.x) > 0.0000000001) { final Coordinate[] newCoords = new Coordinate[geoCoords.length]; int c = 0; for (final Coordinate geoCoord : geoCoords) { newCoords[c++] = new Coordinate(geoCoord.x + modifier.x, geoCoord.y, geoCoord.z); } final Polygon transposedPoly = geometry.getFactory().createPolygon(newCoords); final Geometry adjustedPolyWorldIntersections = world.intersection(transposedPoly); for (int i = 0; i < adjustedPolyWorldIntersections.getNumGeometries(); i++) { final Polygon polyToAdd = (Polygon) adjustedPolyWorldIntersections.getGeometryN(i); if (!polygons.contains(polyToAdd)) { polygons.add(polyToAdd); } } } return polygons; } /** * Make sure the coordinate falls in the range of provided coordinate reference systems's * coordinate system. 'x' coordinate is wrapped around date line. 'y' and 'z' coordinate are * clipped. At some point, this function will be adjusted to project 'y' appropriately. * * @param crs * @param coord * @return the adjusted coordinate */ public static Coordinate adjustCoordinateToFitInRange( final CoordinateReferenceSystem crs, final Coordinate coord) { return new Coordinate( adjustCoordinateDimensionToRange(coord.getX(), crs, 0), clipRange(coord.getY(), crs, 1), clipRange(coord.getZ(), crs, 2)); } /** * @param coord1 * @param coord2 subtracted from coord1 * @return a coordinate the supplies the difference of values for each axis between coord1 and * coord2 */ private static Coordinate diff(final Coordinate coord1, final Coordinate coord2) { return new Coordinate( coord1.getX() - coord2.getX(), coord1.getY() - coord2.getY(), coord1.getZ() - coord2.getZ()); } /** * @param val the value * @param crs * @param axis the coordinate axis * @return */ private static double clipRange( final double val, final CoordinateReferenceSystem crs, final int axis) { final CoordinateSystem coordinateSystem = crs.getCoordinateSystem(); if (coordinateSystem.getDimension() > axis) { final CoordinateSystemAxis coordinateAxis = coordinateSystem.getAxis(axis); if (val < coordinateAxis.getMinimumValue()) { return coordinateAxis.getMinimumValue(); } else if (val > coordinateAxis.getMaximumValue()) { return coordinateAxis.getMaximumValue(); } } return val; } /** * This is perhaps a brain dead approach to do this, but it does handle wrap around cases. Also * supports cases where the wrap around occurs many times. * * @param val the value * @param crs * @param axis the coordinate axis * @return the adjusted coordinate dimension */ public static double adjustCoordinateDimensionToRange( final double val, final CoordinateReferenceSystem crs, final int axis) { final CoordinateSystem coordinateSystem = crs.getCoordinateSystem(); if (coordinateSystem.getDimension() > axis) { final double lowerBound = coordinateSystem.getAxis(axis).getMinimumValue(); final double bound = coordinateSystem.getAxis(axis).getMaximumValue() - lowerBound; final double sign = sign(val); // re-scale to 0 to n, then determine how many times to 'loop // around' final double mult = Math.floor(Math.abs((val + (sign * (-1.0 * lowerBound))) / bound)); return val + (mult * bound * sign * (-1.0)); } return val; } private static double sign(final double val) { return val < 0 ? -1 : 1; } /** * Return a multi-polygon representing the bounded map regions split by the axis * * @param factory * @param crs * @return a world geometry */ public static Geometry world(final GeometryFactory factory, final CoordinateReferenceSystem crs) { return factory.createPolygon(toPolygonCoordinates(crs.getCoordinateSystem())); } private static Coordinate[] toPolygonCoordinates(final CoordinateSystem coordinateSystem) { final Coordinate[] coordinates = new Coordinate[(int) Math.pow(2, coordinateSystem.getDimension()) + 1]; final BitSet greyCode = new BitSet(coordinateSystem.getDimension()); final BitSet mask = getGreyCodeMask(coordinateSystem.getDimension()); for (int i = 0; i < coordinates.length; i++) { coordinates[i] = new Coordinate( getValue(greyCode, coordinateSystem.getAxis(0), 0), getValue(greyCode, coordinateSystem.getAxis(1), 1), coordinateSystem.getDimension() > 2 ? getValue(greyCode, coordinateSystem.getAxis(2), 2) : Double.NaN); grayCode(greyCode, mask); } return coordinates; } private static BitSet getGreyCodeMask(final int dims) { final BitSet mask = new BitSet(dims); for (int i = 0; i < dims; i++) { mask.set(i); } return mask; } private static void grayCode(final BitSet code, final BitSet mask) { BitSetMath.grayCodeInverse(code); BitSetMath.increment(code); code.and(mask); BitSetMath.grayCode(code); } private static double getValue( final BitSet set, final CoordinateSystemAxis axis, final int dimension) { return (set.get(dimension)) ? axis.getMaximumValue() : axis.getMinimumValue(); } /** * Convert meters to decimal degrees based on widest point * * @throws TransformException */ private static double distanceToDegrees( final CoordinateReferenceSystem crs, final Geometry geometry, final double meters) throws TransformException { final GeometryFactory factory = geometry.getFactory(); return (geometry instanceof Point) ? geometry.distance(farthestPoint(crs, (Point) geometry, meters)) : distanceToDegrees( crs, geometry.getEnvelopeInternal(), factory == null ? new GeometryFactory() : factory, meters); } private static double distanceToDegrees( final CoordinateReferenceSystem crs, final Envelope env, final GeometryFactory factory, final double meters) throws TransformException { return Collections.max( Arrays.asList( distanceToDegrees( crs, factory.createPoint(new Coordinate(env.getMaxX(), env.getMaxY())), meters), distanceToDegrees( crs, factory.createPoint(new Coordinate(env.getMaxX(), env.getMinY())), meters), distanceToDegrees( crs, factory.createPoint(new Coordinate(env.getMinX(), env.getMinY())), meters), distanceToDegrees( crs, factory.createPoint(new Coordinate(env.getMinX(), env.getMaxY())), meters))); } /** farther point in longitudinal axis given a latitude */ private static Point farthestPoint( final CoordinateReferenceSystem crs, final Point point, final double meters) { final GeodeticCalculator calc = new GeodeticCalculator(crs); calc.setStartingGeographicPoint(point.getX(), point.getY()); calc.setDirection(90, meters); Point2D dest2D = calc.getDestinationGeographicPoint(); // if this flips over the date line then try the other direction if (dest2D.getX() < point.getX()) { calc.setDirection(-90, meters); dest2D = calc.getDestinationGeographicPoint(); } return point.getFactory().createPoint(new Coordinate(dest2D.getX(), dest2D.getY())); } public static SimpleFeature crsTransform( final SimpleFeature entry, final SimpleFeatureType reprojectedType, final MathTransform transform) { SimpleFeature crsEntry = entry; if (transform != null) { // we can use the transform we have already calculated for this // feature try { // this will clone the feature and retype it to Index CRS crsEntry = SimpleFeatureBuilder.retype(entry, reprojectedType); // this will transform the geometry crsEntry.setDefaultGeometry( JTS.transform((Geometry) entry.getDefaultGeometry(), transform)); } catch (MismatchedDimensionException | TransformException e) { LOGGER.warn( "Unable to perform transform to specified CRS of the index, the feature geometry will remain in its original CRS", e); } } return crsEntry; } }