Kriging with pMode == 1 leaves a black column to the right

[rafaelalmeida]

While doing kriging with pMode == 1 (the mode in which you supply a GridGeometry instead of a collection of points to be interpolated), the output image has a black column to the right:

This is the code that generates the GridGeometry passed to inInterpolationGrid:

GridEnvelope mapGrid = new GridEnvelope2D(0, 0, areaSizeInPixels[0], areaSizeInPixels[1]);
kriging.inInterpolationGrid = new GridGeometry2D(mapGrid, createEnvelopeFromBoundingBox(bounds));

( Where createEnvelopeFromBoundingBox returns an Envelope2D )

The problem seems to arise from a combination of default GeoTools behaviour and this line:

https://github.com/moovida/jgrasstools/blob/56e8fe50de7647cb579a9b28363842353ec20e13/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/modules/statistics/kriging/OmsKriging.java#L568

According to the default Geotools behavior when the pixel anchor is not specified when creating a grid geometry (as is the case above), Geotools will use the OGC default, which is to consider the anchor to be in the middle of the pixel.

In this case, when the coordinate is minimal, Geotools will actually map its resulting x coordinate to -0.5, not 0.5, as I could verify with the code below:

GridEnvelope grid = new GridEnvelope2D(0, 0, 100, 100);
Envelope2D envelope = new Envelope2D(CRS.decode("EPSG:3857"), 10.0, 10.0, 1.0, 1.0);
GridGeometry2D geom = new GridGeometry2D(grid, envelope);

DirectPosition out = geom.getCRSToGrid2D().transform(
        (DirectPosition) new DirectPosition2D(10.0, 10.0), (DirectPosition) new DirectPosition2D());

System.out.println(out.toString());

Which outputs: DirectPosition2D[-0.5000000000001137, 99.5]

Since the line from the OmsKriging module I linked above just casts the coordinates into an int, it drops the decimal cases from -0.5, which becomes -0, which is equivalent to 0. In the second column, the output becomes 0.5, and dropping the case again we have 0.

I don't know why Geotools behaves this way, or if this behaviour is correct.

Should jgrasstools try to compensate for this behavior or let it up to the user to supply a GridGeometry that behaves the way jgrasstools expects it to do?

[moovida]

Hi @rafaelalmeida , this is quite interesting, I haven't seen this problem in a while. Indeed the pixel anchor can lead to issues, but when it does I thought it would always affect one row and one column.

This is a nice catch of yours and I sure would like to fix this.
Obviously it would be great if jgrasstools could compensate for this problem. Do you have a recipe for that? I think different readers have different way to handle that. For example I think asc and tiff might handle that differently, which made merging of these two types hard.

On the bright side, if there was a way to handle that, it would need to be changed only in the OmsRasterReader ...

Thanks for reporting this!

[rafaelalmeida]

Hi @moovida, I'm not sure where the OmsRasterReader would get into the picture.

In this case, I don't get the GridGeometry from a file, since I don't have a file in the first place, just an Envelope and some measurements on it. That is, I'm trying to build the file (in my case, a GeoTIFF) in the first place.

So I think it would have to be fixed directly in OmsKriging, and in other places that rely on GridGeometry to generate points to be processed.

To be honest, I don't know why Geotools behaves that way, or if it's correct. At a first glance, no matter if we are in pixel center or pixel corner mode, the grid is defined in the [0, 100) (0 inclusive, 100 exclusive) interval, and -0.5 is off the grid. But maybe there's a reason it makes sense to return this value. I reported it to the Geotools project as a possible bug, seeking clarification (https://osgeo-org.atlassian.net/browse/GEOT-5177).

I ran a few more tests and may have found a way to correct this in jgrasstools side only. It involves getting a modified transform from the GridGeometry2D object, with pixel orientation set to UPPER_LEFT, for example.

Even then there are still small rounding errors, so we can't just cast to int blindly. Casting to int just drops the decimal, so 0.999999999 would become 0 instead of 1. So I believe it's better to Math.round() before casting.

Here's some code to demonstrate that:

GridEnvelope grid = new GridEnvelope2D(0, 0, 10, 10);
Envelope2D envelope = new Envelope2D(CRS.decode("EPSG:3857"), 0, 0, 1.0, 1.0);
GridGeometry2D geometry = new GridGeometry2D(grid, envelope);

MathTransform2D transform = geometry.getCRSToGrid2D(PixelOrientation.UPPER_LEFT);

int errorCount = 0;
for (int x = 0; x < 10; x++) {
    for (int y = 0; y < 10; y++) {
        DirectPosition P = (DirectPosition) new DirectPosition2D(
                x * .1, 1.0 - y * .1);
        DirectPosition projection = transform.transform(P, null);

        double coordinates[] = projection.getCoordinate();
        int xp = (int) Math.round(coordinates[0]);
        int yp = (int) Math.round(coordinates[1]);

        System.out.println("Expected: x = " + x + ", y = " + y
                + ", got: x = " + xp + ", y = " + yp);

        if (xp != x || yp != y) {
            errorCount += 1;
        }
    }
}

System.out.println("\n" + errorCount + " error(s).");

Which outputs:

(...) (The projected values for each position, for visual conference)
0 error(s).

Should we decide to go this route, I should probably test more coordinate systems to see if it works.

If it does, are you familiar with other places in jgrasstools that use GridGeometry2D? We should probably check those too.