org.opengis.referencing.datum.PixelInCell Java Examples

/**
 * Tests {@link GridDerivation#subgrid(Envelope, double...)} using only the
 * {@link GridExtent} result provided by {@link GridDerivation#getIntersection()}.
 */
@Test
public void testSubExtent() {
    GeneralEnvelope envelope = new GeneralEnvelope(HardCodedCRS.WGS84_3D);
    envelope.setRange(0, -80, 120);
    envelope.setRange(1, -12,  21);
    envelope.setRange(2,  10,  25);
    final MathTransform gridToCRS = MathTransforms.linear(new Matrix4(
            0,   0.5, 0,  -90,
            0.5, 0,   0, -180,
            0,   0,   2,    3,
            0,   0,   0,    1));
    final GridGeometry grid = new GridGeometry(PixelInCell.CELL_CORNER, gridToCRS, envelope, GridRoundingMode.NEAREST);
    assertExtentEquals(
            new long[] {336,  20,  4},
            new long[] {401, 419, 10}, grid.getExtent());
    /*
     * Set the region of interest as a two-dimensional envelope. The vertical dimension is omitted.
     * The result should be that all grid indices in the vertical dimension are kept unchanged.
     */
    envelope = new GeneralEnvelope(HardCodedCRS.WGS84);
    envelope.setRange(0, -70.001, +80.002);
    envelope.setRange(1,   4.997,  15.003);
    assertExtentEquals(new long[] {370,  40,  4},
                       new long[] {389, 339, 10}, grid.derive().subgrid(envelope).getIntersection());
}
 

/**
 * Tests the construction from a geospatial envelope.
 * The "grid to CRS" transform is explicitly given.
 */
@Test
public void testFromGeospatialEnvelope() {
    final GeneralEnvelope envelope = new GeneralEnvelope(HardCodedCRS.WGS84_φλ);
    envelope.setRange(0, -70.001, +80.002);
    envelope.setRange(1,   4.997,  15.003);
    final MathTransform gridToCRS = MathTransforms.linear(new Matrix3(
        0,   0.5, -90,
        0.5, 0,  -180,
        0,   0,     1));
    final GridGeometry grid = new GridGeometry(PixelInCell.CELL_CORNER, gridToCRS, envelope, GridRoundingMode.NEAREST);
    assertExtentEquals(
            new long[] {370, 40},
            new long[] {389, 339}, grid.getExtent());
    assertEnvelopeEquals(new GeneralEnvelope(
            new double[] {-70,  5},
            new double[] {+80, 15}), grid.getEnvelope(), STRICT);
    assertArrayEquals("resolution", new double[] {0.5, 0.5}, grid.getResolution(false), STRICT);
    assertMatrixEquals("gridToCRS", new Matrix3(
            0,   0.5, -89.75,
            0.5, 0,  -179.75,
            0,   0,     1), MathTransforms.getMatrix(grid.getGridToCRS(PixelInCell.CELL_CENTER)), STRICT);
}
 

/**
 * Tests the {@link Types#forCodeName(Class, String, boolean)} method.
 */
@Test
public void testForCodeName() {
    assertSame(ImagingCondition.SEMI_DARKNESS, Types.forCodeName(ImagingCondition.class, "SEMI_DARKNESS", false));
    assertSame(ImagingCondition.SEMI_DARKNESS, Types.forCodeName(ImagingCondition.class, "SEMIDARKNESS",  false));
    assertSame(ImagingCondition.SEMI_DARKNESS, Types.forCodeName(ImagingCondition.class, "semi darkness", false));
    assertSame(ImagingCondition.SEMI_DARKNESS, Types.forCodeName(ImagingCondition.class, "semi-darkness", false));
    assertNull(Types.forCodeName(ImagingCondition.class, "darkness", false));

    assertSame(PixelInCell.CELL_CORNER, Types.forCodeName(PixelInCell.class, "cell corner", false));
    assertSame(PixelInCell.CELL_CORNER, Types.forCodeName(PixelInCell.class, "cellCorner",  false));
    assertSame(PixelInCell.CELL_CENTER, Types.forCodeName(PixelInCell.class, "cell center", false));
    assertSame(PixelInCell.CELL_CENTER, Types.forCodeName(PixelInCell.class, "cellCenter",  false));

    if (PENDING_NEXT_GEOAPI_RELEASE) {
        assertSame(PixelInCell.CELL_CENTER, Types.forCodeName(PixelInCell.class, "cell centre", false));
        assertSame(PixelInCell.CELL_CENTER, Types.forCodeName(PixelInCell.class, "cellCentre",  false));
    }
}
 

/**
 * Tests construction with a non-linear component in the transform.
 */
@Test
public void testNonLinear() {
    final GridExtent extent = new GridExtent(
            new DimensionNameType[] {
                DimensionNameType.COLUMN,
                DimensionNameType.ROW,
                DimensionNameType.VERTICAL,
                DimensionNameType.TIME
            },
            new long[] {  0,   0, 2, 6},
            new long[] {100, 200, 3, 9}, false);
    final MathTransform horizontal = MathTransforms.linear(new Matrix3(
            0.5, 0,    12,
            0,   0.25, -2,
            0,   0,     1));
    final MathTransform vertical  = MathTransforms.interpolate(null, new double[] {1, 2, 4, 10});
    final MathTransform temporal  = MathTransforms.linear(3600, 60);
    final MathTransform gridToCRS = MathTransforms.compound(horizontal, vertical, temporal);
    final GridGeometry  grid      = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, null);
    assertArrayEquals("resolution", new double[] {0.5, 0.25,        6.0, 3600}, grid.getResolution(true),  STRICT);
    assertArrayEquals("resolution", new double[] {0.5, 0.25, Double.NaN, 3600}, grid.getResolution(false), STRICT);
    assertFalse("isConversionLinear", grid.isConversionLinear(0, 1, 2, 3));
    assertTrue ("isConversionLinear", grid.isConversionLinear(0, 1,    3));
}
 

/**
 * Checks that wraparound is well applied when using {@link GridDerivation#slice(DirectPosition)}.
 */
@Test
public void testSliceWithWrapAround() {
    final GridGeometry base = new GridGeometry(
            PixelInCell.CELL_CORNER,
            new AffineTransform2D(-0.02, 0, 0, 0.1, 55, 172),
            new Envelope2D(HardCodedCRS.WGS84_φλ, 42, 172, 13, 51),
            GridRoundingMode.NEAREST);

    final GridGeometry expectedResult = base.derive()
            .slice(new DirectPosition2D(51, 187))
            .build();

    final GridGeometry fromWrapAround = base.derive()
            .slice(new DirectPosition2D(51, -173))
            .build();

    assertEquals("Slice with wrap-around", expectedResult, fromWrapAround);
    assertBetween("Wrapped Y coordinate",
                  base.envelope.getMinimum(1),
                  base.envelope.getMaximum(1),
                  fromWrapAround.envelope.getMedian(1));
}
 

/**
 * Tests the adjustment done for pixel center in {@link GridGeometry#GridGeometry(GridGeometry, GridExtent, MathTransform)}
 * constructor. This test depends on {@link GridGeometry#derive()}, which will (indirectly) invoke the constructor to test.
 * We check envelopes as a more intuitive way to verify consistency than inspecting the math transforms.
 */
@Test
public void testFromOtherDefinedAtCenter() {
    GridExtent extent = new GridExtent(126, 197);
    GridGeometry grid = new GridGeometry(extent, PixelInCell.CELL_CENTER, MathTransforms.identity(2), HardCodedCRS.WGS84);
    GeneralEnvelope expected = new GeneralEnvelope(new double[] {-0.5, -0.5}, new double[] {125.5, 196.5});
    assertEnvelopeEquals(expected, grid.getEnvelope(), STRICT);
    /*
     * Derive a new grid geometry with 10×10 times more cells. The geographic area should be unchanged.
     */
    extent = extent.resize(1260, 1970);
    grid = grid.derive().resize(extent, 0.1, 0.1).build();
    assertEnvelopeEquals(expected, grid.getEnvelope(), STRICT);
    /*
     * If we create a grid geometry with identical properties, the envelope computed by that grid geometry would
     * be different than the envelope computed above if the "grid to CRS" transform is not correctly adjusted.
     */
    final GridGeometry alternative = new GridGeometry(grid.getExtent(), PixelInCell.CELL_CENTER,
             grid.getGridToCRS(PixelInCell.CELL_CENTER), grid.getCoordinateReferenceSystem());
    assertEnvelopeEquals(expected, alternative.getEnvelope(), STRICT);
}
 

/**
 * Tests the {@link GridGeometry#GridGeometry(GridGeometry, GridExtent, MathTransform)} constructor.
 * The math transform used for this test map to pixel corners.
 *
 * @throws TransformException if an error occurred while using the "grid to CRS" transform.
 */
@Test
public void testFromOtherDefinedAtCorner() throws TransformException {
    long[]        low       = new long[] {  1,   3, 2};
    long[]        high      = new long[] {101, 203, 4};
    GridExtent    extent    = new GridExtent(null, low, high, false);
    MathTransform gridToCRS = MathTransforms.translation(5, 7, 8);
    GridGeometry  grid      = new GridGeometry(extent, PixelInCell.CELL_CORNER, gridToCRS, null);

    low    = new long[] { 11,  35, 20};
    high   = new long[] {120, 250, 39};
    extent = new GridExtent(null, low, high, false);
    grid   = new GridGeometry(grid, extent, MathTransforms.scale(2, 1, 3));
    assertSame(extent, grid.getExtent());
    assertMatrixEquals("gridToCRS", new Matrix4(
            2, 0, 0, 5,
            0, 1, 0, 7,     // Combination of above scales (diagonal) and translation (last column).
            0, 0, 3, 8,
            0, 0, 0, 1), MathTransforms.getMatrix(grid.getGridToCRS(PixelInCell.CELL_CORNER)), STRICT);
}
 

/**
 * Tests deriving a grid geometry with an envelope crossing the antimeridian.
 */
@Test
public void testSubgridCrossingAntiMeridian() {
    final GridGeometry grid = new GridGeometry(
            new GridExtent(200, 180), PixelInCell.CELL_CORNER,
            MathTransforms.linear(new Matrix3(
                    1,  0, 80,
                    0, -1, 90,
                    0,  0,  1)), HardCodedCRS.WGS84);

    final GeneralEnvelope areaOfInterest = new GeneralEnvelope(HardCodedCRS.WGS84);
    areaOfInterest.setRange(0, 140, -179);                 // Cross anti-meridian.
    areaOfInterest.setRange(1, -90,   90);

    final GridGeometry subgrid = grid.derive().subgrid(areaOfInterest).build();
    final Envelope subEnv = subgrid.getEnvelope();
    areaOfInterest.setRange(0, 140, 181);
    assertEnvelopeEquals(areaOfInterest, subEnv);
}
 

/**
 * Verifies the exception thrown when we specify an envelope outside the grid extent.
 */
@Test
public void testOutsideDomain() {
    final GridGeometry grid = new GridGeometry(
            new GridExtent(10, 20), PixelInCell.CELL_CORNER,
            MathTransforms.linear(new Matrix3(
                    2, 0, 0,
                    0, 2, 0,
                    0, 0, 1)), HardCodedCRS.WGS84);

    final GeneralEnvelope areaOfInterest = new GeneralEnvelope(HardCodedCRS.WGS84);
    areaOfInterest.setRange(0, 60, 85);
    areaOfInterest.setRange(1, 15, 30);
    try {
        grid.derive().subgrid(areaOfInterest);
        fail("Should not have accepted the given AOI.");
    } catch (DisjointExtentException e) {
        assertNotNull(e.getMessage());
    }
}
 

/**
 * Tests deriving a grid geometry from an area of interest overlapping the grid in such a way
 * that we have to overlap the AOI to the full grid extent. Illustration:
 *
 * {@preformat text
 *                  ┌────────────────────────────────────────────┐
 *                  │             Domain of validity             │
 *                  └────────────────────────────────────────────┘
 *   ┌────────────────────┐                                ┌─────
 *   │  Area of interest  │                                │  AOI
 *   └────────────────────┘                                └─────
 *    ↖………………………………………………………360° period……………………………………………………↗︎
 * }
 */
@Test
public void testAreaOfInterestExpansion() {
    final GridGeometry grid = new GridGeometry(
            new GridExtent(340, 140), PixelInCell.CELL_CORNER,
            MathTransforms.linear(new Matrix3(
                    1,  0, 5,
                    0, -1, 70,
                    0,  0,  1)), HardCodedCRS.WGS84);

    final GeneralEnvelope areaOfInterest = new GeneralEnvelope(HardCodedCRS.WGS84);
    areaOfInterest.setRange(0, -30,  40);
    areaOfInterest.setRange(1, -60,  60);

    final GeneralEnvelope expected = new GeneralEnvelope(HardCodedCRS.WGS84);
    expected.setRange(0,   5, 345);
    expected.setRange(1, -60,  60);

    GridGeometry subgrid = grid.derive().subgrid(areaOfInterest).build();
    assertEnvelopeEquals(expected, subgrid.getEnvelope());
}
 

/**
 * Tests construction with an identity transform mapping pixel center.
 * This results a 0.5 pixel shifts in the "real world" envelope.
 */
@Test
public void testFromPixelCenter() {
    final long[]        low      = new long[] { 0,   0, 2};
    final long[]        high     = new long[] {99, 199, 4};
    final GridExtent    extent   = new GridExtent(null, low, high, true);
    final MathTransform identity = MathTransforms.identity(3);
    final GridGeometry  grid     = new GridGeometry(extent, PixelInCell.CELL_CENTER, identity, null);
    /*
     * Verify properties that should be stored "as-is".
     */
    final MathTransform trCenter = grid.getGridToCRS(PixelInCell.CELL_CENTER);
    assertSame("gridToCRS", identity, trCenter);
    assertExtentEquals(low, high, grid.getExtent());
    /*
     * Verify computed math transform.
     */
    final MathTransform trCorner = grid.getGridToCRS(PixelInCell.CELL_CORNER);
    assertNotSame(trCenter, trCorner);
    assertFalse ("gridToCRS.isIdentity",          trCorner.isIdentity());
    assertEquals("gridToCRS.sourceDimensions", 3, trCorner.getSourceDimensions());
    assertEquals("gridToCRS.targetDimensions", 3, trCorner.getTargetDimensions());
    assertMatrixEquals("gridToCRS", new Matrix4(
            1, 0, 0, -0.5,
            0, 1, 0, -0.5,
            0, 0, 1, -0.5,
            0, 0, 0,  1), MathTransforms.getMatrix(trCorner), STRICT);
    /*
     * Verify the envelope, which should have been computed using the math transform shifted by 0.5.
     */
    assertEnvelopeEquals(new GeneralEnvelope(
            new double[] {-0.5,  -0.5, 1.5},
            new double[] {99.5, 199.5, 4.5}), grid.getEnvelope(), STRICT);
    /*
     * Verify other computed properties.
     */
    assertArrayEquals("resolution", new double[] {1, 1, 1}, grid.getResolution(false), STRICT);
    assertTrue("isConversionLinear", grid.isConversionLinear(0, 1, 2));
}
 

/**
 * Tests {@link GridDerivation#subgrid(Envelope, double...)} with a non-linear "grid to CRS" transform.
 */
@Test
@DependsOnMethod("testSubExtent")
public void testSubExtentNonLinear() {
    final GridExtent extent = new GridExtent(
            new DimensionNameType[] {
                DimensionNameType.COLUMN,
                DimensionNameType.ROW,
                DimensionNameType.VERTICAL
            },
            new long[] {  0,  0, 2},
            new long[] {180, 90, 5}, false);
    final MathTransform linear = MathTransforms.linear(new Matrix4(
            2, 0, 0, -180,
            0, 2, 0,  -90,
            0, 0, 5,   10,
            0, 0, 0,    1));
    final MathTransform latitude  = MathTransforms.interpolate(new double[] {0, 20, 50, 70, 90}, new double[] {-90, -45, 0, 45, 90});
    final MathTransform gridToCRS = MathTransforms.concatenate(linear, MathTransforms.passThrough(1, latitude, 1));
    final GridGeometry  grid      = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, HardCodedCRS.WGS84_3D);
    /*
     * Following tests is similar to the one executed in testSubExtent(). Expected values are only
     * anti-regression values, except the vertical range which is expected to cover all cells. The
     * main purpose of this test is to verify that TransformSeparator has been able to extract the
     * two-dimensional transform despite its non-linear component.
     */
    final GeneralEnvelope envelope = new GeneralEnvelope(HardCodedCRS.WGS84);
    envelope.setRange(0, -70.001, +80.002);
    envelope.setRange(1,  -4.997,  15.003);
    final GridExtent actual = grid.derive().subgrid(envelope).getIntersection();
    assertEquals(extent.getAxisType(0), actual.getAxisType(0));
    assertExtentEquals(new long[] { 56, 69, 2},
                       new long[] {130, 73, 4}, actual);
}
 

/**
 * Tests {@link GridDerivation#slice(DirectPosition)}.
 */
@Test
public void testSlice() {
    final GridGeometry grid = new GridGeometry(
            new GridExtent(null, new long[] {336, 20, 4}, new long[] {401, 419, 10}, true),
            PixelInCell.CELL_CORNER, MathTransforms.linear(new Matrix4(
                    0,   0.5, 0,  -90,
                    0.5, 0,   0, -180,
                    0,   0,   2,    3,
                    0,   0,   0,    1)), HardCodedCRS.WGS84_3D);
    /*
     * There is two ways to ask for a slice. The first way is to set some coordinates to NaN.
     */
    GridGeometry slice = grid.derive().slice(new GeneralDirectPosition(Double.NaN, Double.NaN, 15)).build();
    assertNotSame(grid, slice);
    assertSame("gridToCRS", grid.gridToCRS, slice.gridToCRS);
    final long[] expectedLow  = {336,  20, 6};
    final long[] expectedHigh = {401, 419, 6};
    assertExtentEquals(expectedLow, expectedHigh, slice.getExtent());
    /*
     * Same test, but using a one-dimensional slice point instead than NaN values.
     * Opportunistically use different units for testing conversions.
     */
    GeneralDirectPosition p = new GeneralDirectPosition(HardCodedCRS.ELLIPSOIDAL_HEIGHT_cm);
    p.setOrdinate(0, 1500);
    slice = grid.derive().slice(p).build();
    assertNotSame(grid, slice);
    assertSame("gridToCRS", grid.gridToCRS, slice.gridToCRS);
    assertExtentEquals(expectedLow, expectedHigh, slice.getExtent());
}
 

/**
 * Tests deriving a grid geometry from an area of interest shifted by 360° before or after the grid valid area.
 */
@Test
public void testSubgridFromShiftedAOI() {
    final GridGeometry grid = new GridGeometry(
            new GridExtent(20, 140), PixelInCell.CELL_CORNER,
            MathTransforms.linear(new Matrix3(
                    1,  0, 80,
                    0, -1, 70,
                    0,  0,  1)), HardCodedCRS.WGS84);

    final GeneralEnvelope areaOfInterest = new GeneralEnvelope(HardCodedCRS.WGS84);
    areaOfInterest.setRange(0,  70,  90);
    areaOfInterest.setRange(1, -80,  60);

    final GeneralEnvelope expected = new GeneralEnvelope(HardCodedCRS.WGS84);
    expected.setRange(0,  80,  90);
    expected.setRange(1, -70,  60);
    GridGeometry subgrid;
    /*
     * Area of interest of the left side.
     */
    areaOfInterest.setRange(0, -290, -270);                    // [70 … 90] - 360
    subgrid = grid.derive().subgrid(areaOfInterest).build();
    assertEnvelopeEquals(expected, subgrid.getEnvelope());
    /*
     * Area of interest on the right side.
     */
    areaOfInterest.setRange(0, 430, 450);                      // [70 … 90] + 360
    subgrid = grid.derive().subgrid(areaOfInterest).build();
    assertEnvelopeEquals(expected, subgrid.getEnvelope());
}
 

/**
 * Tests construction with an identity transform mapping pixel corner.
 */
@Test
public void testFromPixelCorner() {
    final long[]         low     = new long[] {100, 300, 3, 6};
    final long[]         high    = new long[] {200, 400, 4, 7};
    final GridExtent    extent   = new GridExtent(null, low, high, true);
    final MathTransform identity = MathTransforms.identity(4);
    final GridGeometry  grid     = new GridGeometry(extent, PixelInCell.CELL_CORNER, identity, null);
    /*
     * Verify properties that should be stored "as-is".
     */
    final MathTransform trCorner = grid.getGridToCRS(PixelInCell.CELL_CORNER);
    assertSame("gridToCRS", identity, trCorner);
    assertExtentEquals(low, high, grid.getExtent());
    /*
     * Verify computed math transform.
     */
    final MathTransform trCenter = grid.getGridToCRS(PixelInCell.CELL_CENTER);
    assertNotSame(trCenter, trCorner);
    assertFalse ("gridToCRS.isIdentity",          trCenter.isIdentity());
    assertEquals("gridToCRS.sourceDimensions", 4, trCenter.getSourceDimensions());
    assertEquals("gridToCRS.targetDimensions", 4, trCenter.getTargetDimensions());
    assertMatrixEquals("gridToCRS", Matrices.create(5, 5, new double[] {
            1, 0, 0, 0, 0.5,
            0, 1, 0, 0, 0.5,
            0, 0, 1, 0, 0.5,
            0, 0, 0, 1, 0.5,
            0, 0, 0, 0, 1}), MathTransforms.getMatrix(trCenter), STRICT);
    /*
     * Verify the envelope, which should have been computed using the given math transform as-is.
     */
    assertEnvelopeEquals(new GeneralEnvelope(
            new double[] {100, 300, 3, 6},
            new double[] {201, 401, 5, 8}), grid.getEnvelope(), STRICT);
    /*
     * Verify other computed properties.
     */
    assertArrayEquals("resolution", new double[] {1, 1, 1, 1}, grid.getResolution(false), STRICT);
    assertTrue("isConversionLinear", grid.isConversionLinear(0, 1, 2, 3));
}
 

/**
 * Tests the construction from grid geometries having a linear "grid to CRS" conversion.
 *
 * @throws TransformException if an error occurred while computing the grid geometry.
 */
@Test
public void testSubgridFromOtherGrid() throws TransformException {
    GridGeometry   source = grid(  10,   -20,  110,  180, 100, -300);     // Envelope x: [1200 … 11300]   y: [-53800 … 6500]
    GridGeometry   target = grid(2000, -1000, 9000, 8000,   2,   -1);     // Envelope x: [4200 … 18202]   y: [ -7501 … 1500]
    GridDerivation change = target.derive().subgrid(source);              // Envelope x: [4200 … 11300]   y: [ -7501 … 1500]
    GridExtent     extent = change.getIntersection();
    GridExtentTest.assertExtentEquals(extent, 0,  2000, 5549);            // Subrange of target extent.
    GridExtentTest.assertExtentEquals(extent, 1, -1000, 8000);
    assertArrayEquals("subsamplings", new int[] {50, 300}, change.getSubsamplings());       // s = scaleSource / scaleTarget
    /*
     * Above (50, 300) subsampling shall be applied and the `gridToCRS` transform adjusted consequently.
     */
    final GridGeometry tg = change.build();
    extent = tg.getExtent();
    GridExtentTest.assertExtentEquals(extent, 0,  40, 110);
    GridExtentTest.assertExtentEquals(extent, 1,  -3,  27);               // NEAREST grid rounding mode.
    assertMatrixEquals("gridToCRS", new Matrix3(
            100,    0, 200,
              0, -300, 600,
              0,    0,   1), MathTransforms.getMatrix(tg.getGridToCRS(PixelInCell.CELL_CORNER)), STRICT);
    /*
     * The envelope is the intersection of the envelopes of `source` and `target` grid geometries, documented above.
     * That intersection should be approximately the same or smaller. Note that without the clipping documented in
     * `GridExtent(GridExtent, int...)` constructor, the envelope would have been larger.
     */
    GeneralEnvelope expected = new GeneralEnvelope(2);
    expected.setRange(0,  4200, 11300);
    expected.setRange(1, -7501,  1500);
    assertEnvelopeEquals(expected, tg.getEnvelope(), STRICT);
}
 

/**
 * Tests construction from a <cite>grid to CRS</cite> having a 0.5 pixel translation.
 * This translation happens in transform mapping <cite>pixel center</cite> when the
 * corresponding <cite>pixel corner</cite> transformation is identity.
 */
@Test
public void testShifted() {
    final long[]        low      = new long[] {100, 300};
    final long[]        high     = new long[] {200, 400};
    final GridExtent    extent   = new GridExtent(null, low, high, true);
    final MathTransform identity = MathTransforms.linear(new Matrix3(
            1, 0, 0.5,
            0, 1, 0.5,
            0, 0, 1));
    final GridGeometry grid = new GridGeometry(extent, PixelInCell.CELL_CENTER, identity, null);
    assertTrue("gridToCRS.isIdentity", grid.getGridToCRS(PixelInCell.CELL_CORNER).isIdentity());
}
 

/**
 * Tests the construction from a geospatial envelope and an extent.
 * The "grid to CRS" transform is inferred.
 *
 * @see GridExtentTest#testCornerToCRS()
 */
@Test
public void testFromExtentAndEnvelope() {
    final GeneralEnvelope aoi = new GeneralEnvelope(HardCodedCRS.WGS84);
    aoi.setRange(0,  40, 55);
    aoi.setRange(1, -10, 70);
    final GridGeometry grid = new GridGeometry(new GridExtent(null, new long[] {-20, -25}, new long[] {10, 15}, false), aoi);
    final Matrix matrix = MathTransforms.getMatrix(grid.getGridToCRS(PixelInCell.CELL_CORNER));
    assertMatrixEquals("cornerToCRS", new Matrix3(
            0.5,  0,   50,
            0,    2,   40,
            0,    0,    1), matrix, STRICT);
}
 

/**
 * Tests {@link GridGeometry#reduce(int...)}.
 */
@Test
public void testReduce() {
    final GridGeometry grid = new GridGeometry(
            new GridExtent(null, new long[] {336, 20, 4}, new long[] {401, 419, 10}, true),
            PixelInCell.CELL_CORNER, MathTransforms.linear(new Matrix4(
                    0,   0.5, 0,  -90,
                    0.5, 0,   0, -180,
                    0,   0,   2,    3,
                    0,   0,   0,    1)), HardCodedCRS.GEOID_3D);
    /*
     * Tests on the two first dimensions.
     */
    GridGeometry reduced = grid.reduce(0, 1);
    assertNotSame(grid, reduced);
    assertExtentEquals(new long[] {336, 20}, new long[] {401, 419}, reduced.getExtent());
    assertSame("CRS", HardCodedCRS.WGS84, reduced.getCoordinateReferenceSystem());
    assertArrayEquals("resolution", new double[] {0.5, 0.5}, reduced.getResolution(false), STRICT);
    assertMatrixEquals("gridToCRS", new Matrix3(
              0, 0.5,  -90,
              0.5, 0, -180,
              0,   0,    1), MathTransforms.getMatrix(reduced.getGridToCRS(PixelInCell.CELL_CORNER)), STRICT);
    /*
     * Tests on the last dimension.
     */
    reduced = grid.reduce(2);
    assertNotSame(grid, reduced);
    assertExtentEquals(new long[] {4}, new long[] {10}, reduced.getExtent());
    assertSame("CRS", HardCodedCRS.GRAVITY_RELATED_HEIGHT, reduced.getCoordinateReferenceSystem());
    assertArrayEquals("resolution", new double[] {2}, reduced.getResolution(false), STRICT);
    assertMatrixEquals("gridToCRS", new Matrix2(
              2, 3,
              0, 1), MathTransforms.getMatrix(reduced.getGridToCRS(PixelInCell.CELL_CORNER)), STRICT);
}
 

/**
 * Verifies consistency of cached transforms when using translations inferred from {@link PixelOrientation}.
 */
@Test
public void testCache2D() {
    assertSame(centerToCorner(2), centerToCorner2D());
    assertSame(PixelTranslation.translate(MathTransforms.identity(2), PixelInCell.CELL_CORNER, PixelInCell.CELL_CENTER),
               PixelTranslation.translate(MathTransforms.identity(2), PixelOrientation.UPPER_LEFT, PixelOrientation.CENTER, 0, 1));
}
 

@Override
public MathTransform getOriginalGridToWorld(
    final String coverageName,
    final PixelInCell pixInCell) {
  // just reuse super class implementation but ensure that we do not use a
  // cached raster2model
  synchronized (this) {
    raster2Model = null;
    return super.getOriginalGridToWorld(coverageName, pixInCell);
  }
}
 

/**
 * Converts a math transform from a "pixel in cell" convention to another "pixel in cell" convention.
 * This method concatenates −½, 0 or +½ translations on <em>all</em> dimensions before the given transform.
 * If the two given conventions are the same, then this method returns the given transform unchanged.
 *
 * <div class="note"><b>Example:</b>
 * if a given {@code gridToCRS} transform was mapping the <em>cell corner</em> to "real world" coordinates, then a call to
 * <code>translate(gridToCRS, {@link PixelInCell#CELL_CORNER CELL_CORNER}, {@link PixelInCell#CELL_CENTER CELL_CENTER})</code>
 * will return a new transform performing the following steps: first convert grid coordinates from <var>cell center</var>
 * convention ({@code desired}) to <var>cell corner</var> convention ({@code current}), then concatenate the given
 * {@code gridToCRS} transform which was designed for the <em>cell corner</em> convention.
 * The above-cited <var>cell center</var> → <var>cell corner</var> conversion is done by translating the grid coordinates
 * by +½, because the grid coordinates (0,0) relative to cell center is (½,½) relative to cell corner.</div>
 *
 * If the given {@code gridToCRS} is null, then this method ignores all other arguments and returns {@code null}.
 * Otherwise {@code current} and {@code desired} arguments must be non-null.
 *
 * @param  gridToCRS  a math transform from <cite>pixel</cite> coordinates to any CRS, or {@code null}.
 * @param  current    the pixel orientation of the given {@code gridToCRS} transform.
 * @param  desired    the pixel orientation of the desired transform.
 * @return the translation from {@code current} to {@code desired}, or {@code null} if {@code gridToCRS} was null.
 * @throws IllegalArgumentException if {@code current} or {@code desired} is not a known code list value.
 */
public static MathTransform translate(final MathTransform gridToCRS, final PixelInCell current, final PixelInCell desired) {
    if (gridToCRS == null || desired.equals(current)) {
        return gridToCRS;
    }
    final int dimension = gridToCRS.getSourceDimensions();
    final double offset = getPixelTranslation(desired) - getPixelTranslation(current);
    final int ci;               // Cache index.
    if (offset == -0.5) {
        ci = 2*dimension - 2;
    } else if (offset == 0.5) {
        ci = 2*dimension - 1;
    } else {
        ci = -1;
    }
    MathTransform mt;
    if (ci < 0 || ci >= translations.length) {
        mt = MathTransforms.uniformTranslation(dimension, offset);
    } else synchronized (translations) {
        mt = translations[ci];
        if (mt == null) {
            mt = MathTransforms.uniformTranslation(dimension, offset);
            translations[ci] = mt;
        }
    }
    return MathTransforms.concatenate(mt, gridToCRS);
}
 

/**
 * Creates a new grid geometry for the given extent.
 * This method should be invoked only when no existing {@link GridGeometry} can be used as template.
 */
GridGeometry createGridCRS(final Variable variable) {
    final List<Dimension> dimensions = variable.getGridDimensions();
    final long[] upper = new long[dimensions.size()];
    for (int i=0; i<upper.length; i++) {
        final int d = (upper.length - 1) - i;           // Convert CRS dimension to netCDF dimension.
        upper[i] = dimensions.get(d).length();
    }
    return new GridGeometry(new GridExtent(null, null, upper, false), PixelInCell.CELL_CENTER, gridToCRS, crs);
}
 

/**
 * Invoked by JAXB only at unmarshalling time.
 *
 * @see #getPixelInCell()
 */
private void setPixelInCell(final PixelInCell value) {
    if (pixelInCell == null) {
        pixelInCell = value;
    } else {
        MetadataUtilities.propertyAlreadySet(DefaultImageDatum.class, "setPixelInCell", "pixelInCell");
    }
}
 

/**
 * Returns the concatenation of all transformation steps from the given source to the given target.
 * The transform maps grid coordinates (not envelopes).
 *
 * @param  source     the source grid geometry.
 * @param  crsChange  the change of coordinate reference system, or {@code null} if none.
 * @param  target     the target grid geometry.
 * @param  anchor     whether we want the transform for cell corner or cell center.
 */
private static MathTransform path(final GridGeometry source, final CoordinateOperation crsChange,
        final GridGeometry target, final PixelInCell anchor) throws NoninvertibleTransformException
{
    MathTransform step1 = source.getGridToCRS(anchor);
    MathTransform step2 = target.getGridToCRS(anchor);
    if (crsChange != null) {
        step1 = MathTransforms.concatenate(step1, crsChange.getMathTransform());
    }
    if (step1.equals(step2)) {                                          // Optimization for a common case.
        return MathTransforms.identity(step1.getSourceDimensions());
    } else {
        return MathTransforms.concatenate(step1, step2.inverse());
    }
}
 

/**
 * Creates an image CRS using a two-dimensional affine or Cartesian coordinate system.
 *
 * @param cartesian {@code true} for a Cartesian coordinate system, or {@code false} for an affine one.
 */
private static DefaultImageCRS create(final boolean cartesian) {
    return new DefaultImageCRS(Collections.singletonMap(DefaultImageCRS.NAME_KEY, "An image CRS"),
            new DefaultImageDatum(Collections.singletonMap(DefaultImageDatum.NAME_KEY, "C1"), PixelInCell.CELL_CENTER),
            cartesian ? HardCodedCS.GRID : new DefaultAffineCS(
                    Collections.singletonMap(DefaultAffineCS.NAME_KEY, "Grid"),
                            HardCodedAxes.COLUMN, HardCodedAxes.ROW));
}
 

/**
 * Creates a grid geometry with the given extent and scale for testing purpose.
 * An arbitrary translation of (2,3) is added to the "grid to CRS" conversion.
 */
private static GridGeometry grid(int xmin, int ymin, int xmax, int ymax, int xScale, int yScale) throws TransformException {
    GridExtent extent = new GridExtent(null, new long[] {xmin, ymin}, new long[] {xmax, ymax}, true);
    Matrix3 gridToCRS = new Matrix3();
    gridToCRS.m00 = xScale;
    gridToCRS.m11 = yScale;
    gridToCRS.m02 = 200;            // Arbitrary translation.
    gridToCRS.m12 = 500;
    return new GridGeometry(extent, PixelInCell.CELL_CORNER, MathTransforms.linear(gridToCRS), null);
}
 

/**
 * Creates a new grid geometry from a grid extent and a mapping from cell coordinates to "real world" coordinates.
 * At least one of {@code extent}, {@code gridToCRS} or {@code crs} arguments shall be non-null.
 * If {@code gridToCRS} is non-null, then {@code anchor} shall be non-null too with one of the following values:
 *
 * <ul>
 *   <li>{@link PixelInCell#CELL_CENTER} if conversions of cell indices by {@code gridToCRS} give "real world"
 *       coordinates close to the center of each cell.</li>
 *   <li>{@link PixelInCell#CELL_CORNER} if conversions of cell indices by {@code gridToCRS} give "real world"
 *       coordinates at the corner of each cell. The cell corner is the one for which all grid indices have the
 *       smallest values (closest to negative infinity).</li>
 * </ul>
 *
 * <div class="note"><b>API note:</b>
 * there is no default value for {@code anchor} because experience shows that images shifted by ½ pixel
 * (with pixels that may be tens of kilometres large) is a recurrent problem. We want to encourage developers
 * to always think about wether their <cite>grid to CRS</cite> transform is mapping pixel corner or center.</div>
 *
 * <div class="warning"><b>Upcoming API generalization:</b>
 * the {@code extent} type of this method may be changed to {@code GridEnvelope} interface in a future Apache SIS version.
 * This is pending <a href="https://github.com/opengeospatial/geoapi/issues/36">GeoAPI update</a>.
 * In addition, the {@code PixelInCell} code list currently defined in the {@code org.opengis.referencing.datum} package
 * may move in another package in a future GeoAPI version because this type is no longer defined by the ISO 19111 standard
 * after the 2018 revision.</div>
 *
 * @param  extent     the valid extent of grid coordinates, or {@code null} if unknown.
 * @param  anchor     {@linkplain PixelInCell#CELL_CENTER Cell center} for OGC conventions or
 *                    {@linkplain PixelInCell#CELL_CORNER cell corner} for Java2D/JAI conventions.
 * @param  gridToCRS  the mapping from grid coordinates to "real world" coordinates, or {@code null} if unknown.
 * @param  crs        the coordinate reference system of the "real world" coordinates, or {@code null} if unknown.
 * @throws NullPointerException if {@code extent}, {@code gridToCRS} and {@code crs} arguments are all null.
 * @throws MismatchedDimensionException if the math transform and the CRS do not have consistent dimensions.
 * @throws IllegalGridGeometryException if the math transform can not compute the geospatial envelope or resolution from the grid extent.
 */
public GridGeometry(final GridExtent extent, final PixelInCell anchor, final MathTransform gridToCRS, final CoordinateReferenceSystem crs) {
    if (gridToCRS != null) {
        ensureDimensionMatches(gridToCRS.getSourceDimensions(), extent);
        ArgumentChecks.ensureDimensionMatches("crs", gridToCRS.getTargetDimensions(), crs);
    } else if (crs == null) {
        ArgumentChecks.ensureNonNull("extent", extent);
    }
    try {
        this.extent      = extent;
        this.gridToCRS   = PixelTranslation.translate(gridToCRS, anchor, PixelInCell.CELL_CENTER);
        this.cornerToCRS = PixelTranslation.translate(gridToCRS, anchor, PixelInCell.CELL_CORNER);
        this.envelope    = computeEnvelope(gridToCRS, crs, null);   // 'gridToCRS' specified by the user, not 'this.gridToCRS'.
        this.resolution  = resolution(gridToCRS, extent);           // 'gridToCRS' or 'cornerToCRS' does not matter here.
        this.nonLinears  = findNonLinearTargets(gridToCRS);
    } catch (TransformException e) {
        throw new IllegalGridGeometryException(e, "gridToCRS");
    }
}
 

/**
 * Returns a transform from center to corner with the given number of dimensions.
 */
private static MathTransform centerToCorner(final int dimension) {
    return PixelTranslation.translate(MathTransforms.identity(dimension), PixelInCell.CELL_CENTER, PixelInCell.CELL_CORNER);
}
 

/**
 * Creates a math transform using the information provided.
 *
 * @return The math transform.
 * @throws IllegalStateException if the grid range or the envelope were not set.
 */
public static MathTransform createTransform(
    final double[] idRangePerDimension,
    final MultiDimensionalNumericData fullBounds) throws IllegalStateException {
  final GridToEnvelopeMapper mapper = new GridToEnvelopeMapper();
  final boolean swapXY = mapper.getSwapXY();
  final boolean[] reverse = mapper.getReverseAxis();
  final PixelInCell gridType = PixelInCell.CELL_CORNER;
  final int dimension = 2;
  /*
   * Setup the multi-dimensional affine transform for use with OpenGIS. According OpenGIS
   * specification, transforms must map pixel center. This is done by adding 0.5 to grid
   * coordinates.
   */
  final double translate;
  if (PixelInCell.CELL_CENTER.equals(gridType)) {
    translate = 0.5;
  } else if (PixelInCell.CELL_CORNER.equals(gridType)) {
    translate = 0.0;
  } else {
    throw new IllegalStateException(
        Errors.format(ErrorKeys.ILLEGAL_ARGUMENT_$2, "gridType", gridType));
  }
  final Matrix matrix = MatrixFactory.create(dimension + 1);
  final double[] minValuesPerDimension = fullBounds.getMinValuesPerDimension();
  final double[] maxValuesPerDimension = fullBounds.getMaxValuesPerDimension();
  for (int i = 0; i < dimension; i++) {
    // NOTE: i is a dimension in the 'gridRange' space (source
    // coordinates).
    // j is a dimension in the 'userRange' space (target coordinates).
    int j = i;
    if (swapXY) {
      j = 1 - j;
    }
    double scale = idRangePerDimension[j];
    double offset;
    if ((reverse == null) || (j >= reverse.length) || !reverse[j]) {
      offset = minValuesPerDimension[j];
    } else {
      scale = -scale;
      offset = maxValuesPerDimension[j];
    }
    offset -= scale * (-translate);
    matrix.setElement(j, j, 0.0);
    matrix.setElement(j, i, scale);
    matrix.setElement(j, dimension, offset);
  }
  return ProjectiveTransform.create(matrix);
}