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

import static nz.org.riskscape.engine.Matchers.*;
import static org.hamcrest.Matchers.*;
import static org.junit.Assert.*;
import static org.mockito.Mockito.*;

import java.nio.file.Files;
import java.nio.file.Path;
import java.util.EnumSet;
import java.util.List;
import java.util.Optional;

import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridEnvelope2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.junit.After;
import org.junit.Before;
import org.junit.Test;
import org.locationtech.jts.geom.Geometry;

import lombok.extern.slf4j.Slf4j;
import nz.org.riskscape.engine.DummyFunction;
import nz.org.riskscape.engine.GeometryProblems;
import nz.org.riskscape.engine.Matchers;
import nz.org.riskscape.engine.ProjectTest;
import nz.org.riskscape.engine.TemporaryDirectoryTestHelper;
import nz.org.riskscape.engine.Tuple;
import nz.org.riskscape.engine.bind.Parameter;
import nz.org.riskscape.engine.data.Bookmark;
import nz.org.riskscape.engine.data.coverage.CoverageFileBookmarkResolver;
import nz.org.riskscape.engine.data.coverage.GridTypedCoverage;
import nz.org.riskscape.engine.gt.GeometryHelper;
import nz.org.riskscape.engine.problem.GeneralProblems;
import nz.org.riskscape.engine.raster.VectorToRaster;
import nz.org.riskscape.engine.resource.CreateHandle;
import nz.org.riskscape.engine.resource.FileCreateHandle;
import nz.org.riskscape.engine.rl.GeometryFunctions;
import nz.org.riskscape.engine.types.Nullable;
import nz.org.riskscape.engine.types.Referenced;
import nz.org.riskscape.engine.types.Struct;
import nz.org.riskscape.engine.types.Type;
import nz.org.riskscape.engine.types.TypeProblems;
import nz.org.riskscape.engine.types.Types;
import nz.org.riskscape.problem.ResultOrProblems;
import nz.org.riskscape.rl.ast.ExpressionProblems;
import nz.org.riskscape.rl.ast.StructDeclaration;

@Slf4j
public class GeoTiffFormatTest extends ProjectTest implements TemporaryDirectoryTestHelper {

  GeometryHelper geomHelper = GeometryHelper.nzTransverseMercatorXYWithOriginShift(project.getSridSet());
  GeometryHelper latLongHelper = GeometryHelper.latLong(project.getSridSet());
  GeometryHelper longLatHelper = GeometryHelper.longLat(project.getSridSet());

  GeoTiffFormat subject = new GeoTiffFormat();

  Struct inputType = Struct.of("value", Types.FLOATING, "geom", Referenced.of(Types.GEOMETRY, geomHelper.getCrs()));
  GeoTiffFormat.Options options = new GeoTiffFormat.Options();

  Path testDirectory;
  CreateHandle handle;
  GridTypedCoverage templateCoverage;

  CoverageFileBookmarkResolver coverageResolver = new CoverageFileBookmarkResolver(engine);
  DummyFunction myBounds = new DummyFunction("my_bounds", List.of()) {
    @Override
    public Type getReturnType() {
      return Types.GEOMETRY;
    }
  };

  @Before
  public void init() throws Exception {
    testDirectory = createTempDirectory("GeoTiffFormatTest");
    handle = new FileCreateHandle(null, Files.createTempFile(testDirectory, "raster", ".tif").toFile());

    GeometryFunctions geometryFunctions = new GeometryFunctions(engine);
    project.getFunctionSet().addAll(geometryFunctions.getFunctions());
    project.getFunctionSet().add(myBounds);
    myBounds.pickledLossValue = geomHelper.box(0, 0, 100, 100);

    options.gridResolution = Optional.of(10D);
    options.bounds = Optional.of(expressionParser.parse("my_bounds()"));

    // set up the coverage80x80_1x2grid function response
    templateCoverage = mock(GridTypedCoverage.class);
    // make an 80x40 env to be different to default 100x100
    ReferencedEnvelope env = geomHelper.envelope(
        geomHelper.toCoordinate(0, 0),
        geomHelper.toCoordinate(80, 40)
    );
    when(templateCoverage.getEnvelope()).thenReturn(Optional.of(env));
    GridEnvelope2D env2D = new GridEnvelope2D(0, 0, 80, 40);
    GridGeometry2D gridGeom2D = mock(GridGeometry2D.class);
    when(gridGeom2D.getGridRange2D()).thenReturn(env2D);
    GridCoverage2D gridCoverage2D = mock(GridCoverage2D.class);
    when(gridCoverage2D.getGridGeometry()).thenReturn(gridGeom2D);
    when(templateCoverage.getCoverage()).thenReturn(gridCoverage2D);
    when(templateCoverage.getCoordinateReferenceSystem()).thenReturn(geomHelper.getCrs());
  }

  @After
  public void cleanup() throws Exception {
    remove(testDirectory);
  }

  @Test
  public void hasExpectedDetails() {
    assertThat(subject.getId(), is("geotiff"));
    assertThat(subject.getExtension(), is("tif"));
    assertThat(subject.getMediaType(), is("image/tiff"));
    assertThat(subject.getCharacteristics(), is(EnumSet.noneOf(Format.Characteristics.class)));
    assertThat(subject.getWriterParameters().getBindsTo(), is(GeoTiffFormat.Options.class));
  }

  @Test
  public void writesFeaturesToRaster() throws Exception {
    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.box(0, 0, 10, 10)));
    writer.write(Tuple.ofValues(inputType, 20D, geomHelper.box(10, 0, 20, 10)));
    // this feature extends outside of the grid
    writer.write(Tuple.ofValues(inputType, 100D, geomHelper.box(90, 90, 110, 110)));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat(coverage.evaluate(geomHelper.point(1, 1)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(10, 1)), is(20D));
    assertThat(coverage.evaluate(geomHelper.point(11, 11)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(10.01, 10.01)), nullValue());

    assertThat(coverage.evaluate(geomHelper.point(99, 99)), is(100D));

    assertThat(coverage.getCoverage().getGridGeometry().getGridRange2D(), allOf(
        hasProperty("maxX", is(10D)),
        hasProperty("maxY", is(10D))
    ));
  }

  @Test
  public void canUseGridCoverageAsTemplate() throws Exception {
    // we remove the referencing from the input type. we don't need it when bounds is a grid coverage template
    inputType = Struct.of("value", Types.FLOATING, "geom", Types.GEOMETRY);
    options.template = Optional.of(templateCoverage);
    options.bounds = Optional.empty();
    options.gridResolution = Optional.empty();

    RiskscapeWriter writer = newWriter().getOrThrow();
    // NB: geometry will be written in template CRS rather than lat,long
    writer.write(Tuple.ofValues(inputType, 10D, latLongHelper.reproject(geomHelper.point(2.2, 2.2))));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat(coverage.evaluate(geomHelper.point(2.01, 2.01)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(2.99, 2.99)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(2.2, 4.1)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(3.2, 2.2)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(2.2, 1.9)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(1.9, 2.2)), nullValue());

    assertThat(coverage.getEnvelope(), isPresent(is(geomHelper.envelope(
        geomHelper.toCoordinate(0, 0),
        geomHelper.toCoordinate(80, 40)
    ))));
    assertThat(coverage.getCoverage().getGridGeometry().getGridRange2D(), allOf(
        hasProperty("maxX", is(80D)),
        hasProperty("maxY", is(40D))
    ));
  }

  @Test
  public void canUseGridCoverageTemplateWithGridResolutionOverride() throws Exception {
    // we remove the referencing from the input type. we don't need it when bounds is a grid coverage template
    inputType = Struct.of("value", Types.FLOATING, "geom", Types.GEOMETRY);
    options.template = Optional.of(templateCoverage);
    options.bounds = Optional.empty();
    // increase pixel size from 1 to 2
    options.gridResolution = Optional.of(2D);

    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 10D, latLongHelper.reproject(geomHelper.point(2.2, 2.2))));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat(coverage.evaluate(geomHelper.point(2.01, 2.01)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(2.99, 2.99)), is(10D));
    // pixel-size=2 so this is still the same pixel
    assertThat(coverage.evaluate(geomHelper.point(3.2, 2.2)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(2.2, 3.2)), is(10D));

    assertThat(coverage.evaluate(geomHelper.point(2.2, 4.1)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(2.2, 1.9)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(1.9, 2.2)), nullValue());

    assertThat(coverage.getEnvelope(), isPresent(is(geomHelper.envelope(
        geomHelper.toCoordinate(0, 0),
        geomHelper.toCoordinate(80, 40)
    ))));
    // pixels have now halved because they're twice as big
    assertThat(coverage.getCoverage().getGridGeometry().getGridRange2D(), allOf(
        hasProperty("maxX", is(40D)),
        hasProperty("maxY", is(20D))
    ));
  }

  @Test
  public void canUseGridCoverageTemplateWithBoundsOverride() throws Exception {
    // we remove the referencing from the input type. we don't need it when bounds is a grid coverage template
    inputType = Struct.of("value", Types.FLOATING, "geom", Types.GEOMETRY);
    options.template = Optional.of(templateCoverage);
    // use default bounds of 100x100, but grid resolution defaults to coverage (1m pixels)
    options.bounds = Optional.of(expressionParser.parse("my_bounds()"));
    options.gridResolution = Optional.empty();

    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 10D, latLongHelper.reproject(geomHelper.point(2.2, 2.2))));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat(coverage.evaluate(geomHelper.point(2.01, 2.01)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(2.99, 2.99)), is(10D));

    assertThat(coverage.evaluate(geomHelper.point(3.2, 2.2)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(2.2, 3.2)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(2.2, 4.1)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(2.2, 1.9)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(1.9, 2.2)), nullValue());

    assertThat(coverage.getEnvelope(), isPresent(is(geomHelper.envelope(
            geomHelper.toCoordinate(0, 0),
            geomHelper.toCoordinate(100, 100)
    ))));
    assertThat(coverage.getCoverage().getGridGeometry().getGridRange2D(), allOf(
            hasProperty("maxX", is(100D)),
            hasProperty("maxY", is(100D))
    ));
  }

  @Test
  public void canUseLongLatGridCoverageAsTemplate() throws Exception {
    // this is a giant coverage, where each pixel is one degree WGS84
    ReferencedEnvelope env = longLatHelper.envelope(
            longLatHelper.toCoordinate(0, 0),
            longLatHelper.toCoordinate(80, 40)
    );
    when(templateCoverage.getEnvelope()).thenReturn(Optional.of(env));
    when(templateCoverage.getCoordinateReferenceSystem()).thenReturn(longLatHelper.getCrs());

    inputType = Struct.of("value", Types.FLOATING, "geom", Types.GEOMETRY);
    options.template = Optional.of(templateCoverage);
    options.bounds = Optional.empty();
    options.gridResolution = Optional.empty();

    RiskscapeWriter writer = newWriter().getOrThrow();
    // NB: geometry will be written in long,lat
    writer.write(Tuple.ofValues(inputType, 10D, latLongHelper.point(12.2, 2.2)));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat(coverage.evaluate(latLongHelper.point(12.01, 2.01)), is(10D));
    assertThat(coverage.evaluate(latLongHelper.point(12.99, 2.99)), is(10D));
    assertThat(coverage.evaluate(longLatHelper.point(2.01, 12.01)), is(10D));
    assertThat(coverage.evaluate(longLatHelper.point(2.99, 12.99)), is(10D));

    assertThat(coverage.evaluate(latLongHelper.point(12.2, 4.1)), nullValue());
    assertThat(coverage.evaluate(latLongHelper.point(13.2, 2.2)), nullValue());
    assertThat(coverage.evaluate(latLongHelper.point(12.2, 1.9)), nullValue());
    assertThat(coverage.evaluate(latLongHelper.point(11.9, 2.2)), nullValue());

    assertThat(coverage.getEnvelope(), isPresent(is(env)));
    assertThat(coverage.getCoverage().getGridGeometry().getGridRange2D(), allOf(
            hasProperty("maxX", is(80D)),
            hasProperty("maxY", is(40D))
    ));
  }

  @Test
  public void canProduceTiffWithLatLongFeatures() throws Exception {
    inputType = Struct.of("value", Types.FLOATING, "geom", Referenced.of(Types.GEOMETRY, latLongHelper.getCrs()));
    myBounds.pickledLossValue = latLongHelper.box(12, 160, 12.1, 160.2);
    options.gridResolution = Optional.of(1112D);  // approx 0.01 degrees

    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 10D, latLongHelper.point(12.033, 160.054)));

    GridTypedCoverage coverage = closeThenRead(writer);

    // this is the one pixel that got burnt
    assertThat(coverage.evaluate(longLatHelper.point(160.054, 12.033)), is(10D));

    // the coverage got written in long/lat because coverages are always x/y
    assertThat(coverage.getEnvelope(), isPresent(is(longLatHelper.envelope(
        longLatHelper.toCoordinate(160, 12),
        longLatHelper.toCoordinate(160.2, 12.1)
    ))));
    assertThat(coverage.getCoverage().getGridGeometry().getGridRange2D(), allOf(
        hasProperty("maxX", is(20D)),
        hasProperty("maxY", is(10D))
    ));
  }

  @Test
  public void pointsOnlyBurnOnePixel1() throws Exception {
    // we want to ensure that if you save a point geometry it can only burn a single pixel
    // use a 1m grid resolution to keep things simple
    options.gridResolution = Optional.of(1D);
    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.point(20, 20)));
    // from the cell vertex put some point in horizontally
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.point(20.5, 20)));
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.point(20.9, 20)));
    // and now from vertex head down
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.point(20, 19.5)));
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.point(20, 19.1)));
    // some other points in the same pixel
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.point(20.5, 19.5)));

    GridTypedCoverage coverage = closeThenRead(writer);
    // this is the one pixel that got burnt
    assertThat(coverage.evaluate(geomHelper.point(20.5, 19.5)), is(10D));

    // not to the left
    assertThat(coverage.evaluate(geomHelper.point(19.5, 18.5)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(19.5, 19.5)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(19.5, 20.5)), nullValue());

    // not above or below
    assertThat(coverage.evaluate(geomHelper.point(20.5, 18.5)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(20.5, 20.5)), nullValue());

    // not to the right
    assertThat(coverage.evaluate(geomHelper.point(21.5, 18.5)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(21.5, 19.5)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(21.5, 20.5)), nullValue());
  }

  @Test
  public void canWritePolygonWithHoles() throws Exception {
    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.box(
        geomHelper.ring(
            geomHelper.toCoordinate(0, 0),
            geomHelper.toCoordinate(100, 0),
            geomHelper.toCoordinate(100, 100),
            geomHelper.toCoordinate(0, 100),
            geomHelper.toCoordinate(0, 0)
        ),
        geomHelper.ring(
            geomHelper.toCoordinate(19, 19),
            geomHelper.toCoordinate(81, 19),
            geomHelper.toCoordinate(81, 81),
            geomHelper.toCoordinate(19, 81),
            geomHelper.toCoordinate(19, 19)
        )
    )));

    GridTypedCoverage coverage = closeThenRead(writer);

    // sample around the shape
    assertThat(coverage.evaluate(geomHelper.point(15, 15)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(15, 85)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(85, 85)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(85, 15)), is(10D));

    // sample parts of the hole
    assertThat(coverage.evaluate(geomHelper.point(25, 25)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(25, 75)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(75, 75)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(75, 25)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(50, 50)), nullValue());
  }

  @Test
  public void boundsCanBeInDifferentCRS() throws Exception {
    myBounds.pickledLossValue = latLongHelper.reproject((Geometry)myBounds.pickledLossValue);
    // use overwrite strategy because otherwise the first 2 polygons overlap the same pixels slightly.
    // This is because the grid alignment has shifted slightly to match the reprojected bounds
    options.pixelStatistic = VectorToRaster.PixelStrategy.OVERWRITE;
    options.gridResolution = Optional.of(1D);

    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.box(0, 0, 10, 10)));
    writer.write(Tuple.ofValues(inputType, 20D, geomHelper.box(10, 0, 20, 10)));
    // this feature extends outside of the grid
    writer.write(Tuple.ofValues(inputType, 100D, geomHelper.box(90, 90, 110, 110)));

    GridTypedCoverage coverage = closeThenRead(writer);

    assertThat(coverage.evaluate(geomHelper.point(1.1, 1.1)), is(10D));
    assertThat(coverage.evaluate(geomHelper.point(10.1, 1.1)), is(20D));
    assertThat(coverage.evaluate(geomHelper.point(99.1, 99.1)), is(100D));
    // grid-resolution=1, so these points still fall outside the polygons, even though the
    // alignment of the coverage has changed slightly due to the bounds being in a different CRS
    assertThat(coverage.evaluate(geomHelper.point(11.1, 11.1)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(21.1, 1.1)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(89.5, 90.5)), nullValue());

    // the output GeoTIFF aligns to the bounds, which was reprojected from nz -> latlong -> nz
    // which shifts the GeoTIFF alignment slightly. So the polygons spill over into the next
    // cell slightly, because they're not perfectly aligned to the GeoTIFF grid anymore
    assertThat(coverage.evaluate(geomHelper.point(10.01, 10.01)), is(20D));

    assertThat(coverage.getCoverage().getGridGeometry().getGridRange2D(), allOf(
        hasProperty("maxX", is(100D)),
        hasProperty("maxY", is(101D))
    ));
  }

  @Test
  public void canSetPixelToOverwriteWithLastFeature() throws Exception {
    // NB: overwrite probably isn't that useful for users, but it's the old VectorToRaster behaviour and
    // covered here for completeness
    options.pixelStatistic = VectorToRaster.PixelStrategy.OVERWRITE;
    RiskscapeWriter writer = newWriter().getOrThrow();
    // resolution=10 so these are all the same pixel
    writer.write(Tuple.ofValues(inputType, 1D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 2D, geomHelper.box(1, 0, 2, 1)));
    writer.write(Tuple.ofValues(inputType, 3D, geomHelper.box(2, 0, 3, 1)));

    GridTypedCoverage coverage = closeThenRead(writer);
    // resolution=10 so these are all the same pixel
    assertThat(coverage.evaluate(geomHelper.point(0.5, 0.5)), is(3D));
    assertThat(coverage.evaluate(geomHelper.point(1.5, 0.5)), is(3D));
    assertThat(coverage.evaluate(geomHelper.point(2.5, 0.5)), is(3D));

    assertThat(coverage.getCoverage().getGridGeometry().getGridRange2D(), allOf(
        hasProperty("maxX", is(10D)),
        hasProperty("maxY", is(10D))
    ));
  }

  @Test
  public void canSetPixelToMaxValue() throws Exception {
    options.pixelStatistic = VectorToRaster.PixelStrategy.MAX;
    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 19.8D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 19.9D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 19.7D, geomHelper.box(0, 0, 1, 1)));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat((Double)coverage.evaluate(geomHelper.point(0.5, 0.5)), closeTo(19.9D, 0.00001));
  }

  @Test
  public void canSetPixelToMinValue() throws Exception {
    options.pixelStatistic = VectorToRaster.PixelStrategy.MIN;
    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 1D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 2D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 3D, geomHelper.box(0, 0, 1, 1)));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat(coverage.evaluate(geomHelper.point(0.5, 0.5)), is(1D));
  }

  @Test
  public void canSetPixelToSumOfValues() throws Exception {
    inputType = inputType.replace("value", Nullable.FLOATING);
    options.pixelStatistic = VectorToRaster.PixelStrategy.SUM;
    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 1.1D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 20.2D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 3.7D, geomHelper.box(0, 0, 1, 1)));
    // null values shouldn't mess up the total
    writer.write(Tuple.ofValues(inputType, null, geomHelper.box(0, 0, 1, 1)));
    // nor should Nan
    writer.write(Tuple.ofValues(inputType, Double.NaN, geomHelper.box(0, 0, 1, 1)));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat((Double) coverage.evaluate(geomHelper.point(0.5, 0.5)), closeTo(25D, 0.00001));
  }

  @Test
  public void canSetPixelToMeanOfValues() throws Exception {
    inputType = inputType.replace("value", Nullable.FLOATING);
    options.pixelStatistic = VectorToRaster.PixelStrategy.MEAN;
    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 1.5D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 19.25D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 3.25D, geomHelper.box(0, 0, 1, 1)));
    // null values should be ignored from the average
    writer.write(Tuple.ofValues(inputType, null, geomHelper.box(0, 0, 1, 1)));
    // nor should NaN
    writer.write(Tuple.ofValues(inputType, Double.NaN, geomHelper.box(0, 0, 1, 1)));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat((Double) coverage.evaluate(geomHelper.point(0.5, 0.5)), closeTo(8D, 0.0001));
  }

  @Test
  public void willWriteRasterWithOneMeterResolution() throws Exception {
    // same features as canSetPixelToOverwriteWithLastFeature but with one metre resolution
    // they won't overwrite
    options.pixelStatistic = VectorToRaster.PixelStrategy.OVERWRITE;
    options.gridResolution = Optional.of(1D);
    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 1D, geomHelper.box(0, 0, 1, 1)));
    writer.write(Tuple.ofValues(inputType, 2D, geomHelper.box(1, 0, 2, 1)));
    writer.write(Tuple.ofValues(inputType, 3D, geomHelper.box(2, 0, 3, 1)));
    writer.write(Tuple.ofValues(inputType, 99D, geomHelper.box(99, 0, 100, 1)));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat(coverage.evaluate(geomHelper.point(0.5, 0.5)), is(1D));
    assertThat(coverage.evaluate(geomHelper.point(1.5, 0.5)), is(2D));
    assertThat(coverage.evaluate(geomHelper.point(2.5, 0.5)), is(3D));

    assertThat(coverage.evaluate(geomHelper.point(98.99, 0.5)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(99.01, 0.5)), is(99D));
    assertThat(coverage.evaluate(geomHelper.point(99.99, 0.5)), is(99D));

    assertThat(coverage.getCoverage().getGridGeometry().getGridRange2D(), allOf(
        hasProperty("maxX", is(100D)),
        hasProperty("maxY", is(100D))
    ));
  }

  @Test
  public void canDrawLineFeatures() throws Exception {
    options.gridResolution = Optional.of(1D);

    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 1D, geomHelper.line(1, 0.5, 3, 0.5)));
    writer.write(Tuple.ofValues(inputType, 2D, geomHelper.box(10, 10, 12, 12)));

    GridTypedCoverage coverage = closeThenRead(writer);
    assertThat(coverage.evaluate(geomHelper.point(0.99, 1)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(1.01, 0.5)), is(1D));
    assertThat(coverage.evaluate(geomHelper.point(2.01, 0.5)), is(1D));
    assertThat(coverage.evaluate(geomHelper.point(2.99, 0.5)), is(1D));
    // NB: the line glances x=3.0,y=0.5 but technically doesn't intersect the pixel
    assertThat(coverage.evaluate(geomHelper.point(3.01, 0.5)), nullValue());

    assertThat(coverage.evaluate(geomHelper.point(10.5, 10.5)), is(2D));
    assertThat(coverage.evaluate(geomHelper.point(11.9, 11.9)), is(2D));
    assertThat(coverage.evaluate(geomHelper.point(10.5, 11.5)), is(2D));
    assertThat(coverage.evaluate(geomHelper.point(11.5, 10.5)), is(2D));

    // the polygon ends at 12,12 so it technically glances these pixels but doesn't intersect
    assertThat(coverage.evaluate(geomHelper.point(11.01, 12.01)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(12.01, 12.01)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(12.01, 11.01)), nullValue());
  }

  @Test
  public void canAddFeaturesAtTheirCentroid() throws Exception {
    // using the centroid when creating the raster could be a general use case
    // so we have a specific test to ensure it works.
    options.geometry = Optional.of(expressionParser.parse("centroid(geom)"));
    options.gridResolution = Optional.of(1D);

    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 1D, geomHelper.box(0, 0, 3, 1)));
    writer.write(Tuple.ofValues(inputType, 2D, geomHelper.box(1, 0, 4, 1)));

    GridTypedCoverage coverage = closeThenRead(writer);
    // these are the centroid coordinates
    assertThat(coverage.evaluate(geomHelper.point(1.5, 0.5)), is(1D));
    assertThat(coverage.evaluate(geomHelper.point(2.5, 0.5)), is(2D));
    // these pixels intersected the original polygons, but not the centroids
    assertThat(coverage.evaluate(geomHelper.point(0.5, 0.5)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(3.5, 0.5)), nullValue());

    // pixels outside both the polygon and centroid
    assertThat(coverage.evaluate(geomHelper.point(1.5, 1.01)), nullValue());
    assertThat(coverage.evaluate(geomHelper.point(2.5, 1.01)), nullValue());
  }

  @Test
  public void warningForFeatureOutOfBounds() throws Exception {
    RiskscapeWriter writer = newWriter().getOrThrow();
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.point(-100, 100)));
    writer.write(Tuple.ofValues(inputType, 10D, geomHelper.point(-200, 100)));

    closeThenRead(writer);

    assertThat(sunkProblems, contains(
        GeoTiffFormat.PROBLEMS.skippedFeaturesOutOfBounds(2, writer.getStoredAt())
    ));
  }

  @Test
  public void warningForFeatureValueIsNullOrNaN() throws Exception {
    RiskscapeWriter writer = newWriter().getOrThrow();
    // we don't count null/NaN values or null geometries
    writer.write(Tuple.ofValues(inputType, Double.NaN, geomHelper.point(100, 100)));
    writer.write(Tuple.ofValues(inputType, null, geomHelper.point(-200, 100)));
    writer.write(Tuple.ofValues(inputType, 10D, null));

    closeThenRead(writer);

    assertThat(sunkProblems, contains(
        GeoTiffFormat.PROBLEMS.skippedFeaturesNullOrNan(3, writer.getStoredAt())
    ));
  }

  @Test
  public void failsToGetWriterWhenMissingGridResolution() {
    options.gridResolution = Optional.empty();

    // NB might come back and fix this after some more deck-chair shuffling
    assertThat(newWriter(), Matchers.failedResult(
        is(GeneralProblems.required("grid-resolution", Parameter.class))
    ));
  }

  @Test
  public void failsToGetWriterWhenMissingBounds() {
    options.bounds = Optional.empty();

    // NB might come back and fix this after some more deck-chair shuffling
    assertThat(newWriter(), Matchers.failedResult(
        is(GeneralProblems.required("bounds", Parameter.class))
    ));
  }

  @Test
  public void failsToGetWriterNoNumericValue() {
    inputType = inputType.replace("value", Types.TEXT);
    assertThat(newWriter(), Matchers.failedResult(
        is(TypeProblems.get().structMustHaveMemberType(Types.FLOATING, inputType))
    ));
  }

  @Test
  public void failsToGetWriterValueExprNotNumeric() {
    inputType = inputType.replace("value", Types.TEXT);
    options.value = Optional.of(expressionParser.parse("value"));
    assertThat(newWriter(), Matchers.failedResult(
        is(TypeProblems.get().requiresOneOf("value", List.of(Types.FLOATING, Types.INTEGER), Types.TEXT))
    ));
  }

  @Test
  public void failsToGetWriterNoGeometryValue() {
    inputType = inputType.replace("geom", Types.TEXT);
    assertThat(newWriter(), Matchers.failedResult(
        is(TypeProblems.get().structMustHaveMemberType(Types.GEOMETRY, inputType))
    ));
  }

  @Test
  public void failsToGetWriterGeometryExprNotGeometry() {
    inputType = inputType.replace("geom", Types.TEXT);
    options.geometry = Optional.of(expressionParser.parse("geom"));
    assertThat(newWriter(), Matchers.failedResult(
        is(TypeProblems.get().mismatch("geometry", Types.GEOMETRY, Types.TEXT))
    ));
  }

  @Test
  public void failsToGetWriterIfGeometryNotReferenced() {
    inputType = inputType.replace("geom", Types.GEOMETRY);
    options.geometry = Optional.of(expressionParser.parse("geom"));
    assertThat(newWriter(), Matchers.failedResult(
        is(GeometryProblems.get().notReferenced(Types.GEOMETRY))
    ));
  }

  @Test
  public void failsToGetWriterIfBoundsNotConstant() {
    options.bounds = Optional.of(expressionParser.parse("the_bounds"));
    assertThat(newWriter(), Matchers.failedResult(
        is(ExpressionProblems.get().constantRequired(options.bounds.get()))
    ));
  }

  @Test
  public void failsToGetWriterIfBoundsNotGeometry() {
    options.bounds = Optional.of(expressionParser.parse("'south island'"));
    assertThat(newWriter(), Matchers.failedResult(
        is(TypeProblems.get().mismatch("bounds", Types.GEOMETRY, Types.TEXT))
    ));
  }

  @Test
  public void failsToGetWriterIfGridResolutionNotProvided() {
    options.gridResolution = Optional.empty();
    assertThat(newWriter(), Matchers.failedResult(
      hasAncestorProblem(is(GeneralProblems.required("grid-resolution", Parameter.class)))
    ));
  }

  @Test
  public void failsToGetWriterIfBoundsNotProvided() {
    options.bounds = Optional.empty();
    assertThat(newWriter(), Matchers.failedResult(
            hasAncestorProblem(is(GeneralProblems.required("bounds", Parameter.class)))
    ));
  }

  @Test
  public void failsToGetWriterIfGridDimensionsTooBig() {
    options.gridResolution = Optional.of(0.5D);
    options.bounds = Optional.of(expressionParser.parse("my_bounds()"));
    myBounds.pickledLossValue = geomHelper.box(0, 0, 1_000_000, 1_000_000);
    assertThat(newWriter(), Matchers.failedResult(
        is(VectorToRaster.PROBLEMS.gridDimensionsTooBig(2_000_000, 2_000_000, Integer.MAX_VALUE)
            .withChildren(GeoTiffFormat.PROBLEMS.dimensionsTip()))
    ));
  }


  private ResultOrProblems<? extends RiskscapeWriter> newWriter() {
    return options.validate(inputType, executionContext).flatMap(validatedOptions -> newWriter(validatedOptions));
  }

  private ResultOrProblems<? extends RiskscapeWriter> newWriter(GeoTiffFormat.Options opts) {
    return subject.getWriterConstructor(opts)
          .get()
          .newWriter(executionContext, inputType, handle, StructDeclaration.EMPTY);
  }

  private GridTypedCoverage closeThenRead(RiskscapeWriter writer) {
    try {
      writer.close();
      Bookmark bm = Bookmark.fromURI(writer.getStoredAt()).withFormat("geotiff");
      log.info("About to read {}", bm.getLocation());
      return coverageResolver.resolve(bm, executionContext)
          .map(resolved -> resolved.getData(GridTypedCoverage.class))
          .map(result -> result.get())
          .orElseThrow();

    } catch (Throwable t) {
      throw new RuntimeException(t);
    }
  }
}
