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

[Jonathan Blower]

Hi Martin,

Yes, I can see how this works - and clearly a replacement for GridGeometry is needed in the Coverage world.  My main question is - what happens if the transformation from real-world to grid coordinates is not affine?  We have quite a few situations like this, and they are becoming more and more common in global and coastal oceanography.  In these cases the grid is complex and distorted, perhaps to avoid the polar singularity and/or to distort around the coasts.

My second question: in the method

>     wmsImage.drawRenderedImage(coverageImage, gridToCRS);

what kind of interpolation is used from coverageImage to the wmsImage?  Is it nearest-neighbour?  So does the algorithm look a bit like this:

foreach pixel in wmsImage:
    find the nearest i,j grid indices in coverageImage
    extract the data at coverageImage(i,j)
    populate the pixel with that data value

or does it extract the whole of coverageImage in one operation and then do the subsetting in memory?

Thanks,
Jon

------------------------------

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


------------------------------

Message: 3
Date: Mon, 19 Apr 2010 16:39:34 +0200
From: Martin Desruisseaux <martin.desruisseaux at geomatys.fr>
Subject: Re: [geotk] Mapping coordinate system axis order to i,j
	indices on	a map images
To: geotoolkit at lists.osgeo.org
Message-ID: <4BCC6B26.6050409 at geomatys.fr>
Content-Type: text/plain; charset=UTF-8; format=flowed

I forget to said that one issue is that ISO 19123 does not have this 
"GridGeometry" object which existed in OGC 01-004. "GridGeometry" and 
"SampleDimensions" are the two legacy OGC 01-004 objects that I'm tempted to 
port to the new ISO 19123 framework. This is similar with what we did for the 
referencing module, where MathTransform has been ported from the legacy OGC 
01-009 to ISO 19111. Please let me known what you think.

	Martin