[geotk] Mapping coordinate system axis order to i,j indices on
a map images
Martin Desruisseaux
martin.desruisseaux at geomatys.fr
Mon Apr 19 10:36:46 EDT 2010
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
More information about the Geotoolkit
mailing list