view client/src/lib/geo.js @ 2264:627bfcbbf631

client: Draw fairway dimensions for all LoS and in more "cut cases".
author Raimund Renkert <raimund.renkert@intevation.de>
date Thu, 14 Feb 2019 13:53:54 +0100
parents ca33ad696594
children 940ae7c20326
line wrap: on
line source

/* This is Free Software under GNU Affero General Public License v >= 3.0
 * without warranty, see README.md and license for details.
 *
 * SPDX-License-Identifier: AGPL-3.0-or-later
 * License-Filename: LICENSES/AGPL-3.0.txt
 *
 * Copyright (C) 2018 by via donau
 *   – Österreichische Wasserstraßen-Gesellschaft mbH
 * Software engineering by Intevation GmbH
 *
 * Author(s):
 * Thomas Junk <thomas.junk@intevation.de>
 */

/**
 *
 * Distance calculations
 * JS transposition of cross.go functions
 *
 */

import { GeoJSON } from "ol/format.js";
import Feature from "ol/Feature";
import distance from "@turf/distance";
import {
  lineString as turfLineString,
  polygon as turfPolygon
} from "@turf/helpers";
import lineSplit from "@turf/line-split";
import lineSlice from "@turf/line-slice";
import lineIntersect from "@turf/line-intersect";
import booleanWithin from "@turf/boolean-within";
import booleanContains from "@turf/boolean-contains";

const EARTHRADIUS = 6378137.0;

/**
 * Converts to radiant
 * @param {degree} d
 */
const deg2rad = d => {
  return (d * Math.PI) / 180.0;
};

/**
 * Calculates the difference between two points in m
 *
 * Points are given with {lat: $lat, lon: $lon}
 *
 * @param {object} P1
 * @param {object} P2
 */
const distanceBetween = (P1, P2) => {
  const dLat = deg2rad(P2.lat - P1.lat);
  let dLng = Math.abs(deg2rad(P2.lon - P1.lon));
  if (dLng > Math.PI) {
    dLng = 2 * Math.PI - dLng;
  }
  const x = dLng * Math.cos(deg2rad((P1.lat + P2.lat) / 2.0));
  return Math.sqrt(dLat * dLat + x * x) * EARTHRADIUS;
};

/**
 * Takes a triple of [lat, long, alt] and generates
 * an object with according attributes
 * {
 *  lat: $lat,
 *  lon: $lon,
 *  alt: $alt
 * }
 *
 * @param {array} coords
 */
const Point = coords => {
  return {
    lon: coords[0],
    lat: coords[1],
    alt: coords[2]
  };
};

/**
 * Has geoJSON as its input and transforms
 * given coordinates into points representing
 * distance from startpoint / altitude information
 *
 * a) extracting the minimum altitude
 * b) extracting the maximum altitude
 * c) calculating the total length of the given profile
 * d) transposes the datapoints relative to a given start point
 *
 * The calculation of total equals the sum of partial distances between points
 *
 * The x-value of a point is equal to the total distance up to this point
 *
 * The distance between the last point of the last segment and the end point is added
 * to the total
 *
 * @param {object} geoJSON, startPoint, endPoint
 */
const transform = ({ geoJSON, startPoint, endPoint }) => {
  const lineSegments = geoJSON.geometry.coordinates;
  let segmentPoints = [];
  let lengthPolyLine = 0;
  let referencePoint = Point(startPoint);
  let minAlt = Math.abs(lineSegments[0][0][2]);
  let maxAlt = Math.abs(lineSegments[0][0][2]);
  let currentPoint = null;
  for (let segment of lineSegments) {
    let points = [];
    for (let coordinateTriplet of segment) {
      currentPoint = Point(coordinateTriplet);
      lengthPolyLine += distanceBetween(referencePoint, currentPoint);
      let y = Math.abs(currentPoint.alt);
      points.push({
        x: lengthPolyLine,
        y: y
      });
      if (y < minAlt) minAlt = y;
      if (y > maxAlt) maxAlt = y;
      referencePoint = currentPoint;
    }
    segmentPoints.push(points);
  }
  lengthPolyLine += distanceBetween(currentPoint, Point(endPoint));
  return { segmentPoints, lengthPolyLine, minAlt, maxAlt };
};

/**
 * Prepare profile takes geoJSON data in form of
 * a MultiLineString, e.g.
 *
 * {
 *   type: "Feature",
 *   geometry: {
 *   type: "MultiLineString",
 *   coordinates: [
 *                 [
 *                  [16.53593398, 48.14694085, -146.52392755]
 *                  ...
 *                  ]]
 *
 * and transforms it to a structure representing the number of sections
 * where data is present with according lengths and the points
 *
 * {
 *  { points:
 *           [
 *            [ { x: 0.005798201616417183, y: -146.52419461 },
 *              { x: 0, y: -146.52394016 }
 *              ...
 *            ]
 *           ]
 *   lengthPolyLine: 160.06814078495722,
 *   minAlt: -146.73122231,
 *   maxAlt: -145.65155866
 *  }
 *
 * @param {object} geoJSON
 */
const prepareProfile = ({ geoJSON, startPoint, endPoint }) => {
  const { segmentPoints, lengthPolyLine, minAlt, maxAlt } = transform({
    geoJSON,
    startPoint,
    endPoint
  });
  return {
    points: segmentPoints,
    lengthPolyLine: lengthPolyLine,
    minAlt: minAlt,
    maxAlt: maxAlt
  };
};

const generateFeatureRequest = (profileLine, bottleneck_id, date_info) => {
  const feature = new Feature({
    geometry: profileLine,
    bottleneck: bottleneck_id,
    date: date_info
  });
  return new GeoJSON({ geometryName: "geometry" }).writeFeature(feature);
};

const calculateFairwayCoordinates = (profileLine, fairwayGeometry, depth) => {
  // both geometries have to be in EPSG:4326
  // uses turfjs distance() function
  let fairwayCoordinates = [];
  var line = turfLineString(profileLine.getCoordinates());
  var polygon = turfPolygon(fairwayGeometry.getCoordinates());
  var intersects = lineIntersect(line, polygon);

  var startPoint, endPoint

  if (intersects.features.length == 2) {
    let line1 = lineSlice(intersects.features[0].geometry, intersects.features[1].geometry, line);
    startPoint = line1.geometry.coordinates[0];
    endPoint = line1.geometry.coordinates[1];
  } else if (booleanWithin(line, polygon)) {
    startPoint = line.geometry.coordinates[0];
    endPoint = line.geometry.coordinates[1];
  } else {
    let split = lineSplit(line, polygon);
    for (let feature of split.features) {
      if (booleanWithin(feature, polygon)){
        startPoint = feature.geometry.coordinates[0];
        endPoint = feature.geometry.coordinates[1];
      }
      if (booleanContains(polygon, feature)){
        startPoint = feature.geometry.coordinates[0];
        endPoint = feature.geometry.coordinates[1];
      }
    }
  }

  let opts = { units: "kilometers" };
  let pStartPoint = profileLine.getCoordinates()[0];
  fairwayCoordinates.push([
    distance(pStartPoint, startPoint, opts) * 1000,
    distance(pStartPoint, endPoint, opts) * 1000,
    depth
  ]);

  return fairwayCoordinates;
};

const featureToFairwayCoordinates = (feature, profileLine) => {
  // need to use EPSG:3857 which is the proj of vectorSource
  var intersectingPolygon = feature
    .getGeometry()
    .clone()
    .transform("EPSG:3857", "EPSG:4326");
  const fairwayCoordinates = calculateFairwayCoordinates(
    profileLine,
    intersectingPolygon,
    feature.get("min_depth") ? feature.get("min_depth") / 100 : 2.5
  );
  return fairwayCoordinates;
};

export {
  generateFeatureRequest,
  prepareProfile,
  calculateFairwayCoordinates,
  featureToFairwayCoordinates
};