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

[Martin Desruisseaux]

Hellp John

Le 20/04/10 11:53, Jonathan Blower a écrit :
>> Actually the most robust approach is to ensure that we never need to known if
>> axes are swapped/flipped or not, especially since non-affine transforms provide
>> no garantee that axes swapping/flipping do not vary with the position.
>
> Agreed.  In ncWMS we have a HorizontalGrid object that specifies the mapping from grid to WGS84 lat/lon coordinates, rather like the ReferenceableGrid in ISO19123.  It's up to subclasses to specify how the mapping is done (in some very difficult cases it's essentially a look-up table!).  Geotoolkit then handles the transformation from lat/lon to the image CRS.  So we do the transformation in two steps, with lat/lon as a "common language" bridging the gap between Geotoolkit and general NetCDF coordinate systems.
>
> In this system, efficiencies can be made (i.e. fewer transformations need to be performed) if grid axes are aligned with the axes of external CRSs, which is what has motivated my original question on this thread.

I understand. It is fine to check for special cases as an optimization (e.g. "if 
the transform is affine, then check axis order and use that optimized rendering 
chain") when we also have, as a fallback, a general mechanism capable to work 
without explicit knowledge of axis order.

One think that may help to gain efficienty would be to make HorizontalGrid 
implements MathTransform (if possible), so that the referencing module can 
handle it like any other MathTransform. In addition of that, the code could also 
check if the HorizontalGrid is actually affine (basically checking that the 
(dx,dy) translation between a pixel and the neighboard pixel is constant for all 
pixels). For example:

    MathTransform gridToCRS;
    if (theGridIsAffine) {
         Matrix m = ... put the (dx,dy) and the translations here.
         gridToCRS = ProjectiveTransform.create(m);
    } else {
         gridToCRS = theHorizontalGrid;
    }

Most parts of the referencing module can work with arbitrary MathTransforms (so 
the HorizontalGrid implementation should work well), but contain optimizations 
for the affine transform cases. Concatenating the above gridToCRS with the 
"lat/long to image CRS" (using for instance ConcatenatedTransform.create(step1, 
step2)) can lead to maximal performance if the referencing module detect that 
all the transformation chain is made of affine transforms.

Optionally, you could make things more formal by defining the CRS of the grid:

    CoordinateReferenceSystem coverageCRS = DefaultGeographicCRS.WGS84; // or 
anything else
    CoordinateReferenceSystem gridCRS; // To be defined below.
    MathTransform gridToCRS = ... // as explained above.

    String crsName = "I suggest to use the filename or some file ID";
    gridCRS = new DefaultDerivedCRS(crsName, coverageCRS, gridToCRS, 
DefaultCartesienCS.GRID);

Then given the "wmsCRS":

    CoordinateReferenceSystem wmsCRS = CRS.decode("EPSG:4326");

(note that the above will have (lat,long) axis order if you are using the JavaDB 
or HSQL embedded EPSG database. You can also force (long,lat) axis order by 
adding the boolean 'true' argument), then the following will take care of the 
full transformation chain from the gridCRS to the wmsCRS, including any axes 
swpapping or flipping which may be needed:

    MathTransform gridToWMS = CRS.findMathTransform(gridCRS, wmsCRS);

In Geotk, we go one step further by going straight from gridCRS to displayCRS, 
but the idea is the same. In some cases the whole chain is optimized as a single 
AffineTransform object.


> However, I'd be really interested to compare our method (in ncWMS) with Geotoolkit's NetCDF-reading code, particularly in terms of performance.  One worry I have about the Geotoolkit method (which is very neat by the way) is that the extraction of data from the NetCDF coverage is handled "under the hood" in Java2D.  Is it possible that the implementation will read the entirety of a high-resolution field in order to extract just a few grid cells for a low-resolution image?

It is possible, depending on the ImageReadParam. The javax.imageio.ImageReader 
API does not allow reading a given set of points. The only thing that we can 
specify are:

   1) The rectangle bounds of the pixels to read
   2) The subsampling (read 1 pixel, skip 2, etc.), which must
      be constant over all the Rectangle area we are reading.

In term of NetCDF API, this is equivalent to performing a call to 
Variable.read(Section) where Section has been created using the following 
constructor:

   Section(int[] origin, int[] size, int[] stride)

This is not as flexible as the technic documented in ncWMS DataReadingStrategy 
(nicely documented by the way :) The illustrations are really convenient).

Note that further performance improvement are planed for the Geotk 
NetcdfImageReader. Current implementation stores the values in the image using a 
JAI WritableRectIter. In a future version, I would like to write the values 
straight into the image DataBuffer.


> Also, is the Geotk method essentially limited to affine transformations or would there be a way of handling arbitrary transformations?  Looking at the Java2D API, it seems that handling non-affined transformations might be difficult, but perhaps I've missed something.

Java2D is restricted to AffineTransform. However Geotk can handle arbitrary 
transformations, including HorizontalGrid if it implements MathTransform. If the 
transform is not affine, you can use the following code (if you have a 
GridCoverage2D object; an alternative approach proposed after):

    RenderedImage rasterData = ...;
    MathTransform gridToCRS = ...;
    CoordinateReferenceSystem coverageCRS = ...;
    String coverageName = "I suggest they layer name, maybe with the date";
    GridCoverage2D coverage = gridCoverageFactory.create(coverageName, 
rasterData, coverageCRS, gridToCRS, null, null, null);

    // Resample now
    coverage = (GridCoverage2D) Operations.DEFAULT.resample(coverage, wmsCRS);

    // If you want more control (for example on the resolution of the output
    // image), you can also use the flavor expecting a GridGeometry argument.

More on the "Resample" operation here:

http://www.geotoolkit.org/apidocs/org/geotoolkit/coverage/processing/operation/Resample.html


As an alternative, you can also make the low-level processing yourself. The key 
point is to create an instance of javax.media.jai.Warp object that represent the 
transform, then use the JAI "Warp" operation. You can wraps your HorizontalGrid 
into a JAI Warp object with the following method:

http://www.geotoolkit.org/apidocs/org/geotoolkit/referencing/operation/transform/WarpTransform2D.html#getWarp%28java.lang.CharSequence,%20org.opengis.referencing.operation.MathTransform2D%29

Then give that Warp to the following JAI method:

http://download.java.net/media/jai/javadoc/1.1.3/jai-apidocs/javax/media/jai/operator/WarpDescriptor.html

Current Geotk implementation of the Warp object returned by WarpTransform2D is 
not very fast (while not extremely slow neither). In a future version, I hope to 
make it faster by leveraging the JAI WarpGrid object:

http://download.java.net/media/jai/javadoc/1.1.3/jai-apidocs/javax/media/jai/WarpGrid.html


> We have found that it's most efficient to strike a balance between reading the whole field in a single operation and reading individual pixels in individual operations: see http://www.resc.rdg.ac.uk/trac/ncWMS/browser/trunk/src/java/uk/ac/rdg/resc/ncwms/cdm/DataReadingStrategy.java.

Yes, I have read the nice DataReadingStrategy javadoc. Sound like a nice 
strategy. It doesn't fit in Java2D architecture, but we don't have to use that 
architecture for everything. It is just one possible implementation among others.


> Where should I go to grab your NetCDF-reading code, and do you have a simple example I could work from?  I've had a quick look on the Geotk website but couldn't find much on NetCDF yet.

We are lacking documentation (:

1) Ensure that the geotk-coverageio-netcdf module and all its
    dependencies are on the classpath. If you are using the geotk
    bundles, pickup the geotk-bundle-coverage.pack.gz bundle.

2) Try the following code:


// Create the reader and set the input. Note that we don't
// specify the format (NetCDF) anywhere - the Java Image I/O
// framework will autodect. The input could be a PNG, JPEG,
// TIFF or other formats too.
GridCoverageReader reader = new ImageCoverageReader();
reader.setInput(new File("myFile.netcdf"));

// Set the geodetic envelope to the region of interest.
// The envelope can be in any CRS - Geotk will perform
// the needed transformation itself.
GridCoverageReadParam param = new GridCoverageReadParam();
param.setEnvelope(...);

// The resolution shall be in the same units than the
// envelope. It doesn't need to be the unit of the grid data.
param.setResolution(...);

// Get the list of variable in the NetCDF file. The let
// assume that we have a variable named "Temperature".
List<String> variables = reader.getCoverageNames();
int imageIndex = variables.indexOf("Temperature");

// Get the coverage
GridCoverage2D coverage = reader.read(imageIndex, param);


The above code work at a high level, in geodetic coordinate space. If you prefer 
to work at a lower level - in pixel space - you can use NetcdfImageReader 
directly. Note that the way to use it is very similar to the above, except that 
the coordinates are in pixels rather than geodetic units:

NetcdfImageReader reader = new NetcdfImageReader();
reader.setInput(new File("myFile.netcdf"));
ImageReadParam param = reader.getDefaultReadParam();
param.setSourceRegion(... the envelope in pixel units ...);
param.setSourceSubsampling(... the resolution in pixel units ...);
RenderedImage image = reader.read(imageIndex, param);

The way to handle dimensions above 2 are a little bit convolved, because the 
Java2D is mostly a 2-dimensional API. However there is a mechanism for that in 
SpatialImageReadParam, which could be the subject of an other email.

	Regards,

		Martin