org.opentripplanner.graph_builder.module.ned.VerticalDatum 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.io.BufferedInputStream;
import java.io.DataInputStream;
import java.io.IOException;
import java.io.InputStream;

/**
 * A Vertical datum specified as a grid of offsets from NAD83 http://vdatum.noaa.gov/dev/gtx_info.html
 *
 * @author novalis
 */
public class VerticalDatum {

  double lowerLeftLatitude;

  double lowerLeftLongitude;

  double deltaLatitude;

  double deltaLongitude;

  float[][] datum;

  public VerticalDatum(
    double lowerLeftLongitude,
    double lowerLeftLatitude,
    double width,
    double height,
    float[][] datum
  ) {
    this.lowerLeftLongitude = lowerLeftLongitude;
    this.lowerLeftLatitude = lowerLeftLatitude;
    this.deltaLongitude = width;
    this.deltaLatitude = height;
    this.datum = datum;
  }

  public static VerticalDatum fromGTX(InputStream inputStream) throws IOException {
    DataInputStream stream = new DataInputStream(new BufferedInputStream(inputStream));
    double lowerLeftLatitude = stream.readDouble();
    double lowerLeftLongitude = stream.readDouble();
    if (lowerLeftLongitude > 180) {
      lowerLeftLongitude -= 360; //convert to standard coordinates
    }
    double deltaLatitude = stream.readDouble();
    double deltaLongitude = stream.readDouble();
    int nRows = stream.readInt();
    int nColumns = stream.readInt();
    float[][] data = new float[nRows][nColumns];
    for (int y = 0; y < nRows; ++y) {
      for (int x = 0; x < nColumns; ++x) {
        data[y][x] = stream.readFloat();
      }
    }
    return new VerticalDatum(
      lowerLeftLongitude,
      lowerLeftLatitude,
      deltaLatitude * nColumns,
      deltaLongitude * nRows,
      data
    );
  }

  double interpolatedHeight(double longitude, double latitude) {
    // because VDatums can cross -180 longitude, there's some complication in interpolating
    // longitude

    double lowerLeftAdjusted = lowerLeftLongitude;
    if (lowerLeftLongitude + deltaLongitude > 180) {
      // then lowerLeftLongitude must be left of 180
      if (longitude < lowerLeftLongitude) {
        // lowerLeftLongitude is right of 180 (probablY)
        lowerLeftAdjusted -= 360; // this transforms lowerLeft so that interpolation will
        // work
        if (longitude < lowerLeftLongitude) {
          throw new RuntimeException("longitude out of range");
        }
      }
    }

    if (longitude > lowerLeftAdjusted + deltaLongitude || longitude < lowerLeftAdjusted) {
      throw new RuntimeException("longitude out of range");
    }
    if (latitude < lowerLeftLatitude || latitude > lowerLeftLatitude + deltaLatitude) {
      throw new RuntimeException("latitude out of range");
    }

    double longitudeNormalized = (longitude - lowerLeftAdjusted) / deltaLongitude;
    double latitudeNormalized = (latitude - lowerLeftLatitude) / deltaLatitude;

    int rows = datum.length;
    int columns = datum[0].length;

    int x1 = (int) Math.floor(longitudeNormalized * columns);
    int y1 = (int) Math.floor(latitudeNormalized * rows);

    double gridXFraction = longitudeNormalized * columns - x1;
    double gridYFraction = latitudeNormalized * rows - y1;
    return (
      datum[y1][x1] *
      gridXFraction *
      gridYFraction +
      datum[y1][x1 + 1] *
      (1 - gridXFraction) *
      gridYFraction +
      datum[y1 + 1][x1] *
      gridXFraction *
      (1 - gridYFraction) +
      datum[y1 + 1][x1 + 1] *
      (1 - gridXFraction) *
      (1 - gridYFraction)
    );
  }

  boolean covers(double longitude, double latitude) {
    double lowerLeftAdjusted = lowerLeftLongitude;
    if (lowerLeftLongitude + deltaLongitude > 180) {
      // then lowerLeftLongitude must be left of 180
      if (longitude < lowerLeftLongitude) {
        // lowerLeftLongitude is right of 180 (probablY)
        lowerLeftAdjusted -= 360; // this transforms lowerLeft so that interpolation will
        // work
        if (longitude < lowerLeftLongitude) {
          return false;
        }
      }
    }

    if (longitude > lowerLeftAdjusted + deltaLongitude || longitude < lowerLeftAdjusted) {
      return false;
    }
    if (latitude < lowerLeftLatitude || latitude > lowerLeftLatitude + deltaLatitude) {
      return false;
    }
    return true;
  }
}