org.opentripplanner.graph_builder.module.ned.UnifiedGridCoverage Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of otp Show documentation
The OpenTripPlanner multimodal journey planning system
There is a newer version: 2.6.0
package org.opentripplanner.graph_builder.module.ned;

import java.util.ArrayList;
import java.util.List;
import javax.media.jai.InterpolationBilinear;
import org.geotools.coverage.AbstractCoverage;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.Interpolator2D;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.index.SpatialIndex;
import org.locationtech.jts.index.strtree.STRtree;
import org.opengis.coverage.CannotEvaluateException;
import org.opengis.coverage.Coverage;
import org.opengis.coverage.PointOutsideCoverageException;
import org.opengis.coverage.SampleDimension;
import org.opengis.geometry.DirectPosition;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

/**
 * Stitches together multiple elevation maps into a single elevation map, hackily.  This is
 * horrible, but the geotools way of doing things is too slow.
 *
 * @author novalis
 */
public class UnifiedGridCoverage extends AbstractCoverage {

  private static final long serialVersionUID = -7798801307087575896L;

  private static final Logger log = LoggerFactory.getLogger(UnifiedGridCoverage.class);

  /**
   * A spatial index of the intersection of all regions and datums. For smaller-scale deployments,
   * this spatial index might perform more slowly than manually iterating over each region and
   * datum. However, in larger regions, this can result in 20+% more calculations being done during
   * the same time compared to the brute-force method. In smaller-scale deployments since the
   * overall time is not as long, we leave this in here for the benefit of larger regions where this
   * will result in much better performance.
   */
  private final SpatialIndex datumRegionIndex;
  private final ArrayList regions;

  /**
   * It would be nice if we could construct this unified coverage with zero sub-coverages and add
   * all sub-coverages in the same way. However, the superclass constructor (AbstractCoverage) needs
   * a coverage to copy properties from. So the first sub-coverage needs to be passed in at
   * construction time.
   */
  protected UnifiedGridCoverage(List regionCoverages, List datums) {
    super("unified", regionCoverages.get(0));
    regions = new ArrayList<>();
    datumRegionIndex = new STRtree();

    // Iterate through region coverages, creating Interpolators for each region and then adding them and the
    // intersected datum to the spatial index.
    for (GridCoverage2D regionCoverage : regionCoverages) {
      // TODO might bicubic interpolation give better results?
      GridCoverage2D regionCoverageInterpolator = Interpolator2D.create(
        regionCoverage,
        new InterpolationBilinear()
      );

      // Iterate over datums to find intersection envelope with each region and add to spatial index.
      for (VerticalDatum datum : datums) {
        Envelope datumEnvelope = new Envelope(
          datum.lowerLeftLongitude,
          datum.lowerLeftLongitude + datum.deltaLongitude,
          datum.lowerLeftLatitude,
          datum.lowerLeftLatitude + datum.deltaLatitude
        );
        ReferencedEnvelope regionEnvelope = new ReferencedEnvelope(
          regionCoverageInterpolator.getEnvelope()
        );
        Envelope intersection = regionEnvelope.intersection(datumEnvelope);
        datumRegionIndex.insert(intersection, new DatumRegion(datum, regionCoverageInterpolator));
      }
      regions.add(regionCoverageInterpolator);
    }
  }

  @Override
  public Object evaluate(DirectPosition point) throws CannotEvaluateException {
    /* we don't use this function, we use evaluate(DirectPosition point, double[] values) */
    return null;
  }

  @Override
  public int getNumSampleDimensions() {
    return regions.get(0).getNumSampleDimensions();
  }

  @Override
  public SampleDimension getSampleDimension(int index) throws IndexOutOfBoundsException {
    return regions.get(0).getSampleDimension(index);
  }

  /**
   * Calculate the elevation at a given point
   */
  public double[] evaluate(DirectPosition point, double[] values) throws CannotEvaluateException {
    double x = point.getOrdinate(0);
    double y = point.getOrdinate(1);
    Coordinate pointCoordinate = new Coordinate(x, y);
    Envelope envelope = new Envelope(pointCoordinate);
    List coverageCandidates = datumRegionIndex.query(envelope);
    if (coverageCandidates.size() > 0) {
      // Found a match for coverage/datum.
      DatumRegion datumRegion = coverageCandidates.get(0);
      double[] result = datumRegion.region.evaluate(point, values);
      result[0] += datumRegion.datum.interpolatedHeight(x, y);
      return result;
    }
    throw new PointOutsideCoverageException("Point not found: " + point);
  }

  public static class DatumRegion {

    public final VerticalDatum datum;
    public final Coverage region;

    public DatumRegion(VerticalDatum datum, Coverage region) {
      this.datum = datum;
      this.region = region;
    }
  }
}