[geotk] Grib file coordinate transformation
Martin Desruisseaux
martin.desruisseaux at geomatys.fr
Mon Jan 30 14:51:55 EST 2012
Hello Pablo
Ideally, the NetCDF reader would have setup the ProjectedCRS automatically. In
practice, this work is in progress right now so a temporary workaround is
needed. Instead of:
> GridCoverage coverage = CoverageIO.read("~/myfile.grb");
we will need something more elaborated:
NetcdfImageReader netcdfReader = new NetcdfImageReader(null);
netcdfReader.setInput("~/myfile.grb");
ImageCoverageReader reader = new ImageCoverageReader();
reader.setInput(netcdfReader);
GridCoverage2D coverage = reader.read(0, null);
reader.dispose();
netcdfReader.dispose();
The above code should produce the same GridCoverage2D than the one produced by
the convenience CoverageIO method. Once we verified that it work, you can apply
the workaround described in the "/How do I declare the actual georeferencing of
an image file?/" FAQ on
http://www.geotoolkit.org/modules/coverage/faq.html#specifyCRS (you can ignore
the discussion about the widget - unless you want to see your metadata - and
ignore all the code example below "/// If the 'grid to CRS' transform need to be
modified.../").
At this point, you are almost done. Replace the current code:
> values = gc2D.evaluate(new Point2D.Float(600, -420), values);
by:
DirectPosition pos = new DirectPosition2D(DefaultGeographicCRS.WGS84, longitude, latitude);
values = gc2D.evaluate(pos, values);
However before doing so, it may be worth to check if your ProjectedCRS is the
expected one. The magnitude of the (600, -420) coordinates make me suspect that
they are kilometres rather than metres, in which case we need to instruct the
CRS so. First, I suggest that you try the following:
MathTransform tr = CRS.findMathTransform(crs, DefaultGeographicCRS.WGS84);
System.out.println(tr.transform(pos));
where 'pos' is the above DirectPosition object. See if the result seems
reasonable. If not because the linear units should be kilometres, we would need
to change the construction of the CartesianCS object.
Regards,
Martin
More information about the Geotoolkit
mailing list