/*
 * 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.data.coverage;


import java.util.Collections;
import java.util.Iterator;
import java.util.LinkedList;
import java.util.List;
import java.util.Optional;
import java.util.function.BiFunction;
import java.util.function.Function;

import org.geotools.api.coverage.PointOutsideCoverageException;
import org.geotools.api.geometry.MismatchedDimensionException;
import org.geotools.api.geometry.Position;
import org.geotools.api.referencing.crs.CoordinateReferenceSystem;
import org.geotools.api.referencing.operation.TransformException;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.util.CoverageUtilities;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.referencing.CRS;
import org.geotools.referencing.CRS.AxisOrder;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.Point;

import it.geosolutions.jaiext.range.NoDataContainer;
import lombok.Getter;
import nz.org.riskscape.engine.SRIDSet;
import nz.org.riskscape.engine.coverage.TypedCoverage;
import nz.org.riskscape.engine.geo.GeometryFamily;
import nz.org.riskscape.engine.geo.GeometryUtils;
import nz.org.riskscape.engine.grid.FeatureGrid;
import nz.org.riskscape.engine.grid.FeatureGridCell;
import nz.org.riskscape.engine.relation.Relation;
import nz.org.riskscape.engine.types.Type;
import nz.org.riskscape.engine.types.Types;
import nz.org.riskscape.engine.util.Pair;

public class GridTypedCoverage implements TypedCoverage {

  private final CoverageValueProviderHelper helper;

  @Getter
  private final GridCoverage2D coverage;
  @Getter
  private final NoDataContainer noData;
  @Getter
  private final SRIDSet sridSet;

  private final AxisOrder axisOrder;

  public GridTypedCoverage(GridCoverage2D coverage, SRIDSet sridSet) {
    this.helper = new CoverageValueProviderHelper(new ReferencedEnvelope(coverage.getEnvelope()), sridSet);
    this.coverage = coverage;
    if (coverage instanceof GridCoverage2D) {
      this.noData = CoverageUtilities.getNoDataProperty(coverage);
    } else {
      this.noData = null;
    }
    this.sridSet = sridSet;
    // memoise the axis order for the CRS, we will need this for any cutting
    this.axisOrder = CRS.getAxisOrder(getCoordinateReferenceSystem());
  }

  @Override
  public Type getType() {
    return Types.FLOATING;
  }

  @Override
  public Object evaluate(Point inputPoint) {
    try {
      Point point = (Point)helper.reprojectIfNecessary(inputPoint);
      Position position = helper.pointToDirectPosition(point);
      if (! helper.contains(position)) {
        return null;
      }
      double[] stored;
      try {
        stored = coverage.evaluate(position, (double[]) null);
      } catch (PointOutsideCoverageException ex) {
        // due to rounding errors, it's possible that a bounds check (above) can succeed but the coverage can
        // still throw a PointOutsideCoverageException
        return null;
      }

      double gridValue = stored[0];
      if (isNoData(gridValue)) {
        return null;
      }

      return stored[0];

    } catch (CoverageValueProviderHelper.EmptyGeometryException e) {
      // You can't convert an empty point to a direct position as it has no position.
      return null;
    }
  }

  private boolean isNoData(double gridValue) {
    //Note that the geotools GeoTiffReader does not set the NODATA property if the nodata value is NaN
    //That's why we return null if noData is null and value is NaN
    return ((noData != null && noData.getAsRange().contains(gridValue))
        || (noData == null && Double.isNaN(gridValue))
        );
  }

  /**
   * Iterates over the grid cells that overlap the lookupGeom envelope and applies the mapper function.
   * Note that 'no data' cells are ignored.
   * @param mapper manipulates the overlapping cell and returns a result, or null to ignore the cell.
   */
  private <T> List<T> iterateCells(Geometry lookupGeom, BiFunction<FeatureGridCell, Double, T> mapper) {

    try {
      final Geometry geom = helper.reprojectIfNecessary(lookupGeom);
      FeatureGrid featureGrid;
      try {
        featureGrid = new FeatureGrid(geom, axisOrder, coverage.getGridGeometry());
      } catch (MismatchedDimensionException | TransformException e) {
        throw new RuntimeException(e);
      }
      List<T> hits = new LinkedList<>();
      Iterator<FeatureGridCell> cellIterator = featureGrid.cellIterator();

      while (cellIterator.hasNext()) {
        FeatureGridCell cell = cellIterator.next();

        // sample first - this is always cheaper than computing an intersection
        double[] sampled;
        try {
          sampled = coverage.evaluate(cell.getGridPosition(), (double[]) null);
        } catch (PointOutsideCoverageException e) {
          // NB need to investigate if this is just happening on the 'border' or more frequently
          continue;
        }

        if (isNoData(sampled[0])) {
          // no data in cell, we don't bother returning that
          continue;
        }
        T result = mapper.apply(cell, sampled[0]);
        if (result != null) {
          // we have a hit!
          hits.add(result);
        }
      }

      return hits;
    } catch (CoverageValueProviderHelper.EmptyGeometryException e) {
      // You can't convert an empty point to a direct position as it has no position.
      return Collections.emptyList();
    }
  }

  @Override
  public Optional<Function<Geometry, List<Pair<Geometry, Object>>>> getEvaluateIntersectionOp() {
    return Optional.of(geom -> evaluateIntersection(geom));
  }

  public List<Pair<Geometry, Object>> evaluateIntersection(Geometry lookupGeom) {
    int featureSrid = lookupGeom.getSRID();
    final GeometryFamily expectedFamily = GeometryFamily.from(lookupGeom);

    return iterateCells(lookupGeom, (cell, sampled) -> {
      // compute the intersection, throw away glancing intersections (they won't be the same family of geometry)
      Geometry geomCutByCell = GeometryUtils.removeNonFamilyMembers(cell.computeIntersection(), expectedFamily);

      // empty intersection means it didn't intersect.
      if (geomCutByCell.isEmpty()) {
        return null;
      }

      int gridSrid = sridSet.get(cell.getCRS());

      // reproject the cut geometry back to the CRS it originally came in with, if needed
      if (featureSrid != gridSrid) {
        geomCutByCell.setSRID(gridSrid);
        geomCutByCell = sridSet.reproject(geomCutByCell, featureSrid);
      }

      return Pair.of(geomCutByCell, sampled);
    });
  }

  @Override
  public Optional<Function<Geometry, List<? extends Object>>> getEvaluateIntersectingValuesOp() {
    return Optional.of(geom -> evaluateIntersectingValues(geom));
  }

  public List<Object> evaluateIntersectingValues(Geometry lookupGeom) {
    return iterateCells(lookupGeom,
            (cell, sampled) -> cell.intersects() ? sampled : null
    );
  }

  @Override
  public CoordinateReferenceSystem getCoordinateReferenceSystem() {
    return getCoverage().getCoordinateReferenceSystem();
  }

  @Override
  public String toString() {
    return String.format("%s[type=%s, coverage=%s]", getClass().getSimpleName(), getType(), coverage);
  }

  @Override
  public Optional<ReferencedEnvelope> getEnvelope() {
    return Optional.of(helper.getEnvelope());
  }

  @Override
  public Optional<Relation> asRelation() {
    if (coverage instanceof GridCoverage2D) {
      return Optional.of(GridCoverageRelation.create(this, coverage));
    } else {
      return TypedCoverage.super.asRelation();
    }
  }
}
