Examples with CRS org.apache.sis.referencing.CRS used on opensource projects

Example 1 with CRS

use of org.apache.sis.referencing.CRS in project geotoolkit by Geomatys.

the class DistanceSpliteratorTest method testPolyline.

private void testPolyline(final CoordinateReferenceSystem crs, final boolean isOrthodromic, final boolean parallel) throws IOException, FactoryException {
    final DistanceSpliterator.Builder builder = DistanceSpliterator.builder().setCrs(crs);
    final Travel t = fromSegments(crs);
    builder.setPolyline(t.polyline);
    final Spliterator.OfDouble ds = isOrthodromic ? builder.buildOrthodromic() : builder.buildLoxodromic();
    final int ptNumber = t.polyline.size();
    final AtomicInteger count = new AtomicInteger();
    final double totalDistance = StreamSupport.doubleStream(ds, parallel).peek(whatever -> count.incrementAndGet()).sum();
    Assert.assertEquals("Number of treated segments", ptNumber - 1, count.get());
    org.junit.Assert.assertEquals("Computed distance", isOrthodromic ? t.orthodromicDistance : t.loxodromicDistance, // meter to kilometer
    totalDistance / 1000, isOrthodromic ? EngineTest.EPSI : EngineTest.LOXODROMIC_ERROR);
}

Example 2 with CRS

use of org.apache.sis.referencing.CRS in project geotoolkit by Geomatys.

the class GeometryTransformer method getCoordinates.

/**
 * Check that given geometry is a primitive geometry (implements {@link WithCoordinates}), and
 * get its points.
 * @param source The geometry to extract points from.
 * @return Found points, never null, but can be empty.
 * @throws UnconvertibleObjectException If the given geometry does not implement {@link WithCoordinates}.
 */
private Spliterator<Coordinate> getCoordinates(final Object source) throws UnconvertibleObjectException {
    List<Double> values = null;
    if (source instanceof WithCoordinates) {
        final Coordinates coords = ((WithCoordinates) source).getCoordinates();
        if (coords != null) {
            /* HACK : In GML 3, coordinates are just a list of decimal values.
                 * The grouping by coordinate is done using "srsDimension" attribute
                 * on parent geometry type. However, with GML 2, there's another
                 * possibility : Coordinates use two distinct separators : one for
                 * decimal value separation, and another for coordinate separation.
                 * To manage both ways, we first check if coordinates object has
                 * succeeded in splitting underlying decimals in coordinates. If
                 * it does (return arrays with more than one element), we use
                 * that approach. Otherwise, we fallback on standard way which
                 * will try to determine manually the number of dimensions.
                 */
            Iterator<double[]> it = coords.points().iterator();
            if (it.hasNext() && it.next().length > 1) {
                return coords.points().map(GeometryTransformer::toCoordinate).spliterator();
            } else {
                values = coords.getValues();
            }
        } else if (source instanceof Point) {
            DirectPosition dp = ((Point) source).getPos();
            if (dp != null) {
                return Stream.of(convertDirectPosition(dp)).spliterator();
            } else
                // recognized object, but no value. Empty geometry
                return Spliterators.emptySpliterator();
        }
    }
    // TODO : below conditions should be removed when proper abstraction is setup on GML geometry definition.
    if (values == null) {
        final boolean isLineString = source instanceof org.geotoolkit.gml.xml.LineString;
        final boolean isLineStringSegment = source instanceof org.geotoolkit.gml.xml.LineStringSegment;
        if (isLineString || isLineStringSegment) {
            final DirectPositionList posList;
            if (isLineString) {
                posList = ((org.geotoolkit.gml.xml.LineString) source).getPosList();
            } else {
                posList = ((org.geotoolkit.gml.xml.LineStringSegment) source).getPosList();
            }
            if (posList != null) {
                values = posList.getValue();
            } else {
                final List<? extends DirectPosition> pList;
                if (isLineString) {
                    pList = ((org.geotoolkit.gml.xml.LineString) source).getPos();
                } else {
                    pList = ((org.geotoolkit.gml.xml.LineStringSegment) source).getPos();
                }
                if (pList != null) {
                    return pList.stream().map(GeometryTransformer::convertDirectPosition).filter(Objects::nonNull).spliterator();
                } else {
                    // We've identified a line, but there's no data in it
                    values = Collections.EMPTY_LIST;
                }
            }
        } else if (source instanceof org.geotoolkit.gml.xml.LinearRing) {
            // Note : do not check "getCoordinates", because it should have been done above.
            values = asDoubles(() -> ((org.geotoolkit.gml.xml.LinearRing) source).getPosList());
        } else if (source instanceof org.geotoolkit.gml.xml.v311.GeodesicStringType) {
            values = asDoubles(() -> ((org.geotoolkit.gml.xml.v311.GeodesicStringType) source).getPosList());
        } else if (source instanceof org.geotoolkit.gml.xml.v321.GeodesicStringType) {
            values = asDoubles(() -> ((org.geotoolkit.gml.xml.v321.GeodesicStringType) source).getPosList());
        } else if (source instanceof Curve) {
            CurveSegmentArrayProperty segments = ((Curve) source).getSegments();
            if (segments != null) {
                final List<? extends AbstractCurveSegment> curveSegments = segments.getAbstractCurveSegment();
                if (curveSegments != null) {
                    return curveSegments.stream().flatMap(seg -> StreamSupport.stream(getCoordinates(seg), false)).spliterator();
                }
            }
            // If we arrive here, we've got an empty curve.
            values = Collections.EMPTY_LIST;
        } else if (source instanceof ArcByCenterPointType) {
            final ArcByCenterPointType arc = (ArcByCenterPointType) source;
            org.opengis.geometry.DirectPosition dp = arc.getPos();
            if (dp == null) {
                PointPropertyType pp = arc.getPointProperty();
                if (pp == null) {
                    pp = arc.getPointRep();
                }
                if (pp == null) {
                    throw new UnconvertibleObjectException("Not enough information to build an arc.");
                }
                final Geometry point = new GeometryTransformer(pp.getPoint(), this).get();
                dp = JTS.toEnvelope(point).getLowerCorner();
            }
            CoordinateReferenceSystem crs = dp.getCoordinateReferenceSystem();
            if (crs == null) {
                crs = getSrsName().map(this::findCRS).orElseThrow(() -> new UnconvertibleObjectException("Cannot create an arc without its coordinate reference system"));
                final GeneralDirectPosition gdp = new GeneralDirectPosition(dp);
                gdp.setCoordinateReferenceSystem(crs);
                dp = gdp;
            }
            try {
                final Measure startAngle, endAngle;
                // If we miss a start or end angle, it means we're in presence of a circle.
                if (arc.getStartAngle() == null || arc.getEndAngle() == null) {
                    startAngle = new Measure(0, Units.DEGREE);
                    endAngle = new Measure(360, Units.DEGREE);
                } else {
                    startAngle = asMeasure(arc.getStartAngle());
                    endAngle = asMeasure(arc.getEndAngle());
                }
                final Coordinate[] coordinates = drawArc(dp, asMeasure(arc.getRadius()), startAngle, endAngle, ARC_PRECISION);
                return Spliterators.spliterator(coordinates, Spliterator.ORDERED);
            } catch (TransformException ex) {
                throw new UnconvertibleObjectException("Cannot draw an arc.", ex);
            }
        } else if (source instanceof org.geotoolkit.gml.xml.v311.ArcByCenterPointType) {
            // TODO : factorize with above case
            final org.geotoolkit.gml.xml.v311.ArcByCenterPointType arc = (org.geotoolkit.gml.xml.v311.ArcByCenterPointType) source;
            org.opengis.geometry.DirectPosition dp = arc.getPos();
            if (dp == null) {
                org.geotoolkit.gml.xml.v311.PointPropertyType pp = arc.getPointProperty();
                if (pp == null) {
                    pp = arc.getPointRep();
                }
                if (pp == null) {
                    throw new UnconvertibleObjectException("Not enough information to build an arc.");
                }
                final Geometry point = new GeometryTransformer(pp.getPoint(), this).get();
                dp = JTS.toEnvelope(point).getLowerCorner();
            }
            CoordinateReferenceSystem crs = dp.getCoordinateReferenceSystem();
            if (crs == null) {
                crs = getSrsName().map(this::findCRS).orElseThrow(() -> new UnconvertibleObjectException("Cannot create an arc without its coordinate reference system"));
                final GeneralDirectPosition gdp = new GeneralDirectPosition(dp);
                gdp.setCoordinateReferenceSystem(crs);
                dp = gdp;
            }
            try {
                final Measure startAngle, endAngle;
                // If we miss a start or end angle, it means we're in presence of a circle.
                if (arc.getStartAngle() == null || arc.getEndAngle() == null) {
                    startAngle = new Measure(0, Units.DEGREE);
                    endAngle = new Measure(360, Units.DEGREE);
                } else {
                    startAngle = asMeasure(arc.getStartAngle());
                    endAngle = asMeasure(arc.getEndAngle());
                }
                final Coordinate[] coordinates = drawArc(dp, asMeasure(arc.getRadius()), startAngle, endAngle, ARC_PRECISION);
                return Spliterators.spliterator(coordinates, Spliterator.ORDERED);
            } catch (TransformException ex) {
                throw new UnconvertibleObjectException("Cannot draw an arc.", ex);
            }
        }
    }
    if (values != null) {
        if (values.isEmpty()) {
            return Spliterators.emptySpliterator();
        }
        return new CoordinateSpliterator(values, getCoordinateDimension());
    }
    throw new UnconvertibleObjectException("Cannot extract coordinates from source geometry.");
}

Example 3 with CRS

use of org.apache.sis.referencing.CRS in project geotoolkit by Geomatys.

the class AbstractCoverageSymbolizerRenderer method extractSlice.

private static GridGeometry extractSlice(GridGeometry fullArea, GridGeometry areaOfInterest, final int[] margin, boolean applyResolution) throws DataStoreException, TransformException, FactoryException {
    double[] resolution = areaOfInterest.getResolution(true);
    if (fullArea.isDefined(GridGeometry.RESOLUTION)) {
        CoordinateReferenceSystem crsarea = areaOfInterest.getCoordinateReferenceSystem();
        CoordinateReferenceSystem crsdata = fullArea.getCoordinateReferenceSystem();
        if (CRS.isHorizontalCRS(crsarea) && CRS.isHorizontalCRS(crsdata)) {
            // we are dealing with simple 2D rendering, preserve the canvas grid geometry.
            if (margin != null && Arrays.stream(margin).anyMatch(value -> value != 0)) {
                try {
                    // try to adjust margin
                    // TODO : we should use a GridCoverageResource.subset with a margin value but this isn't implemented yet
                    Envelope env = fullArea.getEnvelope();
                    double[] est = CoverageUtilities.estimateResolution(env, fullArea.getResolution(true), areaOfInterest.getCoordinateReferenceSystem());
                    margin[0] = (int) Math.ceil(margin[0] * (est[0] / resolution[0]));
                    margin[1] = (int) Math.ceil(margin[1] * (est[1] / resolution[1]));
                    areaOfInterest = areaOfInterest.derive().margin(margin).build();
                    // Force rebuilding envelope. Not sure it is really needed however.
                    areaOfInterest = new GridGeometry(areaOfInterest.getExtent(), PixelInCell.CELL_CENTER, areaOfInterest.getGridToCRS(PixelInCell.CELL_CENTER), areaOfInterest.getCoordinateReferenceSystem());
                    return areaOfInterest;
                } catch (Exception e) {
                    LOGGER.log(Level.WARNING, "Cannot compute adapted margin. Artifacts may appear on tile borders");
                    LOGGER.log(Level.FINE, "Details about margin computation failure", e);
                }
            }
        }
    } else {
    // we have no way to apply margin
    // must wait for GridCoverageResource.subset with a margin
    }
    // HACK : This method cannot manage incomplete grid geometries, so we have to skip
    if (!fullArea.isDefined(GridGeometry.ENVELOPE | GridGeometry.GRID_TO_CRS | GridGeometry.EXTENT)) {
        return areaOfInterest;
    }
    // on displayed area
    Envelope canvasEnv = areaOfInterest.getEnvelope();
    if (!applyResolution)
        resolution = null;
    // ///// HACK FOR 0/360 /////////////////////////////////////////////
    try {
        Map.Entry<Envelope, double[]> entry = solveWrapAround(fullArea, canvasEnv, resolution);
        if (entry != null) {
            canvasEnv = entry.getKey();
            resolution = applyResolution ? entry.getValue() : null;
        }
    } catch (ProjectionException ex) {
        // mays happen when displaying an area partialy outside
        // computation area of coverage crs
        LOGGER.log(Level.INFO, ex.getMessage(), ex);
    }
    // ///// HACK FOR 0/360 /////////////////////////////////////////////
    GridGeometry slice = fullArea;
    try {
        GridDerivation derivation = fullArea.derive().rounding(GridRoundingMode.ENCLOSING);
        if (margin != null && margin.length > 0)
            derivation = derivation.margin(margin);
        slice = derivation.subgrid(canvasEnv, resolution).sliceByRatio(1, 0, 1).build();
    } catch (DisjointExtentException ex) {
        throw new DisjointCoverageDomainException(ex.getMessage(), ex);
    } catch (IllegalGridGeometryException ex) {
        throw new DisjointCoverageDomainException(ex.getMessage(), ex);
    }
    return slice;
}

Example 4 with CRS

use of org.apache.sis.referencing.CRS in project geotoolkit by Geomatys.

the class RenderingRoutines method prepareQuery.

/**
 * Creates an optimal query to send to the datastore, knowing which properties are knowned and
 * the appropriate bounding box to filter.
 */
public static Query prepareQuery(final RenderingContext2D renderingContext, FeatureSet fs, final MapLayer layer, final Set<String> styleRequieredAtts, final List<Rule> rules, double symbolsMargin) throws PortrayalException {
    final FeatureType schema;
    try {
        schema = fs.getType();
    } catch (DataStoreException ex) {
        throw new PortrayalException(ex.getMessage(), ex);
    }
    // Note: do not use layer boundary to define the target bbox, because it can be expensive.
    // Anyway, the target resource will be better to determine clipping between rendering boundaries and its own.
    final Envelope bbox = optimizeBBox(renderingContext, fs, symbolsMargin);
    final CoordinateReferenceSystem layerCRS = FeatureExt.getCRS(schema);
    final RenderingHints hints = renderingContext.getRenderingHints();
    /*
         * To restrict queried values to the rendering area, we must identify what geometries are used by the style.
         * For each applied symbol, there are 3 possible cases:
         * - if the rule uses default geometries, they will be added to the geometry property list after the loop
         * - The geometric expression is a value reference, we can safely register it in geometric properties. The
         *   reference xpath is unwrapped in a set to ensure we won't create any doublon filters.
         * - If the geometry property is a complex expression(Ex: a value computed from non geometric fields), we keep
         *   it as is to apply a filter directly upon it. Note that even if it's an expression derived from geometric
         *   fields, we cannot apply spatial filter on them, because the expression could drastically change topology.
         *   For example, if the expression is 'buffer', the result geometry would be larger than any of its operands.
         *   TODO: such cases are maybe manageable by replacing bbox filter by a distance filter based upon the buffer
         *   distance. But would it do more good than harm ?
         */
    boolean isDefaultGeometryNeeded = rules == null || rules.isEmpty();
    final Set<String> geomProperties = new HashSet<>();
    final Set<Expression> complexProperties = new HashSet<>();
    if (rules != null) {
        for (Rule r : rules) {
            for (Symbolizer s : r.symbolizers()) {
                final Expression expGeom = s.getGeometry();
                if (isNil(expGeom))
                    isDefaultGeometryNeeded = true;
                else if (expGeom instanceof ValueReference)
                    geomProperties.add(((ValueReference) expGeom).getXPath());
                else
                    complexProperties.add(expGeom);
            }
        }
    }
    if (isDefaultGeometryNeeded) {
        try {
            final PropertyType defaultGeometry = FeatureExt.getDefaultGeometry(schema);
            final String geomName = Features.getLinkTarget(defaultGeometry).orElseGet(() -> defaultGeometry.getName().toString());
            geomProperties.add(geomName);
        } catch (PropertyNotFoundException e) {
            throw new PortrayalException("Default geometry cannot be determined. " + "However, it is needed to properly define filtering rules.");
        } catch (IllegalStateException e) {
            // If there's multiple geometric properties, and no primary one, we will use them all
            schema.getProperties(true).stream().filter(p -> !Features.getLinkTarget(p).isPresent()).filter(AttributeConvention::isGeometryAttribute).map(p -> p.getName().toString()).forEach(geomProperties::add);
        }
    }
    if (!complexProperties.isEmpty()) {
        LOGGER.fine("A style rule uses complex geometric properties. It can severly affect performance");
    }
    final Optional<Filter> spatialFilter = Stream.concat(geomProperties.stream().map(FILTER_FACTORY::property), complexProperties.stream()).<Filter>map(expression -> FILTER_FACTORY.bbox(expression, bbox)).reduce(FILTER_FACTORY::or);
    Filter userFilter = null;
    // concatenate geographic filter with data filter if there is one
    if (layer != null) {
        Query query = layer.getQuery();
        if (query instanceof FeatureQuery) {
            userFilter = ((FeatureQuery) query).getSelection();
        }
    }
    Filter filter;
    if (spatialFilter.isPresent()) {
        if (userFilter == null)
            filter = spatialFilter.get();
        else
            // Note: we give priority to the spatial filter here, because it is our main use case: rendering is driven
            // by bounding box.
            filter = FILTER_FACTORY.and(spatialFilter.get(), userFilter);
    } else if (userFilter == null) {
        throw new PortrayalException("No spatial filter can be determined from style rules, and no user filter specified." + "We refuse dataset full-scan. To authorize it, manually specify Filter 'INCLUDE' on your map layer.");
    } else {
        LOGGER.warning("Spatial filter cannot be determined for rendering. However, user has provided a custom filter that we'll use as sole filtering policy");
        filter = userFilter;
    }
    final Set<String> copy = new HashSet<>();
    final FeatureType expected;
    final String[] atts;
    if (styleRequieredAtts == null) {
        // all properties are requiered
        expected = schema;
        atts = null;
    } else {
        final Set<String> attributs = styleRequieredAtts;
        copy.addAll(attributs);
        copy.addAll(geomProperties);
        try {
            // always include the identifier if it exist
            schema.getProperty(AttributeConvention.IDENTIFIER);
            copy.add(AttributeConvention.IDENTIFIER);
        } catch (PropertyNotFoundException ex) {
        // no id, ignore it
        }
        atts = copy.toArray(new String[copy.size()]);
        // then we reduce it to the first parent property.
        for (int i = 0; i < atts.length; i++) {
            String attName = atts[i];
            int index = attName.indexOf('/');
            if (index == 0) {
                // remove all xpath elements
                // remove first slash
                attName = attName.substring(1);
                final Pattern pattern = Pattern.compile("(\\{[^\\{\\}]*\\})|(\\[[^\\[\\]]*\\])|/{1}");
                final Matcher matcher = pattern.matcher(attName);
                final StringBuilder sb = new StringBuilder();
                int position = 0;
                while (matcher.find()) {
                    final String match = matcher.group();
                    sb.append(attName.substring(position, matcher.start()));
                    position = matcher.end();
                    if (match.charAt(0) == '/') {
                        // we don't query precisely sub elements
                        position = attName.length();
                        break;
                    } else if (match.charAt(0) == '{') {
                        sb.append(match);
                    } else if (match.charAt(0) == '[') {
                    // strip indexes or xpath searches
                    }
                }
                sb.append(attName.substring(position));
                atts[i] = sb.toString();
            }
        }
        try {
            expected = new ViewMapper(schema, atts).getMappedType();
        } catch (MismatchedFeatureException ex) {
            throw new PortrayalException(ex);
        }
    }
    // combine the filter with rule filters----------------------------------
    if (rules != null) {
        List<Filter<Object>> rulefilters = new ArrayList<>();
        for (Rule rule : rules) {
            if (rule.isElseFilter()) {
                // we can't append styling filters, an else rule match all features
                rulefilters = null;
                break;
            }
            final Filter rf = rule.getFilter();
            if (rf == null || rf == Filter.include()) {
                // we can't append styling filters, this rule matchs all features.
                rulefilters = null;
                break;
            }
            rulefilters.add(rf);
        }
        if (rulefilters != null) {
            final Filter combined;
            if (rulefilters.size() == 1) {
                // we can optimze here, since we pass the filter on the query, we can remove
                // the filter on the rule.
                final MutableRule mr = StyleUtilities.copy(rules.get(0));
                mr.setFilter(null);
                rules.set(0, mr);
                combined = rulefilters.get(0);
            } else {
                combined = FILTER_FACTORY.or(rulefilters);
            }
            if (filter != Filter.include()) {
                filter = FILTER_FACTORY.and(filter, combined);
            } else {
                filter = combined;
            }
        }
    }
    // optimize the filter---------------------------------------------------
    filter = FilterUtilities.prepare(filter, Feature.class, expected);
    final Hints queryHints = new Hints();
    final org.geotoolkit.storage.feature.query.Query qb = new org.geotoolkit.storage.feature.query.Query();
    qb.setTypeName(schema.getName());
    qb.setSelection(filter);
    qb.setProperties(atts);
    // resampling and ignore flag only works when we know the layer crs
    if (layerCRS != null) {
        // add resampling -------------------------------------------------------
        Boolean resample = (hints == null) ? null : (Boolean) hints.get(GO2Hints.KEY_GENERALIZE);
        if (!Boolean.FALSE.equals(resample)) {
            // we only disable resampling if it is explictly specified
            double[] res = renderingContext.getResolution(layerCRS);
            // adjust with the generalization factor
            final Number n = (hints == null) ? null : (Number) hints.get(GO2Hints.KEY_GENERALIZE_FACTOR);
            final double factor;
            if (n != null) {
                factor = n.doubleValue();
            } else {
                factor = GO2Hints.GENERALIZE_FACTOR_DEFAULT.doubleValue();
            }
            res[0] *= factor;
            res[1] *= factor;
            qb.setResolution(res);
            try {
                res = renderingContext.getResolution(CRS.forCode("EPSG:3395"));
                res[0] *= factor;
                res[1] *= factor;
                qb.setLinearResolution(Quantities.create(res[0], Units.METRE));
            } catch (FactoryException ex) {
                throw new PortrayalException(ex.getMessage(), ex);
            }
        }
    // add ignore flag ------------------------------------------------------
    // TODO this is efficient but erases values, when plenty of then are to be rendered
    // we should find another way to handle this
    // if(!GO2Utilities.visibleMargin(rules, 1.01f, renderingContext)){
    // //style does not expend itself further than the feature geometry
    // //that mean geometries smaller than a pixel will not be renderer or barely visible
    // queryHints.put(Hints.KEY_IGNORE_SMALL_FEATURES, renderingContext.getResolution(layerCRS));
    // }
    }
    // add reprojection -----------------------------------------------------
    // we don't reproject, the reprojection may produce curves but JTS can not represent those.
    // so we generate those curves in java2d shapes by doing the transformation ourself.
    // TODO wait for a new geometry implementation
    // qb.setCRS(renderingContext.getObjectiveCRS2D());
    // set the acumulated hints
    qb.setHints(queryHints);
    return qb;
}

Example 5 with CRS

use of org.apache.sis.referencing.CRS in project geotoolkit by Geomatys.

the class SimpleUVSource method atOrigin.

@Override
public Optional<UVSource.TimeSet> atOrigin(DirectPosition origin) {
    final Instant time = getTime(origin).orElseThrow(() -> new IllegalArgumentException("Given position should contain time value"));
    final GridGeometry gg = getSourceGeometry();
    final GeneralEnvelope env = new GeneralEnvelope(gg.getEnvelope());
    boolean pointCrossData = subEnvelope(env, CRS::getHorizontalComponent).filter(subEnv -> {
        try {
            final DirectPosition geoloc = CRS.findOperation(origin.getCoordinateReferenceSystem(), subEnv.getCoordinateReferenceSystem(), null).getMathTransform().transform(origin, null);
            return subEnv.contains(geoloc);
        } catch (FactoryException | TransformException e) {
            throw new BackingStoreException("Cannot project source point in grid coverage horizontal system", e);
        }
    }).isPresent();
    // Short-circuit : no horizontal matching. Stop now
    if (!pointCrossData) {
        LOGGER.log(Level.FINE, "[SimpleUVSource] -> Short-circuit: Given origin {0} does not horizontally intersect envelope {1} of data {2}", new Object[] { origin, env, source });
        return Optional.empty();
    }
    final GeneralEnvelope timeEnv = subEnvelope(env, CRS::getTemporalComponent).orElseThrow(() -> new IllegalStateException("Source dataset has no temporal axis"));
    final DefaultTemporalCRS timeCrs = DefaultTemporalCRS.castOrCopy((TemporalCRS) timeEnv.getCoordinateReferenceSystem());
    timeEnv.getLowerCorner().setOrdinate(0, timeCrs.toValue(time));
    final GridDerivation subgrid = gg.derive().rounding(GridRoundingMode.ENCLOSING).subgrid(env);
    // Note: If input data has an elevation, we'll try to freeze it on a slice : either the one specified at origin, or an arbitrary one.
    final GridGeometry subGrid = freezeElevation(env, origin).map(subgrid::slice).orElse(subgrid).build();
    return Optional.of(new TimeSet(subGrid));
}