/*
 * RiskScape™ Copyright New Zealand Institute for Earth Science Limited
 * (Earth Sciences New Zealand) is distributed for research purposes only
 * under the terms of AGPLv3.
 *
 * RiskScape™ Copyright 2026 New Zealand Institute for Earth Science
 * Limited (Earth Sciences New Zealand). All rights reserved. Source code
 * available under the AGPLv3.
 * 
 * This program is free software: you can redistribute it and/or modify it under
 *  the terms of the GNU Affero General Public License as published by the Free
 *  Software Foundation, either version 3 of the License, or (at your option) any
 *  later version.
 * 
 * This program is distributed for RESEARCH PURPOSES ONLY, in the hope that it will
 * be useful for research and education initiatives.
 * 
 * If you are not a researcher, or you are a researcher who wishes to use this
 * program on terms other than AGPLv3 (including those who wish to restrict the
 * distribution of any source code created using this program), please contact:
 * https://riskscape.org.nz
 * 
 * This program is distributed WITHOUT ANY WARRANTY; without even the implied
 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Affero General Public License for more details.  You should have received a copy
 * of the GNU Affero General Public License along with this program.  If not, see
 * <http://www.gnu.org/licenses/>.
 * 
 * By way of summary only, under the AGPLv3:
 *     • Permissions of this strongest copyleft license are conditioned
 *       on making available complete source code of licensed works and
 *       modifications, which include larger works using a licensed work,
 *       under the same license.
 *     • Copyright and license notices must be preserved.
 *     • Contributors provide an express grant of patent rights.
 *     • When a modified version is used to provide a service over a
 *       network, the complete source code of the modified version must be made
 *       available.
 */
package nz.org.riskscape.engine.geo;

import static org.geotools.measure.Units.*;

import java.util.ArrayList;
import java.util.List;
import java.util.Optional;
import java.util.function.Consumer;

import javax.measure.Unit;
import javax.measure.UnitConverter;
import javax.measure.quantity.Length;

import org.geotools.geometry.jts.GeometryCollector;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.referencing.CRS;
import org.geotools.measure.Units;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryCollection;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Lineal;
import org.locationtech.jts.geom.Polygonal;
import org.locationtech.jts.geom.Puntal;
import org.locationtech.jts.operation.overlayng.OverlayNG;
import org.geotools.api.referencing.FactoryException;
import org.geotools.api.referencing.crs.CoordinateReferenceSystem;
import org.geotools.api.referencing.operation.TransformException;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.LineSegment;
import org.locationtech.jts.geom.LineString;
import org.locationtech.jts.geom.Point;

import nz.org.riskscape.engine.SRIDSet;
import nz.org.riskscape.engine.types.Nullable;
import nz.org.riskscape.engine.types.Referenced;
import nz.org.riskscape.engine.types.Type;

public class GeometryUtils {

  public static final double METRES_PER_NAUTICAL_MILE = NAUTICAL_MILE.getConverterTo(METRE).convert(1D);

  // A nautical mile is defined as 1/60 of a degree on any longitudinal line, or the equator.
  // So we can calculate how many metres there are in a degree.
  public static final double METRES_PER_DEGREE = METRES_PER_NAUTICAL_MILE * 60;

  /**
   * Applies the processor {@link Consumer} to the {@link Geometry} or to each part if geometry is
   * a {@link GeometryCollection}.
   *
   * @param geom geometry to process
   * @param processor consumer to process geometries
   */
  public static void processPerPart(Geometry geom, Consumer<Geometry> processor) {
    if (geom instanceof GeometryCollection) {
      GeometryCollection collection = (GeometryCollection) geom;
        for (int i = 0; i < collection.getNumGeometries(); i++) {
          Geometry part = collection.getGeometryN(i);
          processor.accept(part);
        }
    } else {
      processor.accept(geom);
    }
  }

  /**
   * Converts linealMetres to the units of the target {@link CoordinateReferenceSystem}.
   *
   * @param linealMetres length in metres
   * @param crs {@link CoordinateReferenceSystem} to convert to
   * @return length in CRS units
   */
  public static double toCrsUnits(double linealMetres, CoordinateReferenceSystem crs) {
    double lineal = linealMetres;
    Unit<?> crsUnit = crs.getCoordinateSystem().getAxis(0).getUnit();

    if (crsUnit.equals(Units.DEGREE_ANGLE)) {
      // If CRS is angular, we need to convert with out METRES_PER_DEGREE constant
      lineal = linealMetres / METRES_PER_DEGREE;
    } else if (!crsUnit.equals(METRE)) {
      // CRS is some non metric unit. We need to convert metres to that unit.
      Unit<Length> lengthUnit = crsUnit.asType(Length.class);
      UnitConverter metresToCrsUnit = METRE.getConverterTo(lengthUnit);
      lineal = metresToCrsUnit.convert(linealMetres);
    }
    return lineal;
  }

  /**
   * Geometry types are considered equal by this method
   * if they are both {@link Puntal}, both {@link Lineal} or both {@link Polygonal}.  This allows comparison of geometry
   * type without having to consider any {@link GeometryCollection} types, i.e. MultiPoint and Point are equivalent.
   */
  public static boolean equivalentTypes(Geometry g1, Geometry g2) {
    return GeometryFamily.fromClass(g1.getClass()) == GeometryFamily.fromClass(g2.getClass());
  }

  /**
   * Converts a distance that is expected to be in metres, to the units of the geom's
   * {@link CoordinateReferenceSystem}.
   *
   * Note that if the geometry's {@link CoordinateReferenceSystem} is not projected (a regular grid) then the
   * returned distance (in CRS units) is an approximation. Because angular units don't result in the same
   * linear distance at all parts of the world.
   *
   * @param distance to convert
   * @param geom geometry to get the CRS from
   * @param sridSet to obtain CRS information from
   * @return distance converted to the correct units for the geom's CRS
   */
  public static double distanceToCrsUnits(double distance, Geometry geom, SRIDSet sridSet) {
    CoordinateReferenceSystem crs = sridSet.get(geom.getSRID());
    return toCrsUnits(distance, crs);
  }

  /**
   * Tests if re-projecting features from source to target {@link Referenced} types is safe and this
   * would be required . To re-project safely the source and target must be in different projections and
   * either target needs to have a global bounds (be EPSG:4326) or the source has bounds and they
   * fit within the the target bounds.
   *
   * @param sourceGeomType
   * @param targetGeomType
   * @return true if safe to re-project, false otherwise or if both types have the same projection
   */
  public static boolean canReprojectSafely(Type sourceGeomType, Type targetGeomType) {
    return Nullable.strip(sourceGeomType).find(Referenced.class)
        .map(sourceReferenced -> {
          return Nullable.strip(targetGeomType).find(Referenced.class).map(toRef -> {
            if (CRS.equalsIgnoreMetadata(sourceReferenced.getCrs(), toRef.getCrs())) {
              // reprojection is not needed if they are already in the same projection.
              return false;
            }
            try {
              if (CRS.lookupIdentifier(toRef.getCrs(), false).equals("EPSG:4326")) {
                // any CRS should be able to be re-projected to WSG84
                return true;
              }
            } catch (FactoryException e) {}

            ReferencedEnvelope sourceBounds = sourceReferenced.getBounds();
            if (sourceBounds == null) {
              return false;
            }
            try {
              ReferencedEnvelope transformed = sourceBounds.transform(toRef.getCrs(), true);
              ReferencedEnvelope targetBounds = toBounds(toRef.getCrs());
              if (targetBounds != null) {
                return targetBounds.contains((org.locationtech.jts.geom.Envelope) transformed);
              }
              return false;
            } catch (TransformException | FactoryException t) {
              // ReferencedEvelope#transform will throw exceptions if parts of the envelope cannot be transformed.
              // That is some points are outside of projected bounds.
              return false;
            }
          }).orElse(false);
        })
        .orElse(false);
  }

  /**
   * Attempt to get the {@link ReferencedEnvelope) that forms the bounds of the projected area.
   *
   * @param crs
   * @return ReferencedEnvelope that forms the bounds of the projected area.
   */
  private static ReferencedEnvelope toBounds(CoordinateReferenceSystem crs) {
    org.geotools.api.geometry.Bounds targetEnv = CRS.getEnvelope(crs);
    if (targetEnv != null) {
      return new ReferencedEnvelope(
          targetEnv.getMinimum(0), targetEnv.getMaximum(0),
          targetEnv.getMinimum(1), targetEnv.getMaximum(1),
          crs);
    }
    // hmm we didn't find it. This could be because crs has force-XY set.
    try {
      // try to find the projection with force-XY turned off.
      CoordinateReferenceSystem noForce = CRS.decode(CRS.lookupIdentifier(crs, false), false);
      targetEnv = CRS.getEnvelope(noForce);
      if (targetEnv != null) {
        // Cool we found the bounds.
        ReferencedEnvelope transformedEnv = new ReferencedEnvelope(
            targetEnv.getMinimum(0), targetEnv.getMaximum(0),
            targetEnv.getMinimum(1), targetEnv.getMaximum(1),
            noForce);
        return transformedEnv.transform(crs, true);
      }
    } catch (FactoryException | TransformException ex) {}
    return null;
  }

  /**
   * Return a {@link Geometry} that contains all the members of geom that are of the expectedGeomFamily. This
   * method should only be called when the expected geometry family is known with certainty. Such as after
   * the difference of two polygons which should result in a polygonal result.
   *
   * If no members of the expected family exist then an empty geometry of the expected family will be returned.
   *
   * This is useful for processing the result of geometry overlay operations (e.g difference)
   * which may not guarantee to return a homogenus geometry. Refer to {@link OverlayNG} for a discussion on
   * strict mode.
   *
   * @param geom {@link Geometry} to remove non-family members from
   * @param expectedGeomFamily {@link GeometryFamily} to allow in the processed result
   * @return {@link Geometry} that contains all the members of geom that are of the expectedGeomFamily.
   */
  public static Geometry removeNonFamilyMembers(Geometry geom, GeometryFamily expectedGeomFamily) {
    if (!(geom instanceof GeometryCollection)) {
      if (expectedGeomFamily.isSameFamily(geom)) {
        // If the result is not a GeometryCollection then it must be homegeneous, lets just return it.
        return geom;
      } else {
        // but not if it isn't in the right family
        return emptyGeom(expectedGeomFamily, geom.getFactory());
      }
    }

    List<Geometry> parts = new ArrayList<>();
    GeometryUtils.processPerPart(geom, part -> {
      if (expectedGeomFamily.isSameFamily(part)) {
        // we only want to keep parts that are from the same geometry family as expected
        parts.add(part);
      }
    });
    if (parts.isEmpty()) {
      // If there are no parts remaining we need to return an empty geometry (point, line, polygon)
      // Otherwise we'd return an anonymous geometry collection with no members and that would be bad.
      return emptyGeom(expectedGeomFamily, geom.getFactory());
    } else if (parts.size() == 1) {
      // If there is only one part remaining then there is no need to wrap it up in a collection.
      return parts.get(0);
    }

    // Now that we know that the parts are homogeneous we can leave it to geotools GeometryCollector
    // to build as of collection of the right type (multipoint, multiline or multipolygon).
    GeometryCollector collector = new GeometryCollector();
    collector.setFactory(geom.getFactory());
    parts.stream().forEach(part -> collector.add(part));
    return collector.collect();
  }

  /**
   * Get the midpoint of the given geometry, or empty if the midpoint is not contained by the geometry.
   *
   * An alternative to {@link Geometry#getCentroid() } with better handling of lines.
   *
   * When given a {@link LineString} the point half way along the line is returned. In most cases this
   * is far more useful than the centroid which is not likely to be contained by the line (except for
   * straight lines).
   *
   * Not finding a midpoint is expected for many polygon shapes. This is because polygons use the geometric
   * centre of the geometry. For an L shaped geometry this is likely to be outside of the geometry.
   *
   * @param geom geometry to find the midpoint of
   * @return midpoint within geom or empty is one was not found
   */
  public static Optional<Point> getMidpoint(Geometry geom) {
    if (geom instanceof Point point) {
      return Optional.of(point);
    }
    if (geom instanceof LineString line) {
      // a line if very unlikely to contains its centroid (unless it's completely straight)
      // so we traverse along the line segments looking of the segemnt that contains the mid length,
      // and return the point from that location.
      double length = line.getLength();
      double midLength = length / 2;
      double travelled = 0;
      for (int i = 0; i < line.getNumPoints() - 1; i++) {
        LineSegment segment = new LineSegment(
            line.getPointN(i).getCoordinate(),
            line.getPointN(i + 1).getCoordinate()
        );
        double segmentLength = segment.getLength();
        if (travelled + segmentLength > midLength) {
          Coordinate mid = segment.pointAlong((midLength - travelled) / segmentLength);
          return Optional.of(new Point(
              geom.getFactory().getCoordinateSequenceFactory().create(new Coordinate[] {mid}),
              geom.getFactory()
          ));
        }

        travelled += segment.getLength();
      }
    }
    Point mid = geom.getCentroid();
    if (geom.contains(mid)) {
      return Optional.of(mid);
    }
    return Optional.empty();
  }

  private static Geometry emptyGeom(GeometryFamily family, GeometryFactory factory) {
    switch (family) {
      case PUNTAL:
        return factory.createPoint();
      case LINEAL:
        return factory.createLineString();
      case POLYGONAL:
        return factory.createPolygon();
      default:
        throw new RuntimeException("Unknown geometry family " + family);
    }
  }

}
