[geotk] Mapping coordinate system axis order to i,j indices on a map images

[Martin Desruisseaux]

Hello John

Le 19/04/10 15:01, Jonathan Blower a écrit :
> I'm not sure whether this is a Geotk question, or one for the WMS/ISO
> 19111 community.  Given a general 2D coordinate system, how do I know
> which axis should map to the horizontal dimension of a map image, and
> which axis should map to the vertical?  For example, EPSG:4326 uses
> lat/lon order, whereas CRS:84 uses lon/lat, but in both cases one would
> reasonably expect longitude to map to the horizontal direction.
> However, I don't know how a computer would make this inference.

The OGC 01-009 specification ("Grid Coverage implementation specification") had 
a rather clear answer. The answer in ISO 19123 is less obvious to me. So my 
reply above concern OGC 01-009; it would be nice if we can map this answer to 
ISO 19123 objects, and if ISO 19123 has no good equivalence decide which OGC 
01-009 objects (if any) we keep.

In OGC 01-009, "DiscreteGridPointCoverage" was called simply "GridCoverage". The 
legacy GridCoverage interface still around. If you can take a look at the 
interface below, I suggest that you ignore everything except only one method: 
getGridGeometry():

http://www.geoapi.org/snapshot/pending/org/opengis/coverage/grid/GridCoverage.html

The interresting part is that the GridGeometry provides the mapping from the 
grid coordinates to the CRS coordinates in the form of a MathTransform object. 
In the vast majority of cases, the MathTransform is actually an AffineTransform.

I strongly suggest that every GridCoverage or DiscreteGridPointCoverage 
implementations manage their "grid to CRS" mapping with AffineTransforms, never 
with Envelopes and conversion formulas coded manually. In the 2D case, the 
standard java.awt.geom.AffineTransform class does a very good job. In the n-D 
case, you can use 
org.geotoolkit.referencing.operation.transform.ProjectiveTransform.

If the axes are not interchanged and the y axis is not upside-down, then the 
affine matrix looks like:

       ┌                    ┐
       │ xRes      0   xMin │
       │    0   yRes   yMin │
       │    0      0      1 │
       └                    ┘

where "xRes" and "yRes" are the pixel resolution along x and y axis, and "xMin", 
"yMin" and the minimal geodetic coordinates.

If the y-axis is upside-down, then the matrix is:

       ┌                    ┐
       │ xRes      0   xMin │
       │    0  -yRes   yMax │
       │    0      0      1 │
       └                    ┘

If the x and y axes are interchanged (as in EPSG:4326), then the matrix is;

       ┌                    ┐
       │    0   xRes   xMin │
       │ yRes      0   yMin │
       │    0      0      1 │
       └                    ┘

As you can see, by using affine transform instead than envelope, there is no 
ambiguity: every axes swapping/flipping cases are clearly described. You can 
easily calculate the geodetic envelope from the GridEnvelope and the gridToCRS 
affine transnform, so the user don't need to specify the envelope.

As a side note, the following may be convenient:

http://www.geotoolkit.org/apidocs/org/geotoolkit/coverage/grid/GeneralGridGeometry.html

You can use the constructor expecting GridEnvelope + Envelope, and it will 
calculate the affine transform for you with some assumptions about axes. Or you 
can use the constructor expecting GridEnvelope + MathTransform and it will 
calculate the envelope for you.

So when we have our GridGeometry object, we can know if axes are swapped / 
flipped or not by inspection of the matrix coefficients. Actually in most cases 
you don't need to inspect the coefficients at all. If you are using Java2D for 
producing your WMS image, you can do:

     RenderedImage coverageImage = ...;
     AffineTransform gridToCRS = ...;
     BufferedImage wmsImage = ...;
     Graphics2D gr = wmsImage.createGraphics();
     gr.transform(... put your 'CRS to target pixels' transform here ...);
     wmsImage.drawRenderedImage(coverageImage, gridToCRS);
     gr.dispose();


Constructing the "grid to CRS" affine transform
-----------------------------------------------
For constructing the affine transform, you can perform the following steps with 
Geotk.

     CoordinateSystem sourceCS = DefaultCartesianCS.GENERIC_2D;
     CoordinateSystem targetCS = coverageCRS.getCoordinateSystem();
     Matrix m = AbstractCS.swapAndScaleAxis(sourceCS, targetCS);

The key method call is "swapAndScaleAxis", which is documented here:

http://www.geotoolkit.org/apidocs/org/geotoolkit/referencing/cs/AbstractCS.html#swapAndScaleAxis%28org.opengis.referencing.cs.CoordinateSystem,%20org.opengis.referencing.cs.CoordinateSystem%29

At this point, you can inspect the matrix yourself, or let Geotk continue the job:

     MathTransform gridToCRS = ProjectiveTransform.create(m);

http://www.geotoolkit.org/apidocs/org/geotoolkit/referencing/operation/transform/ProjectiveTransform.html

In Geotk implementation, the created 'gridToCRS' will actually be an instance of 
java.awt.geom.AffineTransform if the sourceCS and targetCS are two-dimensional.

Finally (assuming 2D case, but you can also use the GeneralXXX flavors):

     GridEnvelope2D ge = new GridEnvelope2D(x, y, width, height);
     GridGeometry gg = new GeneralGridGeometry(ge, PixelOrientation.UPPER_LEFT, 
gridToCRS, coverageCRS, null);

The above asumed that the "gridToCRS" transform map the upper-left corner of 
pixels, which is the case most of the time.

You can check the geodetic envelope with:

     System.out.println(gg.getEnvelope());


Hope it help? Please let me know if the above doesn't suite your need.

	Martin