Java Code Examples for org.opengis.referencing.datum.PixelInCell#CELL_CORNER

/**
 * 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());
}
 

/**
 * 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 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);
}
 

/**
 * 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());
}
 

/**
 * 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 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 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);
}
 

/**
 * 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);
}
 

/**
 * 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 {@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);
}
 

/**
 * 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);
}
 

/**
 * Tests {@link GridDerivation#subgrid(Envelope, double...)}.
 * Contrarily to {@link #testSubExtent()}, this method checks the full {@link GridGeometry}.
 *
 * @throws TransformException if an error occurred during computation.
 */
@Test
@DependsOnMethod("testSubExtent")
public void testSubgridFromEnvelope() throws TransformException {
    final GeneralEnvelope envelope = new GeneralEnvelope(HardCodedCRS.WGS84_φλ);
    envelope.setRange(0, -70, +80);
    envelope.setRange(1,   5,  15);
    final MathTransform gridToCRS = MathTransforms.linear(new Matrix3(
            0,   0.5, -90,
            0.5, 0,  -180,
            0,   0,     1));
    GridGeometry grid = new GridGeometry(PixelInCell.CELL_CORNER, gridToCRS, envelope, GridRoundingMode.NEAREST);
    assertExtentEquals(new long[] {370, 40}, new long[] {389, 339}, grid.getExtent());
    assertEnvelopeEquals(envelope, grid.getEnvelope(), STRICT);
    /*
     * Set a sub-region. The grid extent and "grid to CRS" transform shall be adjusted
     * in such a way that envelope computed from the new grid geometry is the same.
     */
    envelope.setRange(0, -50, +30);
    envelope.setRange(1,   8,  12);
    grid = grid.derive().subgrid(envelope, 1, 2).build();
    assertExtentEquals(new long[] {94, 40}, new long[] {95, 119}, grid.getExtent());
    assertEnvelopeEquals(envelope, grid.getEnvelope(), STRICT);
    assertMatrixEquals("gridToCRS", new Matrix3(
            0, 1,  -90,
            2, 0, -180,
            0, 0,    1), MathTransforms.getMatrix(grid.getGridToCRS(PixelInCell.CELL_CORNER)), STRICT);
    /*
     * A sub-region again but with a requested resolution which is not a divisor of the actual resolution.
     * It will force GridGeometry to adjust the translation term to compensate. We verify that the adustment
     * is correct by verifying that we still get the same envelope.
     */
    grid = grid.derive().subgrid(envelope, 3, 2).build();
    assertExtentEquals(new long[] {94, 13}, new long[] {95, 39}, grid.getExtent());
    assertEnvelopeEquals(envelope, grid.getEnvelope(), STRICT);
    MathTransform cornerToCRS = grid.getGridToCRS(PixelInCell.CELL_CORNER);
    assertMatrixEquals("gridToCRS", new Matrix3(
            0, 3,  -89,
            2, 0, -180,
            0, 0,    1), MathTransforms.getMatrix(cornerToCRS), STRICT);

    DirectPosition2D src = new DirectPosition2D();
    DirectPosition2D tgt = new DirectPosition2D();
    DirectPosition2D exp = new DirectPosition2D();
    src.x = 94; src.y = 13; exp.x = -50; exp.y =  8; assertEquals("Lower corner", exp, cornerToCRS.transform(src, tgt));
    src.x = 96; src.y = 40; exp.x = +31; exp.y = 12; assertEquals("Upper corner", exp, cornerToCRS.transform(src, tgt));
}
 

/**
 * Ensures that slicing on a corner point does not fail, but gives back a grid geometry centered on a pixel corner.
 *
 * @throws TransformException if we cannot build our test point.
 */
@Test
public void testSliceOnCorner() throws TransformException {
    GridGeometry base = grid(-132, 327, 986, 597, 2, 3);
    /*
     * First of all, we'll try by focusing on the last pixel.
     */
    final DirectPosition2D gridUpperCorner = new DirectPosition2D(
            base.extent.getHigh(0),
            base.extent.getLow(1));

    final DirectPosition geoUpperCorner = base.getGridToCRS(PixelInCell.CELL_CENTER)
            .transform(gridUpperCorner, null);

    GridGeometry slice = base.derive()
            .slice(geoUpperCorner)
            .build();

    long[] expectedGridPoint = {(long) gridUpperCorner.x, (long) gridUpperCorner.y};
    assertExtentEquals(expectedGridPoint, expectedGridPoint, slice.extent);
    /*
     * We will now try to focus on a point near the envelope edge. Note that slicing ensures to return a valid grid
     * for any point INSIDE the envelope, it's non-determinist about points perfectly aligned on the edge.
     * So, here we will test a point very near to the envelope edge, but still into it.
     */
    final GeneralEnvelope grid3d = new GeneralEnvelope(3);
    grid3d.setEnvelope(0, 0, 0, 1920, 1080, 4);

    final DefaultCompoundCRS crs3d = new DefaultCompoundCRS(
            Collections.singletonMap("name", "geo3d"),
            HardCodedCRS.WGS84,
            HardCodedCRS.TIME);

    final GeneralEnvelope geo3d = new GeneralEnvelope(crs3d);
    geo3d.setEnvelope(-180, -90, 1865.128, 180, 90, 1865.256);
    base = new GridGeometry(
            PixelInCell.CELL_CORNER,
            MathTransforms.linear(Matrices.createTransform(grid3d, geo3d)),
            geo3d,
            GridRoundingMode.NEAREST);

    final GeneralDirectPosition geo3dUpperCorner = new GeneralDirectPosition(geo3d.getUpperCorner());
    IntStream.range(0, geo3dUpperCorner.getDimension())
            .forEach(idx -> geo3dUpperCorner.coordinates[idx] -= 1e-7);

    slice = base.derive()
            .slice(geo3dUpperCorner)
            .build();

    // Build expected grid point focused after slicing. We expect it to be upper corner.
    expectedGridPoint = DoubleStream.of(grid3d.getUpperCorner().getCoordinate())
            .mapToLong(value -> (long) value)
            .map(exclusiveValue -> exclusiveValue - 1)        // Exclusive to inclusive.
            .toArray();

    assertExtentEquals(expectedGridPoint, expectedGridPoint, slice.extent);
}